unit Maxv;

interface{-------------------------------------------------------------}
  const
    Nx=250;  {net nodes amount on X axe}
    Ny=100;   {net nodes amount on Y axe}
    Nt=10;
    HH=4e-4;    {rectangle hight }
    WW=1;    {rectangle width, usually is taken 1}
    TT=6.3e-3;{1e-3;    {dt=1e-3 s}
    Mr=Nx+1;
    Mz=Ny+1;
    sdStop=3.4e-4;     {stop value for relative deviation}
    rlc=13;          {real output digits amount}
    STau=0.9 {0.01};        {value for relaxation parameter}
    El_potential=1;  {Volts}
    El_current=0.5;{7.871606067e-5;     {Current density,A/sm^2}
    dz=0.005;        {step,sm}
    W=0.6;         {W=1000 Hz}
    G=1/2;         {1/4 < G < 1/2}

    hm=0.1 * 0.04;
    tm=6.3e-2;

    S0=0;
    S1=1e-5;      {--Epidermis--}
    S2=0.1;      {--Pupillary dermis--Subcutaneous fat--}
    S3=0.01;      {--Reticular dermis--}
    S4=5;      {--Corpulent cells--}
    S5=1;      {--pora-Protein--}
    S6=0.5;      {--Nervous cells--}
    S7=1;      {m}
    S8=1e-7;      {k}

    E0=18.85e-12;
    E1=E0 * 15;{25;{3;       {--Epidermis--}
    E2=E0 * 35;{10;      {--Pupillary dermis--Subcutaneous fat--}
    E3=E0 * 30;{20;{10;      {--Reticular dermis--}
    E4=E0 * 5;{2;{5;       {--Corpulent cells--}
    E5=E0 * 20;{5;{20;      {--pora-Protein--}
    E6=E0 * 25;{7;{25;      {--Nervous cells--}
    E7=E0 * 57;{5;{30;      {m}
    E8=E0 * 21;{30;{2;       {k}

    M1=1.26e-6;

    Start_A=50;
    Last_A=60;
    Start_P=175;
    Last_P=200;

  type
    U_Matrix = array[1..Mr,1..Mz,0..Nt] of single;
    GammaArray = array[1..Mr,1..Mz] of single;

  var
    tau : double;      {relaxation parameter}
    {i, j, n,} num : integer;
    u, u1 : U_Matrix;

     {auxillary parameters while computing}
    sr, sz, tau1, upredict, hr, hz : double;
    S, E{, M} : GammaArray;

     {deviation sum control parameters}
    sd, sdp, sab, sad, sadMax, sabMax : double;

  procedure SetParameters;
  procedure Compute;
  procedure Curr;

implementation{--------------------------------------------------------}


procedure SetParameters;

var

  i, j : word;

begin  {----------------------Set parameters----------------------}
  tau:=STau;

   {auxillary parameters while computing:}
  tau1:=1-tau;
  hr:=HH/Nx;
  hz:=WW/Ny;




    {deviation sum control parameters}
  sd:=1e+9;
  sadMax:=Mr*Mz; sabMax:=2*(Mz+Mr);

  for i:=1 to Mr do
    for j:=1 to Mz do

    begin {----------------------gamma----------------------}

      S[i,j]:=S7;
      E[i,j]:=E7;
      if (j<15) then begin
        S[i,j]:=S0;
        E[i,j]:=E0;
      end;
      if (j>45)and(j<=70) then begin
        S[i,j]:=S8;
        E[i,j]:=E8;
      end;
      if (j=15)or(j=100) then begin
        S[i,j]:=S1;
        E[i,j]:=E1;
      end;
      if ((j>=16)and(j<=17))or((j>=98)and(j<=99))or((j>20)and(j<=25))or((j>=90)and(j<95)) then begin
      S[i,j]:=S2;
      E[i,j]:=E2;
      end;
      if ((j>17)and(j<=20))or((j>=95)and(j<98)) then begin
        S[i,j]:=S3;
        E[i,j]:=E3;
      end;
      if ((j>=27)and(j<34))or((j>=80)and(j<87))and ((i>45)and(i<65)and(j>=17)and(j<=18)) then begin
        S[i,j]:=S4;
        E[i,j]:=E4;
      end;
      if ((j>=36)and(j<=39))and ((i>54)and(i<=56)and(j>18)and(j<=21)and((sqr(i-55)+sqr(j-20))<4))then begin
        S[i,j]:=S6;
        E[i,j]:=E6;
      end;

   {   S[i,j]:=S2;
      E[i,j]:=E2;
      M[i,j]:=M1;
      if j=1 then begin
        S[i,j]:=S1;
        E[i,j]:=E1;
        M[i,j]:=M1;
      end;
      if (j=1)and((i=5)or(i=10)or(i=15)) then begin
        S[i,j]:=S5;
        E[i,j]:=E5;
        M[i,j]:=M1;
      end;
      if (j>4)and(j<=9) then begin
        S[i,j]:=S3;
        E[i,j]:=E3;
        M[i,j]:=M1;
      end;
      if (j=3)or(j=8) then begin
        S[i,j]:=S4;
        E[i,j]:=E4;
        M[i,j]:=M1;
      end;
      if (j=6)or
         ((i>6)and(i<=8)and(j>7)and(j<=10)and((sqr(i-7)+sqr(j-9))<4))or
         ((i>6)and(i<=8)and(j>=10)and(j<=12)and((sqr(i-7)+sqr(j-11))<4)or
         (i>=6)and(i<=9)and(j>=9)and(j<=11)) then begin
        S[i,j]:=S6;
        E[i,j]:=E6;
        M[i,j]:=M1;
      end;  }

    end;

end;{----------------------Set parameters----------------------}


procedure FreeBoundary;

var
  {UG, G : single;}
  i, j, n : integer;

begin {-----------------------FreeBoundary----------------------}

   {-----Right Boundary-----}
    for n:=0 to Nt do
   { begin
    {begin
     { for i:=1 to Mr do begin
       u1[i,1,n]:=0;
       u1[i,Mz,n]:=0;
       if (i>=Start_A)and(i<=Last_A) then
        u1[i,24,n]:= El_Potential - sin( w * n );
       if (i>=Start_P)and(i<=Last_P) then
        u1[i,24,n]:=0;
      end; }

      for i:=1 to Mr do
      for j:=1 to Mz do begin
      
        u1[1,j,n+1]:=0;
        u1[Mr,j,n+1]:=0;
        u1[i,1,n+1]:=0;
        u1[i,Mz,n+1]:=0;
      if (i>=Start_A)and(i<=Last_A)and(j=15) then
      u1[i,j,n+1]:= - El_Potential * sin( w * (n+1) ) else
     if (i>=Start_P)and(i<=Last_P)and(j=15) then
      u1[i,j,n+1]:= El_Potential * sin( w * (n+1) ); { else
  {   for j:=1 to Mz do
      begin
        u1[1,j,n+1]:=0;
        u1[Mr,j,n+1]:=0;
      end;
      for j:=1 to Mr do
      begin
        u1[j,1,n+1]:=0;
        u1[j,Mz,n+1]:=0;
      end;   }
    end;

    {-----Left Boundary-----
  for n:=0 to Nt do
  for j:=1 to Mz-1 do
  begin
    u1[1,j,n+1]:=0; u1[Mr,j,n+1]:=0;

  end;

   {-----Down Boundary-----
  for n:=0 to Nt do
  for i:=1 to Mr do
   { if (i>Start_A)and(i<Last_A){and(i<Start_P))or(i>Last_P) then

    begin

    end;

   {-----Upper Boundary-----
  for i:=2 to Mr-1 do
 {  if (i<Start_P)or(i>Last_P){and(i<Start_P))or(i>Last_P) then

  begin

  end;



      {-----Electrode Surface-----

   for n:=0 to Nt do
    for i:=1 to Mr do
    if (i>Start_P)and(i<Last_P) then

    begin

       u1[i,1,n+1]:= El_Potential - sin( w * (n+1) ); { end

      else    begin

      u1[i,1,n]:=0;
    end;      }

end;{------------------------FreeBoundary----------------------}


Procedure Curr;

const
 Tab=chr(9);

var
  Ji, C, Cd, Cur : single;
  k, l, m : integer;
  f, uc, r : text;

begin{--------------------Current Computing-------------------}
  assign(f,'Cur4.prn');
    rewrite(f);
   { readln;}
    assign(uc,'UCur2.prn');
    rewrite(uc);
     assign(r,'R4.prn');
    rewrite(r);
 {   readln;}
    for m:=1 to Nt do
  begin
     {  Cur:=0;}
    

     Cd:=0;
     for k:=Start_A to Last_A do
     begin
       C:=0;
       for l:=16 to 20 do
       begin
        Ji:=S1*(u[k,l-1,m]-u[k,l,m])/hh;
        if Ji > 0 then
         C:=C+Ji;
       end;
       Cd:=Cd+C;
     end;
    Cur:=Cd*pi*sqr(4e-3)/4;
    if Cur<>0 then
    writeln(m,u[55,15,m],Cur,u[55,15,m]/Cur);
    
    write(f,Cur:16);
   { write(f,tab);} writeln(f,'');
    write(uc,u[55,15,m]:16);
    writeln(uc,'');
     if Cur<>0 then    begin
    write(r,u[55,15,m]/Cur:16);
    writeln(r,'');     end;
   end;
   readln;


     close (uc);
  close (f);
end;{---------------------Current Computing--------------------}


procedure Compute;

var
  IterNum, i, j, n : integer;
 { Ut : U_Matrix;}

begin

  IterNum:=0;
    FreeBoundary;

 {---------------------computation itself---------------------}

  while iterNum<1000{sd>sdStop} do
  begin {-----main loop-----}

    {sdp:=sd;}  sad:=0; sab:=0;

     for n:=1 to Nt do

     for i:=2 to Mr-1 do  {-----inner points scanning-----}
    { if (i>=Start_A)and(i<=Last_A) then
      u1[i,15,n+1]:= - El_Potential * sin( w * n ) else
     if (i>=Start_P)and(i<=Last_P) then
      u1[i,15,n+1]:= El_Potential * sin( w * n )  else}
    for j:=2 to Mz-1 do {begin
     if (j<>15)and((i<Start_A)or((i>Last_A)and(i<Start_P))or(i>Last_P)) then}

      begin

    {  if (i>=Start_A)and(i<=Last_A)and(j=15) then
      u1[i,j,n]:= - El_Potential {* sin( w * n ) else
     if (i>=Start_P)and(i<=Last_P)and(j=15) then
      u1[i,j,n]:= El_Potential {* sin( w * n )  else  }

       if (S[i,j]<>S[i-1,j])and(S[i,j]<>S[i,j-1]) then
         upredict:=( S[i-1,j-1] * U[i-1,j-1,n+1] + S[i,j] * U[i+1,j+1,n+1] ) /
                   ( S[i-1,j-1] + S[i,j] ) else

       if (S[i,j]<>S[i-1,j])and(S[i,j]<>S[i,j+1]) then
         upredict:=( S[i-1,j+1] * U[i-1,j+1,n+1] + S[i,j] * U[i+1,j-1,n+1] ) /
                   ( S[i-1,j+1] + S[i,j] ) else

       if (S[i,j]<>S[i+1,j])and(S[i,j]<>S[i,j-1]) then
         upredict:=( S[i,j] * U[i-1,j+1,n+1] + S[i+1,j-1] * U[i+1,j-1,n+1] ) /
                   ( S[i+1,j-1] + S[i,j] ) else

       if (S[i,j]<>S[i+1,j])and(S[i,j]<>S[i,j+1]) then
         upredict:=( S[i,j] * U[i-1,j-1,n+1] + S[i+1,j+1] * U[i+1,j+1,n+1] ) /
                   ( S[i+1,j+1] + S[i,j] ) else


       if S[i,j]<>S[i-1,j] then
         upredict:=( S[i-1,j] * U[i-1,j,n+1] + S[i,j] * U[i+1,j,n+1] ) /
                   ( S[i-1,j] + S[i,j] ) else
       if S[i+1,j]<>S[i,j] then
         upredict:=( S[i,j] * U[i-1,j,n+1] + S[i+1,j] * U[i+1,j,n+1] ) /
                   ( S[i+1,j] + S[i,j] ) else

       if S[i,j]<>S[i,j-1] then
         upredict:=( S[i,j-1] * U[i,j-1,n+1] + S[i,j] * U[i,j+1,n+1] ) /
                   ( S[i,j-1] + S[i,j] ) else
       if S[i,j+1]<>S[i,j] then
         upredict:=( S[i,j] * U[i,j-1,n+1] + S[i,j+1] * U[i,j+1,n+1] ) /
                   ( S[i,j+1] + S[i,j] )

       else

        upredict:= {- ( U[i+1,j,n] + U[i-1,j,n] + U[i,j+1,n] + U[i,j-1,n] ) /
                   ( w * M1 * HH*HH * ( E[i,j] - S[i,j] ) -4 ); }
       { ( U[i,j,n] * ( 8 * G - 4 + ( E[i,j] - S[i,j] * sqrt(+1)) *
                   M1 * w * hh*hh * hm )) - U[i,j,n-1] * 4 * G +
                   ( U[i+1,j,n+1] + U[i-1,j,n+1] + U[i,j+1,n+1] + U[i,j-1,n+1] +
                   U[i+1,j,n-1] + U[i-1,j,n-1] + U[i,j+1,n-1] + U[i,j-1,n-1] - 2 *
                   ( U[i+1,j,n] + U[i-1,j,n] + U[i,j+1,n] + U[i,j-1,n] ) * G +
                   U[i+1,j,n] + U[i-1,j,n] + U[i,j+1,n] + U[i,j-1,n] ) / ( 4 * G );

  {( U[i-1,j-1,n] + U[i+1,j+1,n] + U[i-1,j+1,n] + U[i+1,j-1,n] +
                   4 * ( U[i+1,j,n] + U[i-1,j,n] + U[i,j+1,n] + U[i,j-1,n] )) / 20;

                  {( U[i,j,n] * 4 * ( M1 * E[i,j] * HH*HH - 2 * TT*TT ) +
                    U[i,j,n-1] * M1 * HH*HH * ( S[i,j] * TT - 2 * E[i,j] )+
                    2 * TT*TT * ( U[i+1,j,n] + U[i,j+1,n] + U[i,j-1,n] + U[i-1,j,n] ))
                    /( HH*HH * M1 * ( S[i,j] * TT + 2 * E[i,j] ));

                   (( U[i-1,j,n] + U[i+1,j,n] + U[i,j-1,n] + U[i,j+1,n] ) *
                   2 * tt*tt * ( 1 - 2 * G ) + U[i,j,n] * ( 8 * tt*tt * ( 2 * G - 1 )
                   + 4 * M1 * E[i,j] * hh*hh ) + U[i,j,n-1] *
                   ( M1 * hh*hh * ( S[i,j] * tt - 2 * E[i,j] )
                   - 8 * G * tt*tt ) + ( U[i+1,j,n+1] + U[i-1,j,n+1] +
                   U[i,j+1,n+1] + U[i,j-1,n+1] + U[i+1,j,n-1] + U[i-1,j,n-1] +
                   U[i,j+1,n-1] + U[i,j-1,n-1] ) * 2 * G * tt*tt ) /
                   ( 8 * G * tt*tt + ( S[i,j] * tt + 2 * E[i,j] ) * M1 * hh*hh );}


                   (( U[i-1,j,n] + U[i+1,j,n] + U[i,j-1,n] + U[i,j+1,n] ) *
                   2 * tt*tt * tm*tm * ( 1 - 2 * G ) + U[i,j,n] * ( 8 * tt*tt * tm*tm * ( 2 * G - 1 )
                   + 4 * M1 * E[i,j] * hh*hh * hm ) + U[i,j,n-1] *
                   ( M1 * hh*hh * hm * ( S[i,j] * tt * tm - 2 * E[i,j] )
                   - 8 * G * tt*tt * tm*tm ) + ( U[i+1,j,n+1] + U[i-1,j,n+1] +
                   U[i,j+1,n+1] + U[i,j-1,n+1] + U[i+1,j,n-1] + U[i-1,j,n-1] +
                   U[i,j+1,n-1] + U[i,j-1,n-1] ) * 2 * G * tt*tt *tm*tm ) /
                   ( 8 * G * tt*tt * tm*tm + ( S[i,j] * tt * tm + 2 * E[i,j] ) * M1 * hh*hh * hm );

                 {  (( U[i-1,j,n] + U[i+1,j,n] + U[i,j-1,n] + U[i,j+1,n] ) *
                   2 * tt*tt * hm * ( 1 - 2 * G ) + U[i,j,n] * ( 8 * tt*tt * hm * ( 2 * G - 1 )
                   + 4 * M1 * E[i,j] * hh*hh * tm*tm ) + U[i,j,n-1] *
                   ( M1 * hh*hh * tm * ( S[i,j] * tt - 2 * tm * E[i,j] )
                   - 8 * G * tt*tt * hm ) + ( U[i+1,j,n+1] + U[i-1,j,n+1] +
                   U[i,j+1,n+1] + U[i,j-1,n+1] + U[i+1,j,n-1] + U[i-1,j,n-1] +
                   U[i,j+1,n-1] + U[i,j-1,n-1] ) * 2 * G * tt*tt * hm ) /
                   ( 8 * G * tt*tt * hm + ( S[i,j] * tt + 2 * tm * E[i,j] ) * M1 * hh*hh * tm ); }


        u1[i,j,n+1]:= {u[i,j,n+1] + tau * upredict;}
                      tau * upredict + (1-tau) * u[i,j,n+1];

        sad:=sad+abs(upredict-u[i,j,n+1]);
      end;

    FreeBoundary;

  {  for n:=1 to Nt-2 do
    for i:=1 to Mr do
    if (i>Start_P)and(i<Last_P) then

    begin
{      if S[i,1]<>S[i-1,1] then
        upredict:=( S[i-1,1] * U[i-1,1] + S[i,1] * U[i+1,1] ) /
                  ( S[i-1,1] + S[i,1] )
    {  if S[i+1,1]<>S[i,1] then
        upredict:=( S[i,1] * U[i-1,1] + S[i+1,1] * U[i+1,1] ) /
                  ( S[i+1,1] + S[i,1] )
       u1[i,1,n+1]:= El_Potential {- sin( w * n+1 ); { end

      else    begin

      u1[i,1,n]:=0;
    end;   }

    sd:=sad/sadMax+sab/sabMax;
    IterNum:=IterNum+1;

    u:=u1;

    writeln('sd=',sd:rlc,'  r.sad=',sad/sadMax:rlc,'  r.sab=',{sab/sabMax}u[55,30,10]:rlc,'  iter=',iternum);

  end; {-----main loop-----}

end;{---------------------------Compute--------------------------- }


Begin
End.