unit Fil_XYZ_otlad_kate;

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


    gamma1=1e-5;      {--rogov.cloi--}
    gamma2=0.1;       {--zern.sloi--}
    gamma3=5;         {--shipov.sloi--}
    gamma4=0.01;      {--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--}
   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 (k>8)and(k<12) then gamma[i,j,k]:=gamma2;//зернистый слой//
        if (k>=12)and(k<19) then gamma[i,j,k]:=gamma3;//шиповатый слой//
         if (k>=19)and(k<25) then gamma[i,j,k]:=gamma4;   //базальный слой//
          if (k>=25)and(k<37) then gamma[i,j,k]:=gamma5;  //сосочковый слой//
            if (k>=37)and(k<55) then gamma[i,j,k]:=gamma6;    //сетчатый слой//
              if (k>=55)and(k<60) then gamma[i,j,k]:=gamma7; //жировая клетчатка//
             //сосуды,расположенные в сосочковом слое//

         if ((((i-15)*(i-15)+(j-15)*(j-15))<=100)and((k>=28)and(k<31)))or
         ((j<5)and((k>=27)and(k<28)))or((j>25)and((k>=29)and(k<30)))
          then gamma[i,j,k]:=gamma9;

                                      //  нервные клетки  //

          if ((((i-20)*(i-20)+(j-21)*(j-21))<=1)and(k=55))or
         ((((i-30)*(i-30)+(j-38)*(j-38))<=1)and(k=55))or
         ((((i-5)*(i-5)+(j-32)*(j-32))<=1)and(k=55))or
         ((((i-42)*(i-42)+(j-24)*(j-24))<=1)and(k=55))or
         ((((i-55)*(i-55)+(j-40)*(j-40))<=1)and(k=55))
          then gamma[i,j,k]:=gamma10;
          //                //
            if ((((i-21)*(i-21)+(j-51)*(j-51)+(k-26)*(k-26))<=1))or
         ((((i-25)*(i-25)+(j-31)*(j-31)+(k-26)*(k-26))<=1))or
         ((((i-31)*(i-31)+(j-29)*(j-29)+(k-26)*(k-26))<=1))or
         ((((i-39)*(i-39)+(j-37)*(j-37)+(k-26)*(k-26))<=1))or
         ((((i-45)*(i-45)+(j-31)*(j-31)+(k-26)*(k-26))<=1))
          then gamma[i,j,k]:=gamma10;
      if ((((i-7)*(i-7)+(j-25)*(j-25))<=9)and((k>=50)and(k<58)))or
         ((((i-27)*(i-27)+(j-16)*(j-16))<=9)and((k>=50)and(k<58)))or
         ((((i-35)*(i-35)+(j-41)*(j-41))<=9)and((k>=50)and(k<58)))or
         ((((i-44)*(i-44)+(j-19)*(j-19))<=9)and((k>=50)and(k<58)))
          then gamma[i,j,k]:=gamma10;
    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.7;
        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 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;
   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:=gamma[i+1,j+1,k]*(u[i+3,j+3,k]-u[i,j,k])/(3*0.5e-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;

  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))>=225)and(k>40)and(k<60)and(((i-21)*(i-21)+(j-21)*(j-21))<=256))
   then begin
         if (i<=21)and(j<=21) then
         Ji:=gamma1*(u[i+3,j+3,k]-u[i,j,k])/(3*0.5e-3)
         else if (i>21)and(j<=21) then
         Ji:=gamma1*(u[i-3,j+3,k]-u[i,j,k])/(3*0.5e-3)
         else if (i<=21)and(j>21) then
         Ji:=gamma1*(u[i+3,j-3,k]-u[i,j,k])/(3*0.5e-3)
         else if (i>21)and(j>21) then
         Ji:=gamma1*(u[i-3,j-3,k]-u[i,j,k])/(3*0.5e-3);
       if Ji < 0 then
         C:=C+Ji;
    end;
    Cd:=Cd+C;
  end;

  Cur:=Cd*pi*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<356{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.