Шпунт В.Х. Динамические электрические свойства кожи человека// Медицинская техника.- 1997.- № 4 -С.38.

Dumitrescu J., Constantin D. Acupuncture stuntifica moderna / Iasi, Editura junimea, 1977. - 383 с.

Fere C. Note sur des modifications de la resistance electrique sous e`influence des excitations sensorielles et des emotions, C.R. Hebd. Seances e.Memories Soc. Biol. - 1888, V.40. - p. 217 - 219.

Приложение 1

Текст основной программы.

program PR_Otlad_MY;

{$APPTYPE CONSOLE}

uses

Fil_XYZ_otlad_MY;

const

SaveFileName1='u_z1_ot_MY.prn';

SaveFileName2='u_y1_ot_MY.prn';

SaveFileName3='u_x1_ot_MY.prn';

SaveFileName4='gamma_ot_MY.prn';

SaveFileName5='plotnost_ot_MY.prn';

SaveFileName6=' plotnost2_ot_MY.prn';

procedure SetInitial;

const

Tab=chr(9);

var

i, j, k : word;

uu : real;

f:text;

begin

for k:=1 to Mk do

begin

uu:=(k-1)*(abs(-0.8+0.7)/50)-0.8;

for i:=1 to Mx do

for j:=1 to Mk do

u[i,j,k]:=uu;

end;

{ assign(f,SaveFileName5);

rewrite(f);

for i:=1 to Mx do

begin

for j:=1 to My do

begin

write(f,u[i,j,20]:16);

if i<Mx+1 then write(f,Tab);

end;

writeln(f,'');

end; }

end;{--------------------SetInitial--------------------}

procedure SaveMatrixToFile;

const

Tab=chr(9);

var

i, j, k : word;

f : text;

begin

assign(f,SaveFileName1);

rewrite(f);

for i:=1 to Mx do

begin

for j:=1 to My do

begin

write(f,u[i,j,20]:16);

if i<Mx+1 then write(f,Tab);

end;

writeln(f,'');

end;

close(f);

assign(f,SaveFileName2);

rewrite(f);

for i:=1 to Mx do

begin

for k:=Mk downto 1 do

begin

write(f,u[i,25,k]:16);

if i<Mx+1 then write(f,Tab);

end;

writeln(f,'');

end;

close(f);

assign(f,SaveFileName3);

rewrite(f);

for j:=1 to My do

begin

for k:=Mk downto 1 do

begin

write(f,u[35,j,k]:16);

if j<My+1 then write(f,Tab);

end;

writeln(f,'');

end;

close(f);

assign(f,SaveFileName4);

rewrite(f);

for j:=1 to My do

begin

for k:=1 to Mk do

begin

write(f,gamma[35,j,k]:16);

if j<Mx then write(f,Tab);

end;

writeln(f,'');

end;

close(f);

assign(f,SaveFileName5);

rewrite(f);

for j:=1 to My do

begin

for k:=Mk downto 1 do

begin

write(f,((U[35,j,k]-U[35,j,k+1])*gamma[35,j,k]/10):16);

if j<Mx then write(f,Tab);

end;

writeln(f,'');

end;

close(f);

assign(f,SaveFileName6);

rewrite(f);

for i:=1 to Mx do

begin

for k:=Mk downto 1 do

begin

write(f,((U[i,25,k]-U[i,25,k+1])*gamma[i,25,k]/10):16);

if i<Mx then write(f,Tab);

end;

writeln(f,'');

end;

close(f);

end;{--------------------SaveMatrixToFile--------------------}

Begin

SetParameters;

SetInitial;

Compute;

SaveMatrixToFile;

End.

Текст модуля.

unit Unit_otlad_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=434e-6; {stop value for relative deviation}

rlc=12; {real output digits amount}

STau=0.9; {value for relaxation parameter}

gamma1=1e-5; {--epidermis--}

gamma2=0.1; {--Sosochk.sloi--}

gamma3=5; {--vascular--}

gamma4=0.01; {--setchat.sloi--}

gamma5=5; {--nervous cells--}

gamma6=0.1; {--subcutaneus fat--}

gamma7=1; {--muscle--}

gamma8=1; {--pora--}

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;

implementation

procedure SetParameters;

var

i, j, k : word;

begin

tau:=STau;

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[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>=25)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>=10)and(k<20)))or

((((i-44)*(i-44)+(j-19)*(j-19))<=9)and((k>=10)and(k<20)))

then gamma[i,j,k]:=gamma5;

end;

end;

procedure FreeBoundary;

var

i, j, k : integer;

delta:single;

begin

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;

procedure Compute;

var

IterNum, i, j, k : integer;

begin

IterNum:=0;

while {IterNum<800} sd>sdStop do

begin {-----main loop-----}

sad:=0;

for i:=2 to Mx{-1} do

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+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;

readln;

end;

Begin

End.

Текст модуля для расчета макромодели (модели пальца).

unit Fil_XYZ_Palec;

interface{-------------------------------------------------------------}

const

Nx=45; {net nodes amount on X axe}

Ny=45; {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=1e-6; {stop value for relative deviation}

rlc=13; {real output digits amount}

STau=0.9; {value for relaxation parameter}

El_potential=0.1; {Volts}

El_current1=5e-2; {Current density,A/m^2}

El_current2=1.6e-3; {Current density,A/m^2}

S=2e-5; {m*m}

l=5e-3; {m}

dz=0.5e-3; {step,m}

gamma0=0;

gamma1=1e-7; {--Skin--}

gamma2=0.1; {--Subcutaneous fat--}

gamma3=1; {--m--}

gamma4=1e-7; {--k--}

gamma5=1e-6; {--n--}

gamma10=6e7; {electrode}

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

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;

{deviation sum control parameters}

sd:=1e+9;

sadMax:=Mr*Mz*Mk; sabMax:=2*(Mz+Mr+Mk);

for i:=1 to Mr do

for j:=1 to Mz do

for k:=1 to Mk do

begin {----------------------gamma----------------------}

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))

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)and(k<=133)or

(((i-21)*(i-21)+(j-21)*(j-21)+(k-133)*(k-133))<144))

then gamma[i,j,k]:=gamma3;

if ((((i-21)*(i-21)+(j-21)*(j-21))<=64)and(k<=133)or

(((i-21)*(i-21)+(j-21)*(j-21)+(k-133)*(k-133))<64))

then gamma[i,j,k]:=gamma4;

if ((((i-21)*(i-21)+(j-21)*(j-21))<=225)and(i<20)and(((j-21)*(j-21)+(k-137)*(k-137))<=36)

and(((i-21)*(i-21)+(j-21)*(j-21))>169))

then gamma[i,j,k]:=gamma5;

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))<=400})) 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))<=400))

then gamma[i,j,k]:=gamma10;

end;

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

procedure FreeBoundary;

var

UG, G, Ug2, G2, Ji1, ji2 : single;

i, j, k, p, q : 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,Mk]:=0;

end;

{-----Electrode Surface-----}

G:=0;

G2:=0;

Ug:=0;

Ug2:=0;

p:=1;

q:=1;

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)and(k<60)and(((i-21)*(i-21)+(j-21)*(j-21))<=256))

then

G2:=G2+1;

end;

Ji2:=El_current2/G2;

for i:=1 to Mr do

for j:=1 to Mz do

for k:=1 to Mk do

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 begin

Upredict:=(gamma1 * U[i+p,j+p,k] - Ji2*p*dz) / gamma1;

u1[i,j,k]:=u[i,j,k]*(1-tau)+tau*upredict;

sab:=sab+abs(upredict-u[i,j,k]);

end

else if (i>21)and(j<=21) then begin

Upredict:=(gamma1 * U[i-p,j+p,k] - Ji2*p*dz) / gamma1;

u1[i,j,k]:=u[i,j,k]*(1-tau)+tau*upredict;

sab:=sab+abs(upredict-u[i,j,k]);

end

else if (i<=21)and(j>21) then begin

Upredict:=(gamma1 * U[i+p,j-p,k] - Ji2*p*dz) / gamma1;

u1[i,j,k]:=u[i,j,k]*(1-tau)+tau*upredict;

sab:=sab+abs(upredict-u[i,j,k]);

end

else if (i>21)and(j>21) then begin

Upredict:=(gamma1 * U[i-p,j-p,k] - Ji2*p*dz) / gamma1;

u1[i,j,k]:=u[i,j,k]*(1-tau)+tau*upredict;

sab:=sab+abs(upredict-u[i,j,k]);

end;

Ug2:=Ug2+U1[i,j,k];

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>115)and(k<=125)and(i>=7)and(i<17)and(j<20)and(((i-21)*(i-21)+(j-21)*(j-21))<=256))

then

G:=G+1;

end;

Ji1:=El_current1/G;

for i:=1 to Mr do

for j:=1 to Mz do

for k:=1 to Mk do

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))<=256))

then begin

Upredict:=(gamma1* U[i+q,j+q,k] - Ji1*q*dz) / gamma1;

u1[i,j,k]:=u[i,j,k]*(1-tau)+tau*upredict;

sab:=sab+abs(upredict-u[i,j,k]);

Ug:=Ug+U1[i,j,k];

end;

writeln (Ug/G,' ',Ug2/G2,' ',Ug2/G2-Ug/G);

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

Procedure Current;

var

Ji, Ji2, 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

if ((((i-21)*(i-21)+(j-21)*(j-21))>=225)and(k>115)and(k<=125)and(i>=7)and(i<17)and(j<21)and(((i-21)*(i-21)+(j-21)*(j-21))<=256))

then begin

Ji:=gamma1*(u[i+1,j+1,k]-u[i,j,k])/(1*dz);

Cd:=Cd+Ji;

end;

Cur:=Cd*pi*sqr(5e-3)/4;

write(Cd,Cur);

readln;

Cd:=0;

for i:=1 to Mr do

for j:=1 to Mz do

for k:=1 to Mk do

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

Ji2:=gamma1*(u[i+1,j+1,k]-u[i,j,k])/(1*dz)

else if (i>21)and(j<=21) then

Ji2:=gamma1*(u[i-1,j+1,k]-u[i,j,k])/(1*dz)

else if (i<=21)and(j>21) then

Ji2:=gamma1*(u[i+1,j-1,k]-u[i,j,k])/(1*dz)

else if (i>21)and(j>21) then

Ji2:=gamma1*(u[i-1,j-1,k]-u[i,j,k])/(1*dz);

Cd:=Cd+Ji2;

end;

Cur:=Cd*pi*2e-2*1e-2;

write(Cd,Cur);

readln;

end;{---------------------Current Computing--------------------}

procedure Compute;

var

IterNum, i, j, k : integer;

begin

IterNum:=0;

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

while 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+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.

11.1. Текст программы для расчета распределения электрического поля в декартовой системе координат.

program V_pryam;

{$APPTYPE CONSOLE}

uses

pryam_r10 in 'pryam_r10.pas';

const

SaveFileName1='ux_pryam.prn';

SaveFileName2='uy_pryam.prn';

SaveFileName3='uz_pryam.prn';

SaveFileName4='gammay.prn';

SaveFileName5='gammax.prn';

SaveFileName6='gammaz.prn';

procedure SetInitial;

const

Tab=chr(9);

var

i, j, k : word;

uu : real; f:text;

begin

for i:=1 to Mx Do

for j:=1 to My Do

for k:=1 to Mz Do

Begin

u[i,j,k] := 1;

End;

end;{--------------------SetInitial--------------------}

procedure SaveMatrixToFile;

const

Tab=chr(9);

var

i, j, k : word;

f : text;

begin

assign(f,SaveFileName3);

rewrite(f);

for i:=2 to Mx-1 do

begin

for j:=2 to My-1 do

begin

write(f,u[i,j,50]:16);

if i<Mx-1 then write(f,Tab);

end; writeln(f,'');

end; close(f);

assign(f,SaveFileName2);

rewrite(f);

for i:=2 to Mx-2 do

begin

for k:=2 to Mz-2 do

begin

write(f,u[i,25,k]:16);

if i<Mx-1 then write(f,Tab);

end;

writeln(f,''); end; close(f);

assign(f,SaveFileName1);

rewrite(f);

for j:=2 to My-2 do

begin

for k:=2 to Mz-2 do

begin

write(f,u[25,j,k]:16);

if j<My-1 then write(f,Tab);

end; writeln(f,'');

end; close(f);

assign(f,SaveFileName4);

rewrite(f);

for k:=1 to Mz do

begin

for i:=1 to Mx do

begin

write(f,gamma[i,26,k]:16);

if k<Mz+1 then write(f,Tab);

end;

writeln(f,''); end; close(f);

assign(f,SaveFileName5);

rewrite(f);

for k:=1 to Mz do

begin

for j:=1 to Mx do

begin

write(f,gamma[26,j,k]:16);

if k<Mz+1 then write(f,Tab);

end;

writeln(f,''); end; close(f);

assign(f,SaveFileName6);

rewrite(f);

for j:=1 to My do

begin

for i:=1 to Mx do

begin

write(f,gamma[i,j,40]:16);

if k<Mz then write(f,Tab);

end;

writeln(f,''); end; close(f);

end;{--------------------SaveMatrixToFile--------------------}

Begin

SetParameters;

SetInitial;

Compute;

SaveMatrixToFile;

End.

unit pryam_r10;

interface{------------------------------------------------------}

const

Nx=60;

Ny=60;

Nz=100;

Mx=Nx;

My=Ny;

Mz=Nz;

rlc=12;

STau=0.9;

gamma1=2e-12;

gamma2=10;

type

U_Matrix = array[1..Mx,1..My,1..Mz] of single;

GammaArray = array[1..Mx,1..My,1..Mz] of single;

var

tau : double;

num : integer;

u, u1 : U_Matrix;

upredict: double;

gamma : GammaArray;

sd, sad, sadMax: double;

procedure SetParameters;

procedure Compute;

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

procedure SetParameters;

var

i, j, k : word;

begin {----------------------Set parameters--------------------}

tau:=STau;

sd:=0;

sadMax:=Mx*My*Mz;

for i:=1 to Mx do

for j:=1 to My do

for k:=1 to Mz do

gamma[i,j,k]:=gamma1;

for k:=1 to Mz do

for j:=1 to My do

for i:=1 to Mx do

begin

if ((i-31)*(i-31)+(j-31)*(j-31)<=25*25)and(k>=7)and(k<=94) then gamma[i,j,k]:=gamma2

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

for k:=1 to Mz do

begin

u1[1,j,k]:=0;

u1[Mx,j,k]:=0;

u1[i,1,k]:=0;

u1[i,My,k]:=0;

u1[i,j,Mz]:=0;

u1[i,j,1]:=0;

end;

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

procedure Compute;

var

IterNum, i, j, k : integer;

begin

IterNum:=0;

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

while IterNum<1000 do

begin {-----main loop-----}

sad:=0;

for i:=2 to Mx do {-----inner points scanning-----}

for j:=2 to My do

for k:=2 to Mz 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+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.

11.2. Текст программы для расчета распределения электрического поля в цилиндрической системе координат.

program v_cylindrical;

{$APPTYPE CONSOLE}

uses

cylindrical in 'cylindrical_proga.pas';

const

SaveFileName1='ur_cylindr.prn';

SaveFileName2='ua_cylindr.prn';

SaveFileName3='uz_cylindr.prn';

SaveFileName4='gammaa_cylindr.prn';

SaveFileName5='gammar_cylindr.prn';

SaveFileName6='gammaz_cylindr.prn';

var

Upr:array[1..100,1..100,1..100] of real;

procedure SetInitial;

const

Tab=chr(9);

var

r,a,z : integer{word};

begin

for r:=1 to Mr Do

for a:=1 to Ma Do

for z:=1 to Mz Do

Begin

u[r,a,z]:= 1;

End;

end;{--------------------SetInitial--------------------}

procedure Perevod;

var a,r,z :word;

begin {----------------------Set parameters----------------------}

sd:=0;

for r:=1 to Mr do

for a:=0 to Ma-1 do

for z:=1 to Mz do

Upr[round(r*cos(a)),round(r*sin(a)),z]:=U2[r,a,z];

for r:=1 to Mr do

for a:=0 to Ma-1 do

for z:=1 to Mz do

if (u2[r,a,z]-u1[r,a,z]<0.001) or (u2[r,a,z]-u1[r,a,z]>-0.001) then

u2[r,a,z]:=u[r,a,z];

end;

procedure SaveMatrixToFile;

const

Tab=chr(9);

var

r,a,z{,i, j, k} : word;

f : text;

begin

assign(f,SaveFileName3);

rewrite(f);

for r:=1 to Mr do

begin

for a:=1 to Ma do

begin

write(f,u[r,a,50]:16);

if r<Mr+1 then write(f,Tab);

end; writeln(f,'');

end; close(f);

assign(f,SaveFileName2);

rewrite(f);

for r:=1 to Mr do

begin

for z:=1 to Mz do

begin

write(f,u[r,15,z]:16);

if r<Mr+1 then write(f,Tab);

end;

writeln(f,''); end; close(f);

assign(f,SaveFileName1);

rewrite(f);

for a:=1 to Ma do

begin

for z:=1 to Mz do

begin

write(f,u[15,a,z]:16);

if z<Mz+1 then write(f,Tab);

end; writeln(f,'');

end; close(f);

assign(f,SaveFileName4);

rewrite(f);

for z:=1 to Mz do

begin

for r:=1 to Mr do

begin

write(f,gamma[r,15,z]:16);

if z<Mz+1 then write(f,Tab);

end;

writeln(f,''); end; close(f);

assign(f,SaveFileName5);

rewrite(f);

for z:=1 to Mz do

begin

for a:=1 to Ma do

begin

write(f,gamma[15,a,z]:16);

if z<Mz+1 then write(f,Tab);

end;

writeln(f,''); end; close(f);

assign(f,SaveFileName6);

rewrite(f);

for a:=1 to Ma do

begin

for r:=1 to Mr do

begin

write(f,gamma[r,a,15]:16);

if z<Mz+1 then write(f,Tab);

end;

writeln(f,''); end; close(f);

end;{--------------------SaveMatrixToFile--------------------}

Begin

SetParameters;

SetInitial;

Compute;

perevod;

SaveMatrixToFile;

End.

unit cylindrical_proga;

interface{------------------------------------------------------}

const

Nr=25;

Na=360;

Nz=100;

Mr=Nr;

Ma=Na;

Mz=Nz;

rlc=12;

STau=0.9;

gamma1=2e-10;

gamma2=1;

type

U_Matrix = array[1..Mr,1..Ma,1..Mz] of single;

GammaArray = array[1..Mr,1..Ma,1..Mz] of single;

var

tau : double;

num : integer;

u, u1,u2 : U_Matrix;

upredict: double;

gamma : GammaArray;

sd, sad, sadMax: double;

procedure SetParameters;

procedure Compute;

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

procedure SetParameters;

var

r,z,a : word;

begin {--------------------Set parameters----------------------}

tau:=STau;

sd:=0;

sadMax:=Mr*Ma*Mz;

for a:=1 to Ma do

for r:=0 to Mr do

for z:=1 to Mz do

gamma[r,a,z]:=0;

for z:=1 to Mz do

for r:=1 to Mr do

for a:=1 to Ma do

begin

if ((r<15)and(z>=7)and(z<=94)) then gamma[r,a,z]:=gamma2

else gamma[r,a,z]:=gamma1;

end;

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

procedure FreeBoundary;

var

r, a, z :integer;

begin {-----------------------FreeBoundary----------------------}

for r:=1 to Mr do

for a:=0 to Ma do

for z:=1 to Mz do

begin

u1[r,a,1]:=0;

u1[r,a,Mz]:=0;

u1[Mr,a,z]:=0;

U1[r,a,7]:=10;

U1[r,a,94]:=-10;

end;

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

procedure Compute;

var

IterNum,r, a, z : integer;

h,h1:real;

begin

IterNum:=0;

h:=1;

h1:=1;{--------------------computation itself------------------}

while IterNum<70 do

begin {-----main loop-----}

sad:=0;

for r:=2 to Mr do {-----inner points scanning-----}

for a:=2 to Ma do

for z:=2 to Mz do

begin

if gamma[r+1,a,z]<>gamma[r,a,z] then

upredict:=( gamma[r,a,z] * U[r-1,a,z] + gamma[r+1,a,z] * U[r+1,a,z] ) /

( gamma[r+1,a,z] + gamma[r,a,z] ) else

if gamma[r,a,z]<>gamma[r-1,a,z] then

upredict:=( gamma[r-1,a,z] * U[r-1,a,z] + gamma[r,a,z] * U[r+1,a,z] ) /

( gamma[r-1,a,z] + gamma[r,a,z] ) else

if gamma[r,a,z]<>gamma[r,a,z-1] then

upredict:=( gamma[r,a,z-1] * U[r,a,z-1] + gamma[r,a,z] * U[r,a,z+1] ) /

( gamma[r,a,z-1] + gamma[r,a,z] ) else

if gamma[r,a,z+1]<>gamma[r,a,z] then

upredict:=( gamma[r,a,z] * U[r,a,z-1] + gamma[r,a,z+1] * U[r,a,z+1] ) /

( gamma[r,a,z+1] + gamma[r,a,z] ) else

upredict:= ( r*r*h1*h1*U[r+1,a,z] +r*r*h1*h1* U[r-1,a,z] - r*h*h1* U[r-1,a,z] +h*h* U[r,a+1,z] +

h*h*U[r,a-1,z] + r*r*h1*h1*U[r,a,z+1]+ r*r*h1*h1*U[r,a,z-1] ) /(4*r*r*h1*h1-r*h*h1*h1);

// upredict:= ( sqr(r)*sqr(h1)*U[r+1,a,z] +sqr(r)*sqr(h1)* U[r-1,a,z] - r*h*h1* U[r-1,a,z] +sqr(h)* U[r,a+1,z] +

// sqr(h)*U[r,a-1,z] + sqr(r)*sqr(h1)*U[r,a,z+1]+ sqr(r)*sqr(h1)*U[r,a,z-1] ) / 4*sqr(r)*sqr(h1)-r*h*sqr(h1); //

u1[r,a,z]:=tau*upredict+(1-tau)*u[r,a,z];

sad:=sad+abs(upredict-u[r,a,z]);

end;

FreeBoundary;

sd:=sad/sadMax;

IterNum:=IterNum+1;

u:=u1;

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

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

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

Begin

End.

11.3. Текст модуля для расчета емкостной составляющей ЭИ в декартовой системе координат (для тест-объекта: цилиндр).

unit MainForm;

interface

uses

Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,

Dialogs, StdCtrls;

type

TfmMainForm = class(TForm)

ButtonSet: TButton;

ButtonCount: TButton;

ButtonDraw: TButton;

ButtonClose: TButton;

Function Setka(Num : Real) : Integer;

procedure ButtonSetClick(Sender: TObject);

procedure ButtonCountClick(Sender: TObject);

procedure ButtonDrawClick(Sender: TObject);

procedure ButtonCloseClick(Sender: TObject);

private

{ Private declarations }

public

{ Public declarations }

end;

Point3D = Record

X : Integer;

Y : Integer;

Z : Integer;

End;

Const

DimX = 0.03; //Размеры пространства по оси Х

DimY = 0.03;

DimZ = 0.075;

Step = 0.0005; //Шаг сетки

nX = Trunc(DimX/Step); //Количество шагов по оси X

nY = Trunc(DimY/Step); //Trunc - чтобы целое число было, иначе будет

nZ = Trunc(DimZ/Step); // считаться дробным, даже если поделится нацело.

// Округляется в сторону нуля (в меньшую сторону)

R = 0.01; //Радиус цилиндра

R2 = R + Step; //Радиус окружностей вокруг цилиндра (следующий шаг)

CenterX = 0.015; //Центр цилиндра по оси X

CenterY = 0.015; //Центр цилиндра по оси Y

Z1 = 0.002; //Нижняя граница цилиндра

Z2 = 0.072; //Верхняя граница цилиндра

UPoverh = 2; //Потенциал точек на поверхности цилиндра

UNad = 0.000001; //Потенциал точек непосредственно над поверхностью цилиндра

var

fmMainForm: TfmMainForm;

//Массив векторов напряженностей:

// E : Array[0..nX, 0..nY, 0..nZ] Of Point3D;

//Массив значений потенциалов в каждой точке:

U : Array[0..nX, 0..nY, 0..nZ] Of Double;

k : integer;

implementation

{$R *.dfm}

Function TfmMainForm.Setka(Num : Real) : Integer;

Begin

//Для преобразования обычных координат (например, в миллиметрах) в дискретные

Result := Trunc(Num/Step);

End;

procedure TfmMainForm.ButtonSetClick(Sender: TObject);

Var

CycleX, CycleY, CycleZ : Integer;

begin

If Setka(Z1) <= Setka(0) Then //Если Z1 - не в самом низу нашего пространства

Begin

ShowMessage('Не могу произвести установку значений - нижняя граница цилиндра совпадает с нижней границей пространства!');

Exit;

End;

If Setka(Z2) >= Setka(DimZ) Then //Если Z2 - не в самом верху нашего пространства

Begin

ShowMessage('Не могу произвести установку значений - верхняя граница цилиндра совпадает с верхней границей пространства!');

Exit;

End;

If (CenterX + R2 > DimX) Or (CenterX - R2 < 0) Or

(CenterY + R2 > DimY) Or (CenterY - R2 < 0) Then

Begin

ShowMessage('Не могу произвести установку значений - проверьте размеры пространства и размеры цилиндра!');

Exit;

End;

For CycleX := 0 To nX Do

For CycleY := 0 To nY Do

For CycleZ := 0 To nZ Do

U[CycleX, CycleY, CycleZ] := 0; //Обнуление потенциалов всех точек

k:=0;

//Боковые стенки цилиндра и пространство непосредственно над цилиндром

For CycleX := 0 To nX Do

For CycleY := 0 To nY Do

For CycleZ := Setka(Z1 + Step) To Setka(Z2 - Step) Do

Begin

//Проверка на принадлежность окружности с радиусом R - принадлежность

// интервалу [R - Step/2, R + Step/2)

If (Sqr(CycleX - Setka(CenterX)) + Sqr(CycleY - Setka(CenterY)) >= Sqr(Setka(R - Step/2))) And

(Sqr(CycleX - Setka(CenterX)) + Sqr(CycleY - Setka(CenterY)) < Sqr(Setka(R + Step/2))) Then

Begin

U[CycleX, CycleY, CycleZ] := UPoverh; //Потенциал точек поверхности цилиндра

k:=k+1;

End;

//Проверка на принадлежность окружности с радиусом R2 = R + Step -

// принадлежность интервалу [R2 - Step/2, R2 + Step/2)

If (Sqr(CycleX - Setka(CenterX)) + Sqr(CycleY - Setka(CenterY)) >= Sqr(Setka(R2 - Step/2))) And

(Sqr(CycleX - Setka(CenterX)) + Sqr(CycleY - Setka(CenterY)) < Sqr(Setka(R2 + Step/2)))

Then

U[CycleX, CycleY, CycleZ] := UNad; //Потенциал точек непросредственно над цилиндром

End;

//Основания цилиндра

For CycleX := 0 To nX Do

For CycleY := 0 To nY Do

If (Sqr(CycleX - Setka(CenterX)) + Sqr(CycleY - Setka(CenterY)) < Sqr(Setka(R + Step/2))) Then

Begin

U[CycleX, CycleY, Setka(Z1)] := UPoverh; //Потенциал точек нижнего основания цилиндра

U[CycleX, CycleY, Setka(Z2)] := UPoverh; //Потенциал точек верхнего основания цилиндра

U[CycleX, CycleY, Setka(Z1 - Step)] := UNad; //Потенциал точек под нижним основанием цилиндра

U[CycleX, CycleY, Setka(Z2 + Step)] := UNad; //Потенциал точек над верхним основанием цилиндра

End;

end;

procedure TfmMainForm.ButtonCountClick(Sender: TObject); //расчет заряда

Var

CycleX, CycleY, CycleZ: Integer;

// CycleX2, CycleY2, CycleZ2 : Integer;

h : Double; //Растояние между поверхностью и точками непосредственно над поверхностью

E : Double; //Напряженность

dS : Double; //h*h

QBok : Double; //Суммарный заряд боковой поверхности цилиндра

QOsnov : Double; //Суммарный заряд оснований цилиндра

Q,Q1,C: Double; //Суммарный заряд

begin

{ Алгоритм: бежим по каждой точке, для которой потенциал равен UPoverh и ищем

ближайшую к ней (но не дальше, чем на один шаг) точку с потенциалом UNad.

Высчитываем расстояние h между точками и напряженность E = (UPoverh-UNad) / h}

If Setka(Z1) <= Setka(0) Then //Если Z1 - не в самом низу нашего пространства

Begin

ShowMessage('Не могу произвести расчет - нижняя граница цилиндра совпадает с нижней границей пространства!');

Exit;

End;

If Setka(Z2) >= Setka(DimZ) Then //Если Z2 - не в самом верху нашего пространства

Begin

ShowMessage('Не могу произвести расчет - верхняя граница цилиндра совпадает с верхней границей пространства!');

Exit;

End;

If (CenterX + R2 > DimX) Or (CenterX - R2 < 0) Or

(CenterY + R2 > DimY) Or (CenterY - R2 < 0) Then

Begin

ShowMessage('Не могу произвести расчет - проверьте размеры пространства и размеры цилиндра!');

Exit;

End;

QBok := 0;

For CycleX := 0 To nX Do

For CycleY := 0 To nY Do

For CycleZ := Setka(Z1 + Step) To Setka(Z2 + Step) Do //Бежим по высотам, в пределах которых лежит цилиндр

Begin

If U[CycleX, CycleY, CycleZ] = UPoverh Then

Begin

//Смотрим, есть ли точка на расстоянии Step от текущей со значением потенциала = UNad

If (U[CycleX + Setka(Step), CycleY, CycleZ] = UNad) Or

(U[CycleX - Setka(Step), CycleY, CycleZ] = UNad) Or

(U[CycleX, CycleY + Setka(Step), CycleZ] = UNad) Or

(U[CycleX, CycleY - Setka(Step), CycleZ] = UNad) Then //Если есть

Begin

h := Step;

End

Else //Если нет, то ближайшая точка находится на расстоянии Step * корень из 2

Begin

h := Step * Sqrt(2);

End;

dS := 2*pi*R*(Z2-Z1)/k {h*h};

E := (UPoverh - UNad) / h;

QBok := QBok + E * dS;

End;

End;

h := Step;

dS := 2*pi*R*(Z2-Z1)/k {h*h};

E := (UPoverh - UNad) / h;

{ QOsnov := 0;

For CycleX := 0 To nX Do

For CycleY := 0 To nY Do

If U[CycleX, CycleY, Setka(Z1)] = UPoverh Then

QOsnov := QOsnov + E * dS;

QOsnov := QOsnov * 2;//Два основания с одинаковым зарядом }

Q := {QOsnov + }QBok;

C:=0;

Q1:=Q*8.8542e-12;

C:=Q1/(UPoverh - UNad);

C:=C/50;

ShowMessage('Суммарный заряд: ' + FormatFloat('0.000E+00', Q1) + #13#10 +

'Емкость по теореме Гаусса: ' + FormatFloat('0.000E+00', C){ + #13#10 +

': ' + FormatFloat('0.000E+00', Q)});

end;

11.4. Текст модуля для расчета емкостной составляющей ЭИ в декартовой системе координат (для тест-объекта: шар).

unit Unit_shar;

interface

uses

Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,

Dialogs, StdCtrls;

type

TPoint3D = Record //Точка в трехмерном пространстве

X, Y, Z : Integer;

End;

TfmForm1 = class(TForm)

Button1: TButton;

procedure Button1Click(Sender: TObject);

private

{ Private declarations }

public

{ Public declarations }

end;

var

fmForm1: TfmForm1;

Arr : Array[1..50, 1..50, 1..50] Of Integer; //Массив из 50х50х50 элементов

Napr1 : Array[1..50, 1..50, 1..50] Of TPoint3D;

//50х50х50 (для каждой точки) направляющих векторов длиной = 1 - для

// поиска координат следующей точки

E : Array[1..50, 1..50, 1..50] Of Double; //Напряженность для каждой точки

implementation

{$R *.dfm}

procedure TfmForm1.Button1Click(Sender: TObject);

Var

Cycle, Cycle2, Cycle3,k : Integer; //Для организации циклов

LengthOfNapr : Double; //Длина направляющего вектора Napr

H,dh,r,S,ds : Double; //Длина вектора Napr1, так как она не всегда равна 1

Q1,Q,C : Double; //Общий заряд при dS = h*h = 1*1

Q2,C1,C2 : Double; //Общий заряд при dS = H*H

begin

For Cycle := 1 To 50 Do

For Cycle2 := 1 To 50 Do

For Cycle3 := 1 To 50 Do

Arr[Cycle, Cycle2, Cycle3] := 0;

k:=0;

S:=0;

ds:=0;

r:=10e-3;

For Cycle := 1 To 50 Do

For Cycle2 := 1 To 50 Do

For Cycle3 := 1 To 50 Do

Begin

If ((Cycle-25) * (Cycle-25) + (Cycle2-25) * (Cycle2-25) + (Cycle3-25) * (Cycle3-25) <= 121) And

((Cycle-25) * (Cycle-25) + (Cycle2-25) * (Cycle2-25) + (Cycle3-25) * (Cycle3-25) > 100)Then

//Проверка на принадлежность сфере - все точки в пределах (10, 11]

begin

k:=k+1;

Arr[Cycle, Cycle2, Cycle3] := 1;

end; //Устанавливаем в 1 - как признак того, что точка принадлежит сфере

End;

S:=4*pi*r*r;

ds:=S/k;

dh:=0.0001;

Q1 := 0;

Q2 := 0;

C1 := 0;

C2 := 0;

For Cycle := 1 To 50 Do

For Cycle2 := 1 To 50 Do

For Cycle3 := 1 To 50 Do

If Arr[Cycle, Cycle2, Cycle3] = 1 Then

Begin

LengthOfNapr := Sqrt(Sqr(Cycle - 25) + Sqr(Cycle2 - 25) + Sqr(Cycle3 - 25));

//Корень из суммы квадратов: Sqr - квадрат, Sqrt - корень

//Координаты точек направляющего вектора:

// начало - центр сферы (25, 25, 25)

// конец - точка на сфере (Cycle, Cycle2, Cycle3)

Napr1[Cycle, Cycle2, Cycle3].X := Round((Cycle - 25) / LengthOfNapr);

Napr1[Cycle, Cycle2, Cycle3].Y := Round((Cycle2 - 25) / LengthOfNapr);

Napr1[Cycle, Cycle2, Cycle3].Z := Round((Cycle3 - 25) / LengthOfNapr);

H := Sqrt(Sqr(Napr1[Cycle, Cycle2, Cycle3].X) +

Sqr(Napr1[Cycle, Cycle2, Cycle3].Y) +

Sqr(Napr1[Cycle, Cycle2, Cycle3].Z))*dh;

If H <> 0 Then

Begin

E[Cycle, Cycle2, Cycle3] := 2 / H;

Q1 := Q1 + E[Cycle, Cycle2, Cycle3] * ds;

Q2 := Q2 + E[Cycle, Cycle2, Cycle3]* dh * dh;

End;

End;

{ Q2:=(Q2*24.6176e-6)/dh*dh;}

Q1:=Q1*8.854e-12;

Q2:=Q2*8.854e-12;

C1:=Q1/2;

C2:=Q2/2;

C:=4*pi*r*8.854e-12*1.00058;

ShowMessage('Интегральный заряд сферы: Q1 = ' + FormatFloat('0.000E+00 ', Q1) + #10#13 +

{' Q2 = ' + FormatFloat('0.000E+00 ', Q2)+#10#13 +'По теореме Гаусса (ds=h*h): C1 = ' +FormatFloat('0.000E+00', C1) + #10#13 +}

'По теореме Гаусса (ds=S/n): C = ' + FormatFloat('0.000E+00 ', C2)+'(Ф)'+ #10#13 +

'По формуле (С=4*pi*eps*R): C = ' + FormatFloat('0.000E+00 ', C)+'(Ф)');

end;

end.

Размещено на Allbest.ru

...