unit Fil_XYZ_1u;

interface{-------------------------------------------------------------}
  const
    Nx=40;  {net nodes amount on X axe}
    Ny=40;   {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=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_current1=5e-8; {A} {7.871606067e-5;     {Current density,A/mm^2}
    El_current2=1.6e-9;{A} {7.871606067e-5;     {Current density,A/mm^2}
    S=2e-5;          {m*m}
    l=5e-3;          {m}
    dz=0.5;        {step,mm}

    gamma0=0;
    gamma1=1e-4;      {--Skin--}
    gamma2=0.1;         {--Subcutaneous fat--}
    gamma3=1;         {--m--}
    gamma4=1e-7;      {--k--}
    gamma5=1e-6;      {--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-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))<196))
          then gamma[i,j,k]:=gamma2;
      if ((((i-21)*(i-21)+(j-21)*(j-21))<=144{144})and(k<=133)or
         (((i-21)*(i-21)+(j-21)*(j-21)+(k-133)*(k-133))<144{144}))
         then gamma[i,j,k]:=gamma3;
      if ((((i-21)*(i-21)+(j-21)*(j-21))<=64{64})and(k<=133)or
         (((i-21)*(i-21)+(j-21)*(j-21)+(k-133)*(k-133))<64{64}))
         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-21)*(i-21)+(j-21)*(j-21))<=225){and(k>127130)and(k<=144){and(i>10)and
         (i<30)}and(i<20)and(((j-21)*(j-21)+(k-137)*(k-137))<={64}36)
         and(((i-21)*(i-21)+(j-21)*(j-21))>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-21)*(i-21)+(j-21)*(j-21))>225)and(k>40{40})and(k<60{60})and(((i-21)*(i-21)+(j-21)*(j-21))<=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-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 gamma[i,j,k]:=gamma10;



    end;

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


procedure FreeBoundary;

var
  UG, G : single;
  i, j, k : 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,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}0.45;
     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)}-0.87;

     end;




      {-----Electrode Surface-----
  G:=0;
  UG:=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 G:=G+gamma1; end;

  upredict:=(UG-El_current2*1*dz)/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))<=256))
   then u1[i,j,k]:=upredict; end;
   writeln(upredict);



   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 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-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 G:=G+gamma1; end;

  upredict:=(-UG+El_current1*1*dz)/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 u1[i,j,k]:=upredict;
    end;
   writeln(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, S, Ss  : single;
  i, j, k : integer;
const ll=1; l=1;

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


  Cd:=0;  Ss:=0;
   for i:=1 to Mr do
      for j:=1 to Mz do
      for k:=1 to Mk do begin
    C:=0; S:=0;
   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 begin
         Ji:=gamma1{[i+ll,j+ll,k]}*(u[i+ll,j+ll,k]-u[i,j,k])/(ll*0.5e-3);
       if Ji > 0 then
         C:=C+Ji;    S:=S+1;
    end;
    Cd:=Cd+C;   Ss:=Ss+S;
  end;
  Cur:=Cd*pi*sqr(5e-3)/(4*Ss);
  write(Cd,Cur,1/Cur);
  readln;

  Cd:=0; Ss:=0;
   for i:=1 to Mr do
      for j:=1 to Mz do
      for k:=1 to Mk do begin
    C:=0; S:=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 begin
         if (i<=21)and(j<=21) then
         Ji:=gamma1*(-u[i+l,j+l,k]+u[i,j,k])/(l*0.5e-3)
         else if (i>21)and(j<=21) then
         Ji:=gamma1*(-u[i-l,j+l,k]+u[i,j,k])/(l*0.5e-3)
         else if (i<=21)and(j>21) then
         Ji:=gamma1*(-u[i+l,j-l,k]+u[i,j,k])/(l*0.5e-3)
         else if (i>21)and(j>21) then
         Ji:=gamma1*(-u[i-l,j-l,k]+u[i,j,k])/(l*0.5e-3);
       if Ji > 0 then
          C:=C+Ji;    S:=S+1;
    end;
    Cd:=Cd+C;   Ss:=Ss+S;
  end;
  Cur:=Cd*pi*{75e-4}2e-2*1e-2/Ss;
  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:=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,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;

    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.