Unit newzapoln;
Interface
Uses Math;
Type coeff_array=array [1..4, 1..2] of real;
const   n=31;
        p=0.15;//1000000; //   область рассматриваемого пространства [мкм]
        h=0.005;//1000000;   //шаг сетки [мкм]
        t=round(p/h);  // количество узлов по одной границе [штук]
        rad=h*4; // радиус скругления
        n_rad = 4;
        r1=10000000000000000;//10000000;  //   удельное сопротивление воздуха [кОм*см]
        r2=4500;//10000000;//     удельное сопротивление кожи [kОм*см]
        r3=0.15;//10000000;//     удельное сопротивление ТА [kОм*см]
        //r4=900;//     удельное сопротивление жировой ткани [Ом*см]
        vv=3*h;//     расстояние электрод/кожа [см]
        v=round(vv/h);//расстояние электрод/кожа[узлов]
        dxx=-4*h;
        dx=round(dxx/h);
        dxx2=-9*h;
        dx2=round(dxx2/h);
        zz=12*h; // до перегиба [см]
        z=round(zz/h); // до перегиба [узлов]
        ww=p-vv;//   толщина слоя кожи [см]
        w=round(ww/h);//толщина слоя кожи [узлов]
        uu=4*h; //  диаметр ТА [см]
        u=round(uu/h);//диаметр ТА [узлов]

Procedure InsUpBottom(n,i:Word;h:real; var F:text);
Procedure InsDownBottom(n,i:Word;h:real; var F:text);
Procedure InsLeftBottom(n,i:Word;h:real;var F:text);
Procedure InsRightBottom(n,i:Word;h:real; var F:text);
Procedure InsMiddlePoint(i:Word;h:real; var F:text);
Procedure InsTissueUBPoint(n,i:Word;r1,r2,h:real; var F:text);
Procedure InsAkupUBPoint(n,i:Word;r1,r3,h:real; var F:text);
Procedure InsTissueDBPoint(n,i:Word;r2,r4,h:real; var F:text);
Procedure InsAkupDBPoint(n,i:Word;r3,r4,h:real; var F:text);
Procedure InsTissueAkupPoint(n,i:Word;r2,r3,h:real; var F:Text);
Procedure InsAkupTissuePoint(n,i:Word;r2,r3,h:real; var F:text);
{Procedure AbovePoint(n,i,flag:integer;r1,r3,h:real; var F:text);
Procedure CurvePoint1(n,i,flag:integer;r2,r3,h:real; var F:text);
Procedure CurvePoint2(n,i,flag:integer;r2,r3,h:real; var F:text);
Procedure BelowPoint(n,i,flag:integer;h:real; var F:text);
Procedure InsidePoint(n,i,flag:integer;h:real; var F:text);}
Procedure InsAkupSkrugPoint (n,i:integer;deltay,r1,r3,h:real;var F:text);
Procedure InsSpecMidPoint (n,i:integer;deltax,deltay,r1,r2,h:real;var F1:text);
Procedure InsSpecMidPoint2 (n,i:integer;deltax,r1,r2,h:real;var F1:text);
procedure VertSpecPoint(h,r2,r3:real;i,n:integer; var arr:coeff_array);
procedure HorSpecPoint(h,r2,r3:real;i,n:integer; var arr:coeff_array);
procedure HorSpecPoint2(h,r2,r3,deltax:real;i,n:integer; var arr:coeff_array);
Implementation
Var
 F:text;


  {Верхняя граница}
Procedure InsUpBottom;
  var k:real;
  j:Longint;
  begin
         for j:=1 to (n*n) do
             begin
                  if j=i then
                     k:=1;
                  if (j<>i) then
                     k:=0;
                  writeln(F,k);
             end;
    end;

  {Нижняя граница}
Procedure InsDownBottom;
    var k:real;
    j:Longint;
    begin
         for j:=1 to (n*n) do
             begin
                  if j=i then
                     begin
                       k:=1;
                     end;
                   if (j<>i)  then
                     k:=0;
                  writeln(F,k);
             end;
     end;

  {Левая граница}
Procedure InsLeftBottom;
     var k:real;
     j:Longint;
     begin
          for j:=1 to n*n do
             begin
                  if j=i then
                     k:=1;
                  if (j<>i) then
                     k:=0;
                  writeln(F,k);
             end;
     end;

  {Правая граница}
Procedure InsRightBottom;
     var k:real;
     j:Longint;
     begin
     for j:=1 to (n*n) do
             begin
                  if j=i then
                     k:=1;
                  if (j<>i)  then
                     k:=0;
                  writeln(F,k);
             end;
     end;

  {Средние точки}
Procedure InsMiddlePoint;
var k:real;
j:Longint;
begin
     for j:=1 to (n*n) do
      begin
       if j=i then
        k:=-4/(h*h);
       if j=i+n then
        k:=1/(h*h);
       if j=i-n then
        k:=1/(h*h);
       if j=i+1 then
        k:=1/(h*h);
       if j=i-1 then
        k:=1/(h*h);
       if (j<>i) and (j<>i+n) and (j<>i-n) and (j<>i+1) and (j<>i-1) then
        k:=0;
       writeln(F,k);
      end;
end;

  {Верхняя граница кожи}
Procedure InsTissueUBPoint;
const q=1000;
var k:real;
j:Longint;
begin
   for j:=1 to (n*n) do
    begin
     if j=i then
      k:=(-1/(r2*h)-1/(r1*h));///q;
     if j=i+n then
      k:=(1/(r2*h));///q;
     if j=i-n then
      k:=(1/(r1*h));///q;
     if (j<>i) and (j<>i+n) and (j<>i-n) then
      k:=0;
     writeln(F,k);
    end;
end;

  {Верхняя граница точки}
Procedure InsAkupUBPoint;
   const q=1000;
   var k:real;
   j:Longint;
   begin
   for j:=1 to (n*n) do
   begin
                  if j=i then
                     k:=(-1/(r3*h)-1/(r1*h));///q;
                  if j=i+n+3 then
                     k:=(1/(r3*h));///q;
                  if j=i-n then
                     k:=(1/(r1*h));///q;
                  if (j<>i) and (j<>i+n) and (j<>i-n) then
                     k:=0;
                  writeln(F,k);
             end;

   end;

  {Нижняя граница ткани}
Procedure InsTissueDBPoint;
   const q=1000;
   var zu,k:real;
   j:Longint;
   begin
   for j:=1 to (n*n) do
   begin
                  if j=i then
                     k:=(-1/(r4*h)-1/(r2*h));///q;
                  if j=i+n then
                     k:=(1/(r4*h));///q;
                  if j=i-n then
                     k:=(1/(r2*h));///q;
                  if (j<>i) and (j<>i+n) and (j<>i-n) then
                     k:=0;
                  writeln(F,k);
             end;
   end;

   {Нижняя граница точки}
Procedure InsAkupDBPoint;
   const q=1000;
   var k:real;
   j:Longint;
   begin
   for j:=1 to n*n do
   begin
                  if j=i then
                     k:=(-1/(r4*h)-1/(r3*h));///q;
                  if j=i+n then
                     k:=(1/(r4*h));///q;
                  if j=i-n then
                     k:=(1/(r3*h));///q;
                  if (j<>i) and (j<>i+n) and (j<>i-n) then
                     k:=0;
                  writeln(F,k);
             end;
   end;

{Граница кожа-точка}
Procedure InsTissueAkupPoint;
   const q=1000;
   var k:real;
   j:Longint;
   begin
   for j:=1 to (n*n) do
   begin
                  if j=i then
                     k:=(-1/(r3*h)-1/(r2*h));///q;
                  if j=i+1 then
                     k:=(1/(r3*h));///q;
                  if j=i-1 then
                     k:=(1/(r2*h));///q;
                  if (j<>i) and (j<>i+1) and (j<>i-1) then
                     k:=0;
                  writeln(F,k);
             end;
   end;

{Граница точка-кожа}
Procedure InsAkupTissuePoint;
   const q=1000;
   var k:real;
   j:Longint;
   begin
   for j:=1 to (n*n) do
   begin
                  if j=i then
                     k:=-1/(r2*h)-1/(r3*h);///q;
                  if j=i+1 then
                     k:=(1/(r2*h));///q;
                  if j=i-1 then
                     k:=(1/(r3*h));///q;
                  if (j<>i) and (j<>i+1) and (j<>i-1) then
                     k:=0;
                  writeln(F,k);
             end;
   end;

////////////////////////////////////////////////////////////////////////////////
/////////////////Новые Процедуры

procedure VertSpecPoint(h,r2,r3:real;i,n:integer; var arr:coeff_array);
const q=1000;
var jj,k: integer;
    x1,x2,y1,y2,b: real;
begin
     jj:=(i mod n) - ((n div 2) + 1 - (u div 2) - n_rad);
     if jj>n_rad then
      begin
       jj:=u + 2*n_rad - jj;
       b:=h*sqrt(n_rad*n_rad-jj*jj);
       y1:=(trunc(b/h) + 1)*h-b;
       y2:=b-trunc(b/h)*h;
       x1:=jj*y1*h/b;
       x2:=jj*y2*h/b;
       arr[1,2]:=((1-x1)/h)*r2*sqrt(x2*x2+y2*y2)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
       arr[1,1]:=i;
       arr[2,2]:=(x1/h)*r2*sqrt(x2*x2+y2*y2)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
       arr[2,1]:=i-1;
       arr[3,2]:=(x2/h)*r3*sqrt(x1*x1+y1*y1)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
       arr[3,1]:=i+n+1;
       arr[4,2]:=((1-x2)/h)*r3*sqrt(x1*x1+y1*y1)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
       arr[4,1]:=i+n;
      end
     else
      begin
       b:=h*sqrt(n_rad*n_rad-jj*jj);
       y1:=(trunc(b/h) + 1)*h-b;
       y2:=b-trunc(b/h)*h;
       x1:=jj*y1*h/b;
       x2:=jj*y2*h/b;
       arr[1,2]:=((1-x1)/h)*r2*sqrt(x2*x2+y2*y2)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
       arr[1,1]:=i;
       arr[2,2]:=(x1/h)*r2*sqrt(x2*x2+y2*y2)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
       arr[2,1]:=i+1;
       arr[3,2]:=(x2/h)*r3*sqrt(x1*x1+y1*y1)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
       arr[3,1]:=i+n-1;
       arr[4,2]:=((1-x2)/h)*r3*sqrt(x1*x1+y1*y1)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
       arr[4,1]:=i+n;
      end;
   end;

procedure HorSpecPoint(h,r2,r3:real;i,n:integer; var arr:coeff_array);
const q=1000;
var jj,k: integer;
    x1,x2,y1,y2,b: real;
begin
     jj:=(i mod n) - ((n div 2) + 1 - (u div 2) - n_rad);
     if jj>n_rad then
      begin
       jj:=u + 2*n_rad - jj;
       b:=h*sqrt(n_rad*n_rad-jj*jj);
       if b=0 then b:=h*sqrt(n_rad*n_rad-(jj-1)*(jj-1));
       y1:=(trunc(b/h) + 1)*h-b;
       y2:=b-trunc(b/h)*h;
       x1:=jj*y1*h/b;
       x2:=jj*y2*h/b;
       arr[1,2]:=((1-x1)/h)*r2*sqrt(x2*x2+y2*y2)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
       arr[1,1]:=i-1;
       arr[2,2]:=(x1/h)*r2*sqrt(x2*x2+y2*y2)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
       arr[2,1]:=i-n-1;
       arr[3,2]:=(x2/h)*r3*sqrt(x1*x1+y1*y1)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
       arr[3,1]:=i+n;
       arr[4,2]:=((1-x2)/h)*r3*sqrt(x1*x1+y1*y1)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
       arr[4,1]:=i;
      end
     else
      begin
       b:=h*sqrt(n_rad*n_rad-jj*jj);
       y1:=(trunc(b/h) + 1)*h-b;
       y2:=b-trunc(b/h)*h;
       x1:=jj*y1*h/b;
       x2:=jj*y2*h/b;
       arr[1,2]:=((1-x1)/h)*r2*sqrt(x2*x2+y2*y2)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
       arr[1,1]:=i+1;
       arr[2,2]:=(x1/h)*r2*sqrt(x2*x2+y2*y2)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
       arr[2,1]:=i-n+1;
       arr[3,2]:=(x2/h)*r3*sqrt(x1*x1+y1*y1)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
       arr[3,1]:=i+n;
       arr[4,2]:=((1-x2)/h)*r3*sqrt(x1*x1+y1*y1)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
       arr[4,1]:=i;
      end;
   end;

Procedure InsAkupSkrugPoint (n,i:integer;deltay,r1,r3,h:real;var F:text);
var k: real;
j,q:Longint;
enter:boolean;
arr_coeff: coeff_array;
begin
 VertSpecPoint(h,r1,r3,round(i+(0.5-sign(deltay)/2)*n),n,arr_coeff);
 for j:=1 to n*n do
  begin
   k:=0;
    if j=i then k:=k+(-1/(r3*abs(deltay))-1/(r1*h));
    if j=i-n then k:=k+(1/(r1*h));
    if j=arr_coeff[1,1] then k:=k+(1/(r3*abs(deltay)))*arr_coeff[1,2];
    if j=arr_coeff[2,1] then k:=k+(1/(r3*abs(deltay)))*arr_coeff[2,2];
    if j=arr_coeff[3,1] then k:=k+(1/(r3*abs(deltay)))*arr_coeff[3,2];
    if j=arr_coeff[4,1] then k:=k+(1/(r3*abs(deltay)))*arr_coeff[4,2];
   writeln(F,k);
  end;
end;

Procedure InsSpecMidPoint (n,i:integer;deltax,deltay,r1,r2,h:real;var F1:text);
var k: real;
j,q:Longint;
enter:boolean;
arr_coeff, arr_coeff2: coeff_array;
begin
 enter:=false;
 for j:=1 to n*n do
   begin
    k:=0;
    if (deltax<>0) and (deltay=0) then
     begin
      HorSpecPoint(h, r1, r2, round(i-0.5+0.5*sign(deltax)), n, arr_coeff);
      if j=i then k:=k+(2/(h*abs(deltax))+2/sqr(h));
      if j=i+n then k:=k+1/(h*h);
      if j=i-n then k:=k+1/(h*h);
      if j=i-sign(deltax) then k:=k+2/(h*(h+abs(deltax)));
      if j=arr_coeff[1,1] then k:=k+(2/((h+abs(deltax))*abs(deltax)))*arr_coeff[1,2];
      if j=arr_coeff[2,1] then k:=k+(2/((h+abs(deltax))*abs(deltax)))*arr_coeff[2,2];
      if j=arr_coeff[3,1] then k:=k+(2/((h+abs(deltax))*abs(deltax)))*arr_coeff[3,2];
      if j=arr_coeff[4,1] then k:=k+(2/((h+abs(deltax))*abs(deltax)))*arr_coeff[4,2];
     end;
    if (deltax=0) and (deltay<>0) then
     begin
      k:=0;
      VertSpecPoint(h,r1,r2,round(i+(0.5-sign(deltay)/2)*n),n,arr_coeff);
       if j=i then k:=k+(2/(h*abs(deltay))+2/sqr(h));
       if j=i+1 then k:=k+1/(h*h);
       if j=i-1 then k:=k+1/(h*h);
       if j=i-n*sign(deltay) then k:=k+2/(h*(h+abs(deltay)));
       if j=arr_coeff[1,1] then k:=k+(2/((h+abs(deltay))*abs(deltay)))*arr_coeff[1,2];
       if j=arr_coeff[2,1] then k:=k+(2/((h+abs(deltay))*abs(deltay)))*arr_coeff[2,2];
       if j=arr_coeff[3,1] then k:=k+(2/((h+abs(deltay))*abs(deltay)))*arr_coeff[3,2];
       if j=arr_coeff[4,1] then k:=k+(2/((h+abs(deltay))*abs(deltay)))*arr_coeff[4,2];
     end;
    if (deltax<>0) and (deltay<>0) then
     begin
        VertSpecPoint(h,r1,r2,round(i+(0.5-sign(deltay)/2)*n),n,arr_coeff);
        HorSpecPoint(h,r1,r2,round(i-(0.5-sign(deltax)/2)*n),n,arr_coeff2);
        if j=i then k:=k+(2/(h*abs(deltay))+2/(h*abs(deltax)));
        if j=i-n*sign(deltay) then k:=k+2/(h*(h+abs(deltay)));
        if j=i-sign(deltax) then k:=k+2/(h*(h+abs(deltax)));
        if j=arr_coeff[1,1] then k:=k+(2/((h+abs(deltay))*abs(deltay)))*arr_coeff[1,2];
        if j=arr_coeff[2,1] then k:=k+(2/((h+abs(deltay))*abs(deltay)))*arr_coeff[2,2];
        if j=arr_coeff[3,1] then k:=k+(2/((h+abs(deltay))*abs(deltay)))*arr_coeff[3,2];
        if j=arr_coeff[4,1] then k:=k+(2/((h+abs(deltay))*abs(deltay)))*arr_coeff[4,2];
        if j=arr_coeff2[1,1] then k:=k+(2/((h+abs(deltay))*abs(deltay)))*arr_coeff2[1,2];
        if j=arr_coeff2[2,1] then k:=k+(2/((h+abs(deltay))*abs(deltay)))*arr_coeff2[2,2];
        if j=arr_coeff2[3,1] then k:=k+(2/((h+abs(deltay))*abs(deltay)))*arr_coeff2[3,2];
        if j=arr_coeff2[4,1] then k:=k+(2/((h+abs(deltay))*abs(deltay)))*arr_coeff2[4,2];
     end;
    writeln(F1,k);
   end;
end;

Procedure HorSpecPoint2;
var j,k: integer;
    x1,x2,y1,y2,b: real;
begin
 j:=i div n;
 y1:=(dx/z)*(h-deltax);
 y2:=(dx/z)*deltax;
 x1:=deltax;
 x2:=h-deltax;
 arr[1,2]:=((1-x1)/h)*r2*sqrt(x2*x2+y2*y2)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
 arr[1,1]:=i+1;
 arr[2,2]:=(x1/h)*r2*sqrt(x2*x2+y2*y2)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
 arr[2,1]:=i-n+1;
 arr[3,2]:=(x2/h)*r3*sqrt(x1*x1+y1*y1)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
 arr[3,1]:=i+n;
 arr[4,2]:=((1-x2)/h)*r3*sqrt(x1*x1+y1*y1)/(r2*sqrt(x2*x2+y2*y2)+r3*sqrt(x1*x1+y1*y1));///q;
 arr[4,1]:=i;
end;

Procedure InsSpecMidPoint2;
var k: real;
j,q:Longint;
enter:boolean;
arr_coeff: coeff_array;
begin
 enter:=false;
 HorSpecPoint2 (h,r2,r3,deltax,i,n,arr_coeff);
 for j:=1 to n*n do
   begin
    k:=0;
      if j=i then k:=k+(2/(h*abs(deltax))+2/sqr(h));
      if j=i+n then k:=k+1/(h*h);
      if j=i-n then k:=k+1/(h*h);
      if j=i-sign(deltax) then k:=k+2/(h*(h+abs(deltax)));
      if j=arr_coeff[1,1] then k:=k+(2/((h+abs(deltax))*abs(deltax)))*arr_coeff[1,2];
      if j=arr_coeff[2,1] then k:=k+(2/((h+abs(deltax))*abs(deltax)))*arr_coeff[2,2];
      if j=arr_coeff[3,1] then k:=k+(2/((h+abs(deltax))*abs(deltax)))*arr_coeff[3,2];
      if j=arr_coeff[4,1] then k:=k+(2/((h+abs(deltax))*abs(deltax)))*arr_coeff[4,2];
      writeln(F1,k);
   end;
end;
END.