unit LIZA;

interface{-------------------------------------------------------------}
  const
    Nx=45;  {net nodes amount on X axe}
    Ny=45;   {net nodes amount on Y axe}
    Nk=150;
    HH=1;    {rectangle hight }
    WW=1;    {rectangle width, usually is taken 1}
    KK=1;
    Mr=Nx+1;
    Mz=Ny+1;
    Mk=Nk+1;
    sdStop={5e-4}3e-3;     {stop value for relative deviation}
    rlc=13;          {real output digits amount}
    STau=0.9;        {value for relaxation parameter}
    El_potential=0.1;  {Volts}
    El_current1=5e-2;   {Current density,A/m^2}
    El_current2=1.6e-3; {Current density,A/m^2}
    S=2e-5;          {m*m}
    l=5e-3;          {m}
    dz=0.5e-3;        {step,m}

    gamma0=0;
    gamma1=1e-7;      {--Skin--}
    gamma2=0.1;         {--Subcutaneous fat--}
    gamma3=1;         {--m--}
    gamma4=1e-7;      {--k--}
    gamma5=1e-6;      {--n--}
    gamma10=6e7;      {electrode}


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

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

     {auxillary parameters while computing}
    sr, sz, tau1, upredict, hr, hz, hk : double;
    gamma : GammaArray;

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

  procedure SetParameters;
  procedure Compute;
  procedure Current;
  Procedure Emkost;
implementation{--------------------------------------------------------}


procedure SetParameters;

var
   i, j, k : word;

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

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

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

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

    begin {----------------------gamma----------------------}
        gamma[i,j,k]:=gamma0;
      if ((((i-21)*(i-21)+(j-21)*(j-21))<=225)and(k<=133)or
         (((i-21)*(i-21)+(j-21)*(j-21)+(k-133)*(k-133))<=225))
          then gamma[i,j,k]:=gamma1;
      if ((((i-21)*(i-21)+(j-21)*(j-21))<196)and(k<=133)or
         (((i-21)*(i-21)+(j-21)*(j-21)+(k-133)*(k-133))<196))
          then gamma[i,j,k]:=gamma2;
      if ((((i-21)*(i-21)+(j-21)*(j-21))<=144)and(k<=133)or
         (((i-21)*(i-21)+(j-21)*(j-21)+(k-133)*(k-133))<144))
         then gamma[i,j,k]:=gamma3;
      if ((((i-21)*(i-21)+(j-21)*(j-21))<=64)and(k<=133)or
         (((i-21)*(i-21)+(j-21)*(j-21)+(k-133)*(k-133))<64))
         then gamma[i,j,k]:=gamma4;
      if ((((i-21)*(i-21)+(j-21)*(j-21))<=225)and(i<20)and(((j-21)*(j-21)+(k-137)*(k-137))<=36)
         and(((i-21)*(i-21)+(j-21)*(j-21))>169))
         then gamma[i,j,k]:=gamma5;

       if ((((i-21)*(i-21)+(j-21)*(j-21))>225)and(k>40)and(k<60)and(((i-21)*(i-21)+(j-21)*(j-21))<=400)) or
          ((((i-21)*(i-21)+(j-21)*(j-21))>225)and(k>115)and(k<=125)and(i>=7)and(i<17)and(j<20)and(((i-21)*(i-21)+(j-21)*(j-21))<=400))
         then gamma[i,j,k]:=gamma10;
       end;

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


procedure FreeBoundary;

var
  UG, G, Ug2, G2, Ji1, ji2 : single;
  i, j, k, p, q : integer;

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


    for i:=2 to Mr-1 do
   for j:=2 to Mz-1 do
   begin
    upredict:=(U[i+1,j,1] + U[i-1,j,1] + U[i,j+1,1] + U[i,j-1,1] + U[i,j,2])/5;
    u1[i,j,1]:=u[i,j,1]*(1-tau)+tau*upredict;
    sab:=sab+abs(upredict-u[i,j,1]);
   end;


   for i:=1 to Mr do
      for j:=1 to Mz do
      for k:=1 to Mk do begin

        u1[1,j,k]:=0;
        u1[Mr,j,k]:=0;
        u1[i,1,k]:=0;
        u1[i,Mz,k]:=0;
        u1[i,j,Mk]:=0;

     end;

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

   G:=0;
   G2:=0;
   Ug:=0;
   Ug2:=0;
   p:=1;
   q:=1;


    for i:=1 to Mr do
      for j:=1 to Mz do
      for k:=1 to Mk do begin
    if ((((i-21)*(i-21)+(j-21)*(j-21))>=225)and(k>40)and(k<60)and(((i-21)*(i-21)+(j-21)*(j-21))<=256))
   then
    G2:=G2+1;
    end;
    Ji2:=El_current2/G2;

    for i:=1 to Mr do
      for j:=1 to Mz do
      for k:=1 to Mk do
    if ((((i-21)*(i-21)+(j-21)*(j-21))>=225)and(k>40{60})and(k<60{80})and(((i-21)*(i-21)+(j-21)*(j-21))<=256))
   then begin
         if (i<=21)and(j<=21) then begin
         Upredict:=(gamma1 * U[i+p,j+p,k] - Ji2*p*dz) / gamma1;
         u1[i,j,k]:=u[i,j,k]*(1-tau)+tau*upredict;
         sab:=sab+abs(upredict-u[i,j,k]);
         end
         else if (i>21)and(j<=21) then begin
         Upredict:=(gamma1 * U[i-p,j+p,k] - Ji2*p*dz) / gamma1;
         u1[i,j,k]:=u[i,j,k]*(1-tau)+tau*upredict;
         sab:=sab+abs(upredict-u[i,j,k]);
         end
         else if (i<=21)and(j>21) then begin
         Upredict:=(gamma1 * U[i+p,j-p,k] - Ji2*p*dz) / gamma1;
         u1[i,j,k]:=u[i,j,k]*(1-tau)+tau*upredict;
         sab:=sab+abs(upredict-u[i,j,k]);
         end
         else if (i>21)and(j>21) then begin
         Upredict:=(gamma1 * U[i-p,j-p,k] - Ji2*p*dz) / gamma1;
         u1[i,j,k]:=u[i,j,k]*(1-tau)+tau*upredict;
         sab:=sab+abs(upredict-u[i,j,k]);
         end;
         Ug2:=Ug2+U1[i,j,k];
         end;


  for i:=1 to Mr do
      for j:=1 to Mz do
      for k:=1 to Mk do begin
   if ((((i-21)*(i-21)+(j-21)*(j-21))>=225)and(k>115)and(k<=125)and(i>=7)and(i<17)and(j<20)and(((i-21)*(i-21)+(j-21)*(j-21))<=256))
   then
   G:=G+1;
    end;
    Ji1:=El_current1/G;
  for i:=1 to Mr do
      for j:=1 to Mz do
      for k:=1 to Mk do
   if ((((i-21)*(i-21)+(j-21)*(j-21))>=225)and(k>115)and(k<=125)and(i>=7)and(i<17)and(j<20)and(((i-21)*(i-21)+(j-21)*(j-21))<=256))
   then begin
    Upredict:=(gamma1* U[i+q,j+q,k] - Ji1*q*dz) / gamma1;
    u1[i,j,k]:=u[i,j,k]*(1-tau)+tau*upredict;
         sab:=sab+abs(upredict-u[i,j,k]);
         Ug:=Ug+U1[i,j,k];
         end;

  writeln (Ug/G,' ',Ug2/G2,' ',Ug2/G2-Ug/G);


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


Procedure Current;

var
  Ji, Ji2, Cd, Cur : single;
  i, j, k : integer;

begin{--------------------Current Computing-------------------}


  Cd:=0;
   for i:=1 to Mr do
      for j:=1 to Mz do
      for k:=1 to Mk do
   if ((((i-21)*(i-21)+(j-21)*(j-21))>=225)and(k>115)and(k<=125)and(i>=7)and(i<17)and(j<21)and(((i-21)*(i-21)+(j-21)*(j-21))<=256))
   then begin
         Ji:=gamma1*(u[i+1,j+1,k]-u[i,j,k])/(1*dz);
        Cd:=Cd+Ji;
  end;
  Cur:=Cd*pi*sqr(5e-3)/4;
  write(Cd,', cur1=',Cur);
  readln;

  Cd:=0;
   for i:=1 to Mr do
      for j:=1 to Mz do
      for k:=1 to Mk do
   if
   ((((i-21)*(i-21)+(j-21)*(j-21))>=225)and(k>40)and(k<60)and(((i-21)*(i-21)+(j-21)*(j-21))<=256))
   then begin
         if (i<=21)and(j<=21) then
         Ji2:=gamma1*(u[i+1,j+1,k]-u[i,j,k])/(1*dz)
         else if (i>21)and(j<=21) then
         Ji2:=gamma1*(u[i-1,j+1,k]-u[i,j,k])/(1*dz)
         else if (i<=21)and(j>21) then
         Ji2:=gamma1*(u[i+1,j-1,k]-u[i,j,k])/(1*dz)
         else if (i>21)and(j>21) then
         Ji2:=gamma1*(u[i-1,j-1,k]-u[i,j,k])/(1*dz);
      Cd:=Cd+Ji2;
  end;
  Cur:=Cd*pi*2e-2*1e-2;
  write(Cd,', cur2=',Cur);
  readln;
end;{---------------------Current Computing--------------------}

procedure Emkost;       {---------------------Emkost Computing--------------------}

var
   Ei,Ej,Ek,E,Esi,Esj,Esk,C,Q:real;
   i,j,k:integer;
begin
  Ei:=0;
  Ej:=0;
  Ek:=0;
  E:=0;
  Esi:=0;
  Esj:=0;
  Esk:=0;
  for j:=1 to Mz do
      for k:=1 to Mk do
      for i:=1 to Mr do
       begin
         if ((((i-21)*(i-21)+(j-21)*(j-21))<=225)and(k<=133)or
         (((i-21)*(i-21)+(j-21)*(j-21)+(k-133)*(k-133))<=225)and(k>133))
          then
           begin
            Esi:=(U[i+1,j,k]-U[i,j,k])/dz;
            Ei:=Ei+Esi;
           end;
       end;
  for i:=1 to Mr do
      for k:=1 to Mk do
      for j:=1 to Mz do
       begin
         if ((((i-21)*(i-21)+(j-21)*(j-21))<=225)and(k<=133)or
         (((i-21)*(i-21)+(j-21)*(j-21)+(k-133)*(k-133))<=225)and(k>133))
          then
           begin
            Esj:=(U[i,j+1,k]-U[i,j,k])/dz;
            Ej:=Ej+Esj;
           end;
       end;
  for i:=1 to Mr do
      for j:=1 to Mz do
      for k:=1 to Mk do
       begin
         if ((((i-21)*(i-21)+(j-21)*(j-21))<=225)and(k<=133)or
         (((i-21)*(i-21)+(j-21)*(j-21)+(k-133)*(k-133))<=225)and(k>133))
          then
           begin
            Esk:=(U[i,j,k+1]-U[i,j,k])/dz;
            Ek:=Ek+Esk;
           end;
       end;
  E:=sqrt(Ei*Ei+Ej*Ej+Ek*Ek);
  Q:=E*0.00361632;
  C:=Q/1.32;
  writeln(' q=',Q,' C=',C);
  readln;
end;
procedure Compute;

var
  IterNum, i, j, k : integer;

begin
  IterNum:=0;

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

  while sd>sdStop do
  begin {-----main loop-----}

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

    for i:=2 to Mr-1 do  {-----inner points scanning-----}
      for j:=2 to Mz-1 do
        for k:=2 to Mk-1 do

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

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

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

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

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

       else

        upredict:= ( U[i+1,j,k] + U[i-1,j,k] + U[i,j+1,k] + U[i,j-1,k] +
                   U[i,j,k+1] + U[i,j,k-1] ) / 6;



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

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

    FreeBoundary;

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

    u:=u1;

    writeln('sd=',sd:rlc,{'  r.sad=',sad/sadMax:rlc,'  r.sab=',sab/sabMax:rlc,}'  iter=',iternum);

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

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


Begin
End.