unit Fil_XYZ_2;

interface{-------------------------------------------------------------}
  const
    Nx=50;  {net nodes amount on X axe}
    Ny=50;   {net nodes amount on Y axe}
    Nk=50;
    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=10;  {Volts}
    El_current=0.5;{7.871606067e-5;     {Current density,A/sm^2}
    dz=0.5;        {step,sm}

    gamma0=0;
    gamma1=1e-7;      {--Epidermis--}
    gamma2=1e-1;      {--Pupillary dermis--Subcutaneous fat--}
    gamma3=1e-2;      {--Reticular dermis--}
    gamma4=2;      {--Corpulent cells--}
    gamma5=1e-2;      {--pora-Protein--}
    gamma6=5e-3;      {--Nervous cells--}


    gamma7=1;         {--m--}
    gamma8=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----------------------}

      gamma[i,j,k]:=gamma2;
      if j<=2 then gamma[i,j,k]:=gamma1;
      if (j>6)and(j<=12) then gamma[i,j,k]:=gamma3;
      if (j>17)and(j<=27) then gamma[i,j,k]:=gamma7;
      if (j>27)and(j<=Mz) then gamma[i,j,k]:=gamma8;

      if (j>20)and(j<=23) then gamma[i,j,k]:=gamma4;
      if (j=6)then gamma[i,j,k]:=gamma4;

      if (j=10) then gamma[i,j,k]:=gamma6;
      if (j>25)and(j<=27) then gamma[i,j,k]:=gamma6;
      if (((i-15)*(i-15)+(j-12)*(j-12)+(k-35)*(k-35))<=7) then gamma[i,j,k]:=gamma6;
      if (((i-15)*(i-15)+(j-12)*(j-12)+(k-15)*(k-15))<=7) then gamma[i,j,k]:=gamma6;
      if (((i-35)*(i-35)+(j-12)*(j-12)+(k-35)*(k-35))<=7) then gamma[i,j,k]:=gamma6;
      if (((i-35)*(i-35)+(j-12)*(j-12)+(k-15)*(k-15))<=7) then gamma[i,j,k]:=gamma6;
      if (i=15)and(k=35)and(j>=12)and(j<27) then gamma[i,j,k]:=gamma6;
      if (i=15)and(k=15)and(j>=12)and(j<27) then gamma[i,j,k]:=gamma6;
      if (i=35)and(k=35)and(j>=12)and(j<27) then gamma[i,j,k]:=gamma6;
      if (i=35)and(k=15)and(j>=12)and(j<27) then gamma[i,j,k]:=gamma6;
      if (((i-12)*(i-12)+(j-4)*(j-4)+(k-12)*(k-12))<=2) then gamma[i,j,k]:=gamma6;
      if (((i-24)*(i-24)+(j-4)*(j-4)+(k-24)*(k-24))<=2) then gamma[i,j,k]:=gamma6;
      if (((i-36)*(i-36)+(j-4)*(j-4)+(k-36)*(k-36))<=2) then gamma[i,j,k]:=gamma6;
      if (((i-12)*(i-12)+(j-4)*(j-4)+(k-24)*(k-24))<=2) then gamma[i,j,k]:=gamma6;
      if (((i-24)*(i-24)+(j-4)*(j-4)+(k-12)*(k-12))<=2) then gamma[i,j,k]:=gamma6;
      if (((i-36)*(i-36)+(j-4)*(j-4)+(k-12)*(k-12))<=2) then gamma[i,j,k]:=gamma6;
      if (((i-12)*(i-12)+(j-4)*(j-4)+(k-36)*(k-36))<=2) then gamma[i,j,k]:=gamma6;
      if (((i-24)*(i-24)+(j-4)*(j-4)+(k-36)*(k-36))<=2) then gamma[i,j,k]:=gamma6;
      if (((i-36)*(i-36)+(j-4)*(j-4)+(k-24)*(k-24))<=2) then gamma[i,j,k]:=gamma6;
      if (i=12)and(k=12)and(j>=4)and(j<11) then gamma[i,j,k]:=gamma6;
      if (i=12)and(k=24)and(j>=4)and(j<11) then gamma[i,j,k]:=gamma6;
      if (i=12)and(k=36)and(j>=4)and(j<11) then gamma[i,j,k]:=gamma6;
      if (i=24)and(k=12)and(j>=4)and(j<11) then gamma[i,j,k]:=gamma6;
      if (i=24)and(k=24)and(j>=4)and(j<11) then gamma[i,j,k]:=gamma6;
      if (i=24)and(k=36)and(j>=4)and(j<11) then gamma[i,j,k]:=gamma6;
      if (i=36)and(k=12)and(j>=4)and(j<11) then gamma[i,j,k]:=gamma6;
      if (i=36)and(k=24)and(j>=4)and(j<11) then gamma[i,j,k]:=gamma6;
      if (i=36)and(k=36)and(j>=4)and(j<11) then gamma[i,j,k]:=gamma6;

      if (((i-25)*(i-25)+(j-15)*(j-15)+(k-25)*(k-25))<=7) then gamma[i,j,k]:=gamma5;
      if (((i-25)*(i-25)+(j-15)*(j-15)+(k-10)*(k-10))<=7) then gamma[i,j,k]:=gamma5;
      if (((i-25)*(i-25)+(j-15)*(j-15)+(k-38)*(k-38))<=7) then gamma[i,j,k]:=gamma5;
      if (((i-10)*(i-10)+(j-15)*(j-15)+(k-25)*(k-25))<=7) then gamma[i,j,k]:=gamma5;
      if (((i-38)*(i-38)+(j-15)*(j-15)+(k-25)*(k-25))<=7) then gamma[i,j,k]:=gamma5;
      if (i>=25)and(i<=26)and(k>=25)and(k<=26)and(j<=15) then gamma[i,j,k]:=gamma5;
      if (i>=25)and(i<=26)and(k>=10)and(k<=11)and(j<=15) then gamma[i,j,k]:=gamma5;
      if (i>=25)and(i<=26)and(k>=37)and(k<=38)and(j<=15) then gamma[i,j,k]:=gamma5;
      if (i>=10)and(i<=11)and(k>=25)and(k<=26)and(j<=15) then gamma[i,j,k]:=gamma5;
      if (i>=37)and(i<=38)and(k>=25)and(k<=26)and(j<=15) then gamma[i,j,k]:=gamma5;
      
    end;

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


procedure FreeBoundary;

var
  UG, G : single;
  i, j, k : integer;

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

  {-----Right Boundary-----}
  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]);

    upredict:=(U[i+1,j,Mk] + U[i-1,j,Mk] + U[i,j+1,Mk] + U[i,j-1,Mk] + U[i,j,Mk-1])/5;
    u1[i,j,Mk]:=u[i,j,Mk]*(1-tau)+tau*upredict;
    sab:=sab+abs(upredict-u[i,j,Mk]);

  end;


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

{   if ((i-9)*(i-9)+(k-9)*(k-9))>100
      ((i<5)or(i>13))and((k<5)or(k>13)) then }

   begin
  {  upredict:=(U[i+1,1,k] + U[i-1,1,k] + U[i,2,k] + U[i,1,k+1] + U[i,1,k-1])/5;}
    u1[i,1,k]:=-1{u[i,1,k]*(1-tau)+tau*upredict};
{    sab:=sab+abs(upredict-u[i,1,k]);            }

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

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

    upredict:=(U[Mr-1,j,k] + U[Mr,j+1,k] + U[Mr,j-1,k] + U[Mr,j,k+1] + U[Mr,j,k-1])/5;
    u1[Mr,j,k]:=u[Mr,j,k]*(1-tau)+tau*upredict;
    sab:=sab+abs(upredict-u[Mr,j,k]);
  end;





   {-----Left Boundary-----}
  for i:=2 to Mr-1 do
  begin
  {  upredict:=( U[i+1,1,1] + U[i-1,1,1] + U[i,2,1] + U[i,1,2] ) / 4;}
    u1[i,1,1]:=-1{u[i,1,1]*(1-tau)+tau*upredict;
    sab:=sab+abs(upredict-u[i,1,1])};

{    upredict:=( U[i+1,1,Mk] + U[i-1,1,Mk] + U[i,2,Mk] + U[i,1,Mk-1] ) / 4;}
    u1[i,1,Mk]:=-1{u[i,1,Mk]*(1-tau)+tau*upredict;
    sab:=sab+abs(upredict-u[i,1,Mk])};

    upredict:=( U[i+1,Mz,1] + U[i-1,Mz,1] + U[i,Mz-1,1] + U[i,Mz,2] ) / 4;
    u1[i,Mz,1]:=u[i,Mz,1]*(1-tau)+tau*upredict;
    sab:=sab+abs(upredict-u[i,Mz,1]);

    upredict:=( U[i+1,Mz,Mk] + U[i-1,Mz,Mk] + U[i,Mz-1,Mk] + U[i,Mz,Mk-1] ) / 4;
    u1[i,Mz,Mk]:=u[i,Mz,Mk]*(1-tau)+tau*upredict;
    sab:=sab+abs(upredict-u[i,Mz,Mk]);
  end;


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

    upredict:=( U[Mr-1,j,1] + U[Mr,j+1,1] + U[Mr,j-1,1] + U[Mr,j,2] ) / 4;
    u1[Mr,j,1]:=u[Mr,j,1]*(1-tau)+tau*upredict;
    sab:=sab+abs(upredict-u[Mr,j,1]);

    upredict:=( U[2,j,Mk] + U[1,j+1,Mk] + U[1,j-1,Mk] + U[1,j,Mk-1] ) / 4;
    u1[1,j,Mk]:=u[1,j,Mk]*(1-tau)+tau*upredict;
    sab:=sab+abs(upredict-u[1,j,Mk]);

    upredict:=( U[Mr-1,j,Mk] + U[Mr,j+1,Mk] + U[Mr,j-1,Mk] + U[Mr,j,Mk-1] ) / 4;
    u1[Mr,j,Mk]:=u[Mr,j,Mk]*(1-tau)+tau*upredict;
    sab:=sab+abs(upredict-u[Mr,j,Mk]);

  end;

   for k:=2 to Mk-1 do
  begin
   { upredict:=( U[2,1,k] + U[1,2,k] + U[1,1,k+1] + U[1,1,k-1] ) / 4;}
    u1[1,1,k]:=-1{u[1,1,k]*(1-tau)+tau*upredict;
    sab:=sab+abs(upredict-u[1,1,k])};

    upredict:=( U[2,Mz,k] + U[1,Mz-1,k] + U[1,Mz,k+1] + U[1,Mz,k-1] ) / 4;
    u1[1,Mz,k]:=u[1,Mz,k]*(1-tau)+tau*upredict;
    sab:=sab+abs(upredict-u[1,Mz,k]);

{    upredict:= ( U[Mr-1,1,k] + U[Mr,2,k] + U[Mr,1,k+1] + U[Mr,1,k-1] ) / 4;}
    u1[Mr,1,k]:=-1{u[Mr,1,k]*(1-tau)+tau*upredict;
    sab:=sab+abs(upredict-u[Mr,1,k])};

    upredict:= ( U[Mr-1,Mz,k] + U[Mr,Mz-1,k] + U[Mr,Mz,k+1] + U[Mr,Mz,k-1] )/4;
    u1[Mr,Mz,k]:=u[Mr,Mz,k]*(1-tau)+tau*upredict;
    sab:=sab+abs(upredict-u[Mr,Mz,k]);

  end;



     {-----left down corner----}
{  upredict:=( U[2,1,1] + U[1,2,1] + U[1,1,2] ) / 3;}
  u1[1,1,1]:=-1{u[1,1,1]*(1-tau)+tau*upredict;
  sab:=sab+abs(upredict-u[1,1,1])};

{  upredict:=( U[2,1,Mk] + U[1,2,Mk] + U[1,1,Mk-1] ) / 3;}
  u1[1,1,Mk]:=-1{u[1,1,Mk]*(1-tau)+tau*upredict;
  sab:=sab+abs(upredict-u[1,1,Mk])};

  upredict:=( U[2,Mz,1] + U[1,Mz-1,1] + U[1,Mz,2] ) / 3;
  u1[1,Mz,1]:=u[1,Mz,1]*(1-tau)+tau*upredict;
  sab:=sab+abs(upredict-u[1,Mz,1]);

{  upredict:=( U[Mr-1,1,1] + U[Mr,2,1] + U[Mr,1,2] ) / 3;}
  u1[Mr,1,1]:=-1{u[Mr,1,1]*(1-tau)+tau*upredict;
  sab:=sab+abs(upredict-u[Mr,1,1])};

  upredict:=( U[2,Mz,Mk] + U[1,Mz-1,Mk] + U[1,Mz,Mk-1] ) / 3;
  u1[1,Mz,Mk]:=u[1,Mz,Mk]*(1-tau)+tau*upredict;
  sab:=sab+abs(upredict-u[1,Mz,Mk]);

{  upredict:=( U[Mr-1,1,Mk] + U[Mr,2,Mk] + U[Mr,1,Mk-1] ) / 3;}
  u1[Mr,1,Mk]:=-1{u[Mr,1,Mk]*(1-tau)+tau*upredict;
  sab:=sab+abs(upredict-u[Mr,1,Mk])};

  upredict:=( U[Mr-1,Mz,Mk] + U[Mr,Mz-1,Mk] + U[Mr,Mz,Mk-1] ) / 3;
  u1[Mr,Mz,Mk]:=u[Mr,Mz,Mk]*(1-tau)+tau*upredict;
  sab:=sab+abs(upredict-u[Mr,Mz,Mk]);

  upredict:=( U[Mr-1,Mz,1] + U[Mr,Mz-1,1] + U[Mr,Mz,2] ) / 3;
  u1[Mr,Mz,1]:=u[Mr,Mz,1]*(1-tau)+tau*upredict;
  sab:=sab+abs(upredict-u[Mr,Mz,1]);

 {1 G:=0;
  UG:=0;

   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-9)*(i-9)+(k-9)*(k-9))<100
        (i>21)and(i<40)and(k>21)and(k<40) then begin
     {u1[i,5,k]:=20;     u1[i,Mz-6,k]:=-20;
     UG:=UG+U[i,6,k]*gamma[i,5,k];
     G:=G+gamma[i,5,k];
      end;
     end;
  upredict:=(UG+El_current*dz)/G;
      for i:=1 to Mr do
      for j:=1 to Mz do
      for k:=1 to Mk do begin

     if (i>21)and(i<40)and(k>21)and(k<40) then  begin
     u1[i,5,k]:=-upredict;
     u1[i,Mz-6,k]:=upredict;
     end;
     end; 1}

  {2  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-9)*(i-9)+(k-9)*(k-9))<100
        (i>21)and(i<40)and(k>21)and(k<40)and((j=5)or(j=Mz-6)) then begin
      upr
      u1[i,5,k]:=-(El_current*dz/(gamma[i,5,k]*20*20*dz)+U[i,6,k]);
      u1[i,Mz-6,k]:=0{(El_current*dz/(gamma[i,Mz-6,k]*20*20*dz)+U[i,Mz-7,k]);
     end;
     end;  2}

 {      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]:=-1;
        u1[i,Mz,k]:=0;
        u1[i,j,1]:=0;
        u1[i,j,Mk]:=0;
  {   if (i>21)and(i<40)and(k>21)and(k<40)and((j=2)or(j=Mz-3))
      then begin
      u1[i,2,k]:=-5;
      u1[i,Mz-3,k]:=0;
     end;
     end;


     {-----Electrode Surface-----
  G:=0;
  UG:=0;

  for i:=Start_A to Last_A do
    UG:=UG+U[i,2]*gamma[i,1];
  for i:=Start_A to Last_A do
    G:=G+gamma[i,1];

  upredict:=(UG-El_current*dz)/G;
  for i:=Start_A to Last_A do
    u1[i,1]:=upredict;


 { 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, C, Cd, Cur : single;

begin{--------------------Current Computing-------------------
  write(u1[40,1]);
  readln;
  Cd:=0;
  for i:=Start_A to Last_A do
  begin
    C:=0;
    for j:=1 to 1 do
    begin
       Ji:=gamma[i,j]*(u[i,j+1]-u[i,j])/0.005;
       if Ji > 0 then
         C:=C+Ji;
    end;
    Cd:=Cd+C;
  end;
  Cur:=Cd*pi*sqr(4e-2)/4;
  write(Cd,Cur,5/Cur);
  readln;
end;{---------------------Current Computing--------------------}


Procedure Current;

var
  Ji, 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;
   begin
         Ji:=gamma[i,1,k]*(u[i,4,k]-u[i,1,k])/(3*0.1e-3);
       if Ji > 0 then
         C:=C+Ji;
    end;
    Cd:=Cd+C;
  end;
  Cur:=Cd*pi*sqr(5e-3)/4;
  write(Cd,Cur,1/Cur);
  readln;
end;{---------------------Current Computing--------------------}


procedure Compute;

var
  IterNum, i, j, k : integer;

begin
  IterNum:=0;

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

  while IterNum<1500{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-----}

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


Begin
End.