unit Fil_XYZ_MY;

interface{-------------------------------------------------------------}
  const
    Nx=50;   {net nodes amount on X axe}
    Ny=50;   {net nodes amount on Y axe}
    Nk=50;
    Mx=Nx+1;
    My=Ny+1;
    Mk=Nk+1;
    sdStop=2e-6;     {stop value for relative deviation}
    rlc=9;           {real output digits amount}
    STau=0.9;        {value for relaxation parameter}


    gamma1=1e-5;      {--epidermis--}
    gamma2=0.1;       {--sosochk.sloi--}
    gamma3=5;         {--vaskular--}
    gamma4=0.01;      {--setchat.sloi--}
    gamma5=5e-3;       {--Nervous cells--}
    gamma6=0.1;       {--Subcutaneous fat--}
    gamma7=1;         {Muscle}
    gamma8=1e-2;      {--pora-Protein--}

  type
    U_Matrix = array[1..Mx,1..My,1..Mk] of single;
    GammaArray = array[1..Mx,1..My,1..Mk] of single;

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

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

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

  procedure SetParameters;
  procedure Compute;
  {procedure Current; }

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


procedure SetParameters;

var
    i, j, k : word;

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

   {auxillary parameters while computing:}
  {tau1:=1-tau;}

   {deviation sum control parameters}
  sd:=1e+9;
  sadMax:=Mx*My*Mk;
  for k:=1 to Mk do
   for i:=1 to Mx do
    for j:=1 to My do
     begin      {----------------------gamma----------------------}
       gamma[i,j,k]:=gamma1;
      if (k>2)and(k<7) then gamma[i,j,k]:=gamma2;

      if ((((i-25)*(i-25)+(j-25)*(j-25))<=100)and((k>=3)and(k<6)))or((j<15)
      and((k>=4)and(k<5)))or((j>35)and((k>=4)and(k<5)))
           then gamma[i,j,k]:=gamma3;

      if ((k>=12)and(k<14))or((k>=29)and(k<33))
           then gamma[i,j,k]:=gamma3;
      if ((k>=7)and(k<8))or((k>9)and(k<12))or((k>=14)and(k<=15))
           then gamma[i,j,k]:=gamma4;
      if ((k>=8)and(k<=9))or((k>=40)and(k<42))
           then gamma[i,j,k]:=gamma5;
      if (k>15)and(k<25) then gamma[i,j,k]:=gamma6;
      if ((k>=15)and(k<29))or((k>=33)and(k<40))or((k>=42)and(k<Nk))
           then gamma[i,j,k]:=gamma7;

      if ((((i-22)*(i-22)+(j-27)*(j-27))<=1)and(k<18)or
         (((i-22)*(i-22)+(j-27)*(j-27)+(k-20)*(k-20))<=4))or
          ((((i-7)*(i-7)+(j-7)*(j-7))<=1)and(k<18)or
         (((i-7)*(i-7)+(j-7)*(j-7)+(k-20)*(k-20))<=4))or
         ((((i-35)*(i-35)+(j-9)*(j-9))<=1)and(k<18)or
         (((i-35)*(i-35)+(j-9)*(j-9)+(k-20)*(k-20))<=4))or
         ((((i-43)*(i-43)+(j-39)*(j-39))<=1)and(k<18)or
         (((i-43)*(i-43)+(j-39)*(j-39)+(k-20)*(k-20))<=4))
          then gamma[i,j,k]:=gamma8;
      if ((((i-14)*(i-14)+(j-15)*(j-15))<=1)and(k=15))or
         ((((i-25)*(i-25)+(j-33)*(j-33))<=1)and(k=15))or
         ((((i-6)*(i-6)+(j-33)*(j-33))<=1)and(k=15))or
         ((((i-32)*(i-32)+(j-14)*(j-14))<=1)and(k=15))or
         ((((i-45)*(i-45)+(j-29)*(j-29))<=1)and(k=15))
          then gamma[i,j,k]:=gamma5;
      if ((((i-11)*(i-11)+(j-41)*(j-41)+(k-3)*(k-3))<=1))or
         ((((i-15)*(i-15)+(j-21)*(j-21)+(k-3)*(k-3))<=1))or
         ((((i-18)*(i-18)+(j-30)*(j-30)+(k-3)*(k-3))<=1))or
         ((((i-21)*(i-21)+(j-19)*(j-19)+(k-3)*(k-3))<=1))or
         ((((i-29)*(i-29)+(j-27)*(j-27)+(k-3)*(k-3))<=1))or
         ((((i-33)*(i-33)+(j-33)*(j-33)+(k-3)*(k-3))<=1))or
         ((((i-35)*(i-35)+(j-21)*(j-21)+(k-3)*(k-3))<=1))
          then gamma[i,j,k]:=gamma5;
      if ((((i-7)*(i-7)+(j-25)*(j-25))<9)and((k>=10)and(k<20)))or
         ((((i-27)*(i-27)+(j-16)*(j-16))<9)and((k>=10)and(k<20)))or
         ((((i-35)*(i-35)+(j-41)*(j-41))<9)and((k>=12)and(k<22)))or
         ((((i-44)*(i-44)+(j-19)*(j-19))<9)and((k>=11)and(k<21)))
          then gamma[i,j,k]:=gamma5;
     end;

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


procedure FreeBoundary;

var
   i, j, k : integer;
   delta:single;
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 Mx do
    for j:=1 to My do
     begin
      u1[i,j,1]:=-0.8;
      u1[i,j,Mk]:=-0.7;
     end;
   delta:=abs(u1[1,1,1]-u1[1,1,Mk])/Mk;

   for i:=1 to Mx do
    for k:=1 to Mk do
     begin
      u1[i,1,k+1]:=u1[i,1,k]+delta;
      u1[i,My,k+1]:=u1[i,My,k]+delta;
     end;
    for j:=1 to My do
    for k:=1 to Mk do
     begin
      u1[1,j,k+1]:=u1[1,j,k]+delta;
      u1[Mx,j,k+1]:=u1[Mx,j,k]+delta;
     end;


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


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 Mx do  {-----inner points scanning-----}
      for j:=1 to My do
        for k:=1 to Mk 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.