unit Fil_XYZ_5;

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

    gamma0={1e-6}0;
    gamma1=5e-7{5};      {--Skin--}
    gamma2=0.1{0001};         {--Subcutaneous fat--}
    gamma3=1{0.00000};         {--m--}
    gamma4=1e-7{10};      {--k--}
    gamma5=1e-7;      {--n--}
    gamma6=1e-2;      {--pora-Protein--}
    gamma7=5e-3;      {--Nervous cells--}
    gamma10=6e7;      {electrode}

    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}

    Start_A=32;
    Last_A=53;
    Start_P=100;
    Last_P=150;

  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;

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


procedure SetParameters;

var
  {hr2, hz2, }sh : double;
  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;

 { hz2:=sqr(hz);
  hr2:=sqr(hr);}
  sh:=1/(sqr(hr)+sqr(hz));
  sr:=sqr(hz)*sh;
  sz:=sqr(hr)*sh;

    {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
      for k:=1 to Mk do

    begin {----------------------gamma----------------------}
      {amma[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){or
         (((i-20)*(i-20)+(j-20)*(j-20)+(k-130)*(k-130))>=196)   )
          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))<169))
          then gamma[i,j,k]:=gamma2;
      if ((((i-21)*(i-21)+(j-21)*(j-21))<=100)and(k<=133)or
         (((i-21)*(i-21)+(j-21)*(j-21)+(k-133)*(k-133))<100))
         then gamma[i,j,k]:=gamma3;
      if ((((i-21)*(i-21)+(j-21)*(j-21))<=225){and(k>127130)and(k<=144){and(i>10)and
         (i<30)        and(j<20)and(((i-21)*(i-21)+(k-137)*(k-137))<=36)
         and(((i-21)*(i-21)+(j-21)*(j-21))>196))
         then gamma[i,j,k]:=gamma4;   }

       gamma[i,j,k]:=gamma0;
      if ((((i-22)*(i-22)+(j-22)*(j-22))<=400{225})and(k<=140{133})or
         ((((i-22)*(i-22)+(j-22)*(j-22)+(k-136{40}{133})*(k-136{40}{133}))<=400{225})and(k>140)){or
         (((i-20)*(i-20)+(j-20)*(j-20)+(k-130)*(k-130))>=196)})
          then gamma[i,j,k]:=gamma1;
      if ((((i-22)*(i-22)+(j-22)*(j-22))<289{196})and(k<=140{133})or
         (((i-22)*(i-22)+(j-22)*(j-22)+(k-{140}136)*(k-{140}136))<289{196})and(k>140))
          then gamma[i,j,k]:=gamma2;
      if ((((i-22)*(i-22)+(j-22)*(j-22))<{=}225{144})and(k<=140{133})or
         (((i-22)*(i-22)+(j-22)*(j-22)+(k-{140}136)*(k-{140}136))<225{144})and(k>140))
         then gamma[i,j,k]:=gamma3;
      if ((((i-22)*(i-22)+(j-22)*(j-22))<{=}100{64})and(k<=140{133})or
         (((i-22)*(i-22)+(j-22)*(j-22)+(k-{140}136)*(k-{140}136{137{133}))<100{64})and(k>140))
         then gamma[i,j,k]:=gamma4;
      {if ((((i-21)*(i-21)+(j-21)*(j-21))<=225){and(k>127130)and(k<=144){and(i>10)and
         (i<30)and(j<20)and(((i-21)*(i-21)+(k-137)*(k-137))<={6436)
         and(((i-21)*(i-21)+(j-21)*(j-21))>169{144{196)) }
      if ((((i-22)*(i-22)+(j-22)*(j-22))<=400{225}){and(k>127130)and(k<=144){and(i>10)and
         (i<30)}and(i<20)and(((j-22)*(j-22)+(k-140{137})*(k-140{137}))<={64}36)
         and(((i-22)*(i-22)+(j-22)*(j-22))>324{169}{144{196}))
         then gamma[i,j,k]:=gamma5;

{      if ((((i-21)*(i-21)+(j-21)*(j-21))>=225)and(k>60)and(k<80)and(((i-21)*(i-21)+(j-21)*(j-21))<=361)) 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))<=361))}
       if ((((i-22)*(i-22)+(j-22)*(j-22))>400{225})and(k>40{40})and(k<60{60}){and(((i-22)*(i-22)+(j-22)*(j-22))<=484{400{361)}) or
         { ((((i-21)*(i-21)+(j-21)*(j-21))>=225)and(k>115)and(k<=125)and(i>=5{7)and(i<10{17)and(j>=18)and(j<=22{20)and(((i-21)*(i-21)+(j-21)*(j-21))<=400{361))  }
           ((((i-22)*(i-22)+(j-22)*(j-22))>400{225})and(k>115{3})and(k<=125{7})and(j>=17{5})and(j<27{9})and(i<20){and(((i-22)*(i-22)+(j-22)*(j-22))<=484{400{256)})

    {     ((((i-21)*(i-21)+(j-21)*(j-21))>225)and(k>115)and(k<=125)and{(i>=7)and(i<17)and(j<12)and(((i-21)*(i-21)+(j-21)*(j-21))<=256))}

         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:=1 to Mr do
      for j:=1 to Mz do
      for k:=1 to Mk do begin
      if ((((i-22)*(i-22)+(j-22)*(j-22))>529)and(k>141)or
         (((i-22)*(i-22)+(j-22)*(j-22)+(k-140)*(k-140))>529)) then

        u1[i,j,k]:=0;  end;}


        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,1]:=0;}
        u1[i,j,Mk]:=0;

   

 { 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
     u1[i,j,k]:={0...1;
     if {      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))<=361))
((((i-21)*(i-21)+(j-21)*(j-21))>=225)and(k>115)and(k<=125)and(i>=5{7)and(i<10{17)and(j>=18)and(j<=22{20)and(((i-21)*(i-21)+(j-21)*(j-21))<=256))...
 ((((i-20)*(i-20)+(j-20)*(j-20))>=225)and(k>115)and(k<=125)and(i>=7)and(i<17)and(j<20)and(((i-20)*(i-20)+(j-20)*(j-20))<=256))
        then
     u1[i,j,k]:={-El_Current*l/(S*gamma10)-1;   }

     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))>=400{225)and(k>40{60)and(k<60{80)and(((i-21)*(i-21)+(j-21)*(j-21))<=441{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-22)*(i-22)+(j-22)*(j-22))>=400{225})and(k>40{60})and(k<60{80})and(((i-22)*(i-22)+(j-22)*(j-22))<=441{256}))
   then begin  u1[i,j,k]:=0.52;
       {  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-22)*(i-22)+(j-22)*(j-22))>=400{225})and(k>115)and(k<=125)and(j>=17{5})and(j<27{9})and(i<20)and(((i-22)*(i-22)+(j-22)*(j-22))<=441{256}))
   then u1[i,j,k]:=-0.8;
  { 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))>=400{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))<=441{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);


{  G:=0;
  UG:=0;
  G2:=0;
  Ug2:=0;

   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{60...)and(k<60{80...)and(((i-21)*(i-21)+(j-21)*(j-21))<=256))
   then begin
         if (i<=21)and(j<=21) then
         UG:=UG+U[i+1,j+1,k]*gamma1
         else if (i>21)and(j<=21) then
         UG:=UG+U[i-1,j+1,k]*gamma1
         else if (i<=21)and(j>21) then
         UG:=UG+U[i+1,j-1,k]*gamma1
         else if (i>21)and(j>21) then
         UG:=UG+U[i-1,j-1,k]*gamma1; 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>40{60...)and(k<60{80...)and(((i-21)*(i-21)+(j-21)*(j-21))<=256))
   then begin G2:=G2+1; G:=G+gamma1; end; end;
       upredict:=(UG-El_current2*1*dz{*G2...)/G;
   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{60...)and(k<60{80...)and(((i-21)*(i-21)+(j-21)*(j-21))<=400{256...))
   then begin u1[i,j,k]:=u[i,j,k]*(1-tau)+tau*upredict;
    sab:=sab+abs(upredict-u[i,j,k]);  end;
     end;
   writeln(u1[35,35,50]{upredict...);

  G:=0;
  UG:=0;
  G2:=0;

   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-20)*(i-20)+(j-20)*(j-20))<=256))
   then UG:=UG+U[i+1,j+1,k]*gamma1; 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-20)*(i-20)+(j-20)*(j-20))<=256))
   then begin G2:=G2+1; G:=G+gamma1; end; end;

  upredict:=(-UG+El_current1*1*dz{*G2...)/Gamma1;
   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-20)*(i-20)+(j-20)*(j-20))<=400{256...))
   then begin u1[i,j,k]:=u[i,j,k]*(1-tau)+tau*upredict;
    sab:=sab+abs(upredict-u[i,j,k]); end;
    end;
   writeln(u1[12,12,120]{upredict...);  }





 {   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{60)and(k<60{80)and(((i-21)*(i-21)+(j-21)*(j-21))<=256))
   then begin
         if (i<=21)and(j<=21) then
         U1[i,j,k]:=(gamma1 * U[i+3,j+3,k] - El_current2*3*dz) / gamma1
         else if (i>21)and(j<=21) then
         U1[i,j,k]:=(gamma1 * U[i-3,j+3,k] - El_current2*3*dz) / gamma1
         else if (i<=21)and(j>21) then
         U1[i,j,k]:=(gamma1 * U[i+3,j-3,k] - El_current2*3*dz) / gamma1
         else if (i>21)and(j>21) then
         U1[i,j,k]:=(gamma1 * U[i-3,j-3,k] - El_current2*3*dz) / gamma1;
    Ug2:=Ug2+U1[i,j,k];
    G2:=G2+1;
    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>115{3)and(k<=125{7)and(i>=7{5)and(i<17{9)and(j<20)and(((i-21)*(i-21)+(j-21)*(j-21))<=256))

  {    ((((i-21)*(i-21)+(j-21)*(j-21))>=225)and(k>115)and(k<=125)and{(i>=7)and(i<17)and(j<12)and(((i-21)*(i-21)+(j-21)*(j-21))<=256))...
    then begin U1[i,j,k]:=(gamma1{[i+1,j+1,k+1]... * U[i+3,j+3,k] - El_current1*3*dz) / gamma1{[i,j,k]...;
    Ug:=Ug+U1[i,j,k];
    G:=G+1;
    end; end;

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

 {  for i:=1 to Mr do
      for j:=1 to Mz do
      for k:=1 to Mk do begin
   if ((((i-20)*(i-20)+(j-20)*(j-20))>=225)and(k>115)and(k<=125)and(i>=7)and(i<17)and(j<20)and(((i-20)*(i-20)+(j-20)*(j-20))<=256))
   then  write(u1[i,j,k], ','); end;    }




 { for i:=Start_A to Last_A do u1[i,1]:=-El_potential;
  for i:=Start_P to Last_P do u1[i,1]:=El_Potential; }

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


Procedure Current;

var
  Ji, Ji2, C, 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{ begin
    C:=0; }
   if ((((i-21)*(i-21)+(j-21)*(j-21))>=400{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))<=441{256}))
   then begin
         Ji:=gamma1*(u[i+1,j+1,k]-u[i,j,k])/(1*dz);
      { if Ji > 0 then}
       {  C:=C+Ji;
    end; }
    Cd:=Cd+Ji{C};
  end;
  Cur:=Cd*pi*sqr(5e-3)/4;
  write(Cd,Cur{,1/Cur});
  readln;

  Cd:=0;
   for i:=1 to Mr do
      for j:=1 to Mz do
      for k:=1 to Mk do{ begin
    C:=0;   }
   if
   ((((i-21)*(i-21)+(j-21)*(j-21))>=40{225})and(k>40{60})and(k<60{80})and(((i-21)*(i-21)+(j-21)*(j-21))<=441{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);
     {  if Ji2 > 0 then}
     {    C:=C+Ji;
    end; }
    Cd:=Cd+Ji2{C};
  end;
  Cur:=Cd*pi*2e-2*1e-2{75e-4*1e-2};
  write(Cd,Cur{,1/Cur});
  readln;
end;{---------------------Current Computing--------------------}


procedure Compute;

var
  IterNum, i, j, k : integer;

begin
  IterNum:=0;

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

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

    {sdp:=sd;}  sad:=0; sab:=0;
    for i:=1 to Mr do
      for j:=1 to Mz do
      for k:=1 to Mk do begin
      if ((((i-22)*(i-22)+(j-22)*(j-22))>529{225})and(k>141{133})or
         (((i-22)*(i-22)+(j-22)*(j-22)+(k-{140}136)*(k-{140}136))>529{225})) then

        u1[i,j,k]:=0;  end;
        {FreeBoundary;}

    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 ((((i-22)*(i-22)+(j-22)*(j-22))<=400{529}{225})and(k<=141{133})or
         (((i-22)*(i-22)+(j-22)*(j-22)+(k-{140}136)*(k-{140}136))<=400{529}{225})) then 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,k-1]) then
         upredict:=( gamma[i-1,j,k-1] * U[i,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+1,j,k])and(gamma[i,j,k]<>gamma[i,j+1,k]) then
         upredict:=( gamma[i,j,k] * U[i-1,j-1,k] + gamma[i+1,j+1,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,j,k] * U[i,j,k-1] + gamma[i+1,j,k+1] * 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,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] )    }

       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;      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.