unit Fil_XYZ_o_k_s_d_e;

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


    gamma1=2e-5;      {--rogov.cloi--}
    gamma2=0.1;       {--zern.sloi--}
    gamma3=0.5;       {--shipov.sloi--}
    gamma4=0.9;       {--bazal.sloi--}
    gamma5=5;         {--sosochk.sloi--}
    gamma6=0.1;       {--setchat.sloi--}
    gamma7=1;         {--subcutaneus fat--}
    gamma8=1;         {--geleza--}
    gamma9=5;         {--vascular--}
    gamma10=5;        {--nervous cells--}
    gamma11=1;        {--dielectric--}
    gamma12=0.8;      {--electrode--}
    gamma13=1e-15;    {--air--}
   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}
    upredict: double;
    gamma : GammaArray;

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

  procedure SetParameters;
  procedure Compute;
 { procedure Current; }

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


procedure SetParameters;

var
    i, j, k : word;

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

   {deviation sum control parameters}
  sd:=1e+9;
  sadMax:=Mx*My*Mk;

    for i:=1 to Mx do
      for j:=1 to My do
        for k:=1 to Mk do
    begin {----------------------gamma----------------------}
    gamma[i,j,k]:=gamma1;    //роговый слой//
      if (j<20)and(k>=8)and(k<12) then gamma[i,j,k]:=gamma2;//зернистый слой//

        if (j<20)and(k>=12)and(k<19) then gamma[i,j,k]:=gamma3;//шиповатый слой//
         if (j<20)and(k>=19)and(k<25) then gamma[i,j,k]:=gamma4;   //базальный слой//
          if (j<20)and(k>=25)and(k<37) then gamma[i,j,k]:=gamma5;  //сосочковый слой//
            if (j<20)and(k>=37)and(k<55) then gamma[i,j,k]:=gamma6;    //сетчатый слой//
              if (j<20)and(k>=55)and(k<60) then gamma[i,j,k]:=gamma7; //жировая клетчатка//
               //ступенька//
             {  if (j>=20)and (j<25) then gamma[i,j,k]:=gamma11};// слой диэлектрика//
                {if (j>=25) then gamma[i,j,k]:=gamma12;//диэлектрик//  }
         if (j>=20)and(k>8)and(k<16) then gamma[i,j,k]:=gamma1;    //роговый слой//
      if (j>=20)and(k>=16)and(k<20) then gamma[i,j,k]:=gamma2;//зернистый слой//
       if (j>=20)and(k>=20)and(k<27) then gamma[i,j,k]:=gamma3;//шиповатый слой//
         if (j>=20)and(k>=27)and(k<30) then gamma[i,j,k]:=gamma4;   //базальный слой//
          if (j>=20)and(k>=30)and(k<37) then gamma[i,j,k]:=gamma5;  //сосочковый слой//
            if (j>=20)and(k>=37)and(k<55) then gamma[i,j,k]:=gamma6;    //сетчатый слой//
              if (j>=20)and(k>=55)and(k<60) then gamma[i,j,k]:=gamma7; //жировая клетчатка//


           {  //сосуды,расположенные в сосочковом слое//

        if ((((i-25)*(i-25)+(j-25)*(j-25))<=100)and((k>=26)and(k<36)))or
         ((j<15)and((k>=29)and(k<33)))or((j>33)and((k>=36)and(k<39)))
          then gamma[i,j,k]:=gamma9;   }

         {                             //  нервные клетки  //
         if (j<20)and( k>=45) and (k<=46))or
         (j>=20)and (k>=46)and (k>=47) then gamma[i,j,k]:=gamma10;
                         //в сетчатом слое//
          if ((((i-10)*(i-10)+(j-10)*(j-10))<=1)and(k=41))or
         ((((i-30)*(i-30)+(j-39)*(j-39))<=1)and(k=41))or
         ((((i-5)*(i-5)+(j-32)*(j-32))<=1)and(k=41))or
         ((((i-42)*(i-42)+(j-21)*(j-21))<=1)and(k=42))or
         ((((i-55)*(i-55)+(j-50)*(j-50))<=1)and(k=43))
          then gamma[i,j,k]:=gamma10;
          //            в сосочковом слое    //
            if ((((i-21)*(i-21)+(j-8)*(j-8)+(k-26)*(k-26))<=1))or
         ((((i-25)*(i-25)+(j-18)*(j-18)+(k-26)*(k-26))<=1))or
         ((((i-31)*(i-31)+(j-32)*(j-32)+(k-26)*(k-26))<=1))or
         ((((i-39)*(i-39)+(j-41)*(j-41)+(k-35)*(k-35))<=1))or
         ((((i-45)*(i-45)+(j-48)*(j-48)+(k-36)*(k-36))<=1))
          then gamma[i,j,k]:=gamma10;
      if ((((i-7)*(i-7)+(j-9)*(j-9))<=9)and((k>=50)and(k<59)))or
         ((((i-27)*(i-27)+(j-20)*(j-20))<=9)and((k>=50)and(k<59)))or
         ((((i-35)*(i-35)+(j-39)*(j-39))<=9)and((k>=50)and(k<59)))or
         ((((i-45)*(i-45)+(j-50)*(j-50))<=9)and((k>=50)and(k<59)))
          then gamma[i,j,k]:=gamma10;   }
        {      //поры//
              if ((((i-5)*(i-5)+(j-5)*(j-5))<=1)and(k<48)or
         (((i-5)*(i-5)+(j-5)*(j-5)+(k-50)*(k-50))<=4))or

         ((((i-11)*(i-11)+(j-15)*(j-15))<=1)and(k<50)or
         (((i-11)*(i-11)+(j-15)*(j-15)+(k-52)*(k-52))<=4))or

         ((((i-21)*(i-21)+(j-25)*(j-25))<=1)and(k<48)or
         (((i-21)*(i-21)+(j-25)*(j-25)+(k-50)*(k-50))<=4))or
         
         ((((i-39)*(i-39)+(j-35)*(j-35))<=1)and(k<50)or
         (((i-39)*(i-39)+(j-35)*(j-35)+(k-52)*(k-52))<=4))or

             ((((i-48)*(i-48)+(j-45)*(j-45))<=1)and(k<49)and (k>16)or
         (((i-48)*(i-48)+(j-45)*(j-45)+(k-51)*(k-51))<=4)) or

           ((((i-51)*(i-51)+(j-55)*(j-55))<=1)and(k<48)and (k>16) or
         (((i-51)*(i-51)+(j-55)*(j-55)+(k-50)*(k-50))<=4))
          then gamma[i,j,k]:=gamma8;   }
    end;


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


procedure FreeBoundary;

var
   i, j, k : integer;
   delta:single;
begin {-----------------------FreeBoundary----------------------}

     for i:=1 to Mx do
      for j:=1 to My do
        begin
       {  u1[i,j,1]:=-0.8; }
         u1[i,j,Mk]:=-0.73;
                   end;

       for i:=1 to Mx do
      for j:=25 to My do
        begin
         u1[i,j,8]:=-0.8;
        { u1[i,j,Mk]:=-0.73;}
       end;
      {        for i:=1 to Mx do
      for j:=20 to 25 do
        begin
         gamma[i,j,k]*(U[i+1,j,k]-U[i,j,k])/0.05=0;     end;   }

              for i:=1 to Mx do
              for k:=8 to Mk do
        begin
         gamma1[i,j,k]*(U[i,j+1,k]-U[i,j,k])/0.05=0;     end;
             for i:=1 to Mx do
              for j:=0 to 20 do
        begin
         gamma1[i,j,k]*(U[i+1,j,k]-U[i,j,k])/0.05=0;     end;



        { u1[i,j,Mk]:=-0.73;
       end;                 }
  {delta:=abs(u1[1,1,1]-u1[1,1,Mk])/Nk;
   for i:=1 to Mx do
      for k:=1 to Mk do
        begin
         u1[i,1,k]:=-0.8+delta*(k-1);
         u1[i,My,k]:=-0.8+delta*(k-1);
        end;
   for j:=1 to My do
      for k:=1 to Mk do
        begin
         u1[1,j,k]:=-0.8+delta*(k-1);
         u1[Mx,j,k]:=-0.8+delta*(k-1);
        end;       }

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


procedure Compute;

var
  IterNum, i, j, k : integer;

begin
  IterNum:=0;

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

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

      sad:=0;

    for i:=2 to Mx{-1} do  {-----inner points scanning-----}
      for j:=2 to My{-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;
    IterNum:=IterNum+1;

    u:=u1;

    writeln('sd=',sd:rlc,'  iter=',iternum);

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

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


Begin
End.