getSiFunction := siExp;
end;
apPWCon : begin
writeln('PIECEWISE CONST');
getSiFunction := siPWConst;
end;
apPWLin : begin
writeln('PIECEWISE LINEAR');
getSiFunction := siPWLinear;
end;
apHypTg : begin
writeln('HYPERBOLIC TANGENT');
getSiFunction := siHyperTg;
end;
end;
end;
procedure setApproximationData( var SIG ; nVal : byte );
var
Sigma : Parameters absolute SIG; i:byte;
begin
appCount := nVal;
for i:=1 to nVal do appSigma[ i ]:=Sigma[ i ];
if( appType = apSpline )then Spline( appCount, appSigma, siB, siC, siD);
end;
procedure getApproximationData( var SIG ; var N : byte );
var
Sigma : Parameters absolute SIG; i:byte;
begin
N := appCount;
for i:=1 to appCount do Sigma[ i ]:=appSigma[ i ];
end;
procedure setApproximationItem( SIG:real ; N : byte );
begin
appSigma[ N ] := SIG;
if( appType = apSpline )then Spline( appCount, appSigma, siB, siC, siD);
end;
procedure functionFi( x:real ; var result:complex );{get boundary conditions function value}
var
beta : array[ 1..maxPAR ]of real;
q : array[ 1..maxPAR ]of complex;
fi : array[ 0..maxPAR ]of complex;
th , z1 , z2 , z3 , z4 , z5 , z6 , z7 : complex;
i : byte;
begin
mkComp( 0, 0, fi[0] );
with cInfo do
for i:=1 to Nlay do
begin
beta[i]:=R*sqrt( w*mu0*sigma[i] ); {calculation of beta}
mkComp( sqr(x), sqr(beta[i])*m[i], z7 ); {calculation of q, z7=q^2}
SqrtC( z7, q[i] );
mulComp( q[i], b[i], z6 ); {calculation of th,z6=q*b}
tanH( z6, th ); {th=tanH(q*b)}
mkComp( sqr(m[i])*sqr(x), 0, z6 ); {z6=m2x2}
SubC( z6, z7, z5); {z5=m2x2-q2}
AddC( z6, z7, z4); {z4=m2x2+q2}
MulC( z5, th, z1); {z1=z5*th}
MulC( z4, th, z2); {z2=z4*th}
mulComp( q[i], 2*x*m[i], z3 ); {z3=2xmq}
SubC( z2, z3, z4 );
MulC( z4, fi[i-1], z5 );
SubC( z1, z5, z6 ); {z6=high}
AddC( z2, z3, z4 );
MulC( z1, fi[i-1], z5 );
SubC( z4, z5, th ); {th=low}
DivC( z6, th, fi[i] );
end;
eqlComp( result, fi[ cInfo.Nlay ] );
end;
procedure funcSimple( x:real; var result:complex );{intergrand function value}
var
z : complex;
begin
with cInfo do
begin
functionFi( x, result );
mulComp( result, exp( -x*H ), z );
mulComp( z, J1( x*Kr ), result );
mulComp( result, J1( x/Kr ), z );
eqlComp( result, z );
end;
end;
procedure funcMax( x:real; var result:complex );{max value; When Fi[Nlay]=1}
var
z1, z2 : complex;
begin
with cInfo do
begin
mkComp( 1,0,z1 );
mulComp( z1, exp(-x*H), z2 );
mulComp( z2, J1( x*Kr ), z1 );
mulComp( z1, J1( x/Kr ), result );
end;
end;
procedure integralBS( func:procFunc ; var result:complex );{integral by Simpson}
var
z , y , tmp : complex;
hh : real;
i, iLast : word;
begin
with cInfo do
begin
hh:=xMax/steps;
iLast:=steps div 2;
nullComp(tmp);
func( 0, z );
eqlComp( y, z );
for i:=1 to iLast do
begin
func( ( 2*i-1 )*hh , z );
deltaComp( tmp, z );
end;
mulComp( tmp, 4, z );
deltaComp( y, z );
nullComp( tmp );
iLast:=iLast-1;
for i:=1 to iLast do
begin
func( 2*i*hh ,z );
deltaComp( tmp, z );
end;
mulComp( tmp, 2, z );
deltaComp( y, z );
func( xmax, z );
deltaComp(y,z);
mulComp( y, hh/3, z );
eqlComp( result, z );
end;
end;{I = h/3 * [ F0 + 4*sum(F2k-1) + 2*sum(F2k) + Fn ]}
procedure integral( F:procFunc; var result:complex );{integral with given error}
var
e , e15 : real;
flag : boolean;
delta , integ1H , integ2H : complex;
begin
with cInfo do
begin
e15 :=eps*15; I2h-I1h
steps:=20;
flag :=false;
integralBS( F, integ2H );
repeat
eqlComp( integ1H, integ2H );
steps:=steps*2;
integralBS( F, integ2H );
SubC( integ2H, integ1H, delta );
e:=Leng( delta );
if( e<e15 )then flag:=true;
until( ( flag ) OR (steps>maxsteps) );
if( flag )then
begin
eqlComp( result, integ2H );
end
else
begin
writeln('Error: Too big number of integration steps.');
halt(1);
end;
end;
end;
procedure initConst( par1, par2 : integer; par3, par4 : real; par5:boolean );
var
i : byte;
bThick, dl, x : real;
const
Ri=0.02; hi=0.005; { radius and lift-off of excitation coil}
Rm=0.02; hm=0.005; { radius and lift-off of measuring coil}
begin
with cInfo do
begin
Nlay :=par1;
xMax :=par3;
maxsteps:=par2;
R :=sqrt( Ri*Rm );
H :=( hi+hm )/R;
Kr :=sqrt( Ri/Rm );
eps :=par4;
bThick :=hThick*0.002/R; {2*b/R [m]}
for i:=1 to Nlay do m[i]:= mu[i];
if par5 then
begin
bThick:=bThick/NLay;
for i:=1 to Nlay do b[i]:=bThick;
dl:=1/NLay;
x:=dl/2; {x grows up from bottom of slab to the top}
for i:=1 to Nlay do
begin
zCentre[i]:=x;
x:=x + dl;
end;
end
else
for i:=1 to Nlay do
begin
b[i]:=( zLayer[i+1]-zLayer[i] )*bThick;
zCentre[i]:=( zLayer[i+1]+zLayer[i] )/2;
end;
end;
end;
procedure init( f:real );{get current approach of conductivity values}
var
i : byte;
begin
with cInfo do
begin
w:=PI23*f;
for i:=1 to Nlay do sigma[i]:=getSiFunction( zCentre[i] )*1E6;
end;
end;
procedure getVoltage( freq : real ; var ur,ui : real );
var
U, U0, Uvn, tmp : complex;
begin
init( freq );
integral( funcSimple, U ); { U =Uvn }
integral( funcMax , U0 ); { U0=Uvn max }
divComp( U, Leng(U0), Uvn ); U0
mkComp( 0, 1, tmp ); { tmp=( 0+j1 ) }
MulC( tmp, Uvn, U ); { U= j*Uvn = Uvn* }
ur := U.re;
ui := U.im;
end;
END.
Ï1.8 Ìîäóëü ðåøåíèÿ îáðàòíîé çàäà÷è ÂÒÊ äëÿ ÍÂÒÏ EMinimum
Unit EMinimum;
INTERFACE
Uses EData, Crt, EFile, EDirect;
procedure doMinimization;
IMPLEMENTATION
procedure getFunctional( Reg : byte );
var
ur, ui, dur, dui, Rgt : real;
ur2, ui2: Functionals;
i, j, k : byte;
begin
getApproximationData( si , k );
setApproximationData( Rg, mCur );
case Reg of
0 : for i:=1 to nFreqs do {get functional F value}
begin
getVoltage( freqs[i], ur, ui );
Uer[ i ]:=ur; {we need it for dU}
Uei[ i ]:=ui;
Fh[1,i] := SQR( ur-Umr[i] ) + SQR( ui-Umi[i] );
end;
{Right:U'(i)= (U(i+1)-U(i))/h}
1 : for i:=1 to mCur do {get dF/dSI[i] value}
begin
Rgt:=Rg[i]*( 1+incVal ); {si[i]=si[i]+dsi[i]}
setApproximationItem( Rgt, i ); {set new si[i] value}
for j:=1 to nFreqs do
begin {get dUr/dSI,dUi/dSI}
getVoltage( freqs[ j ], ur, ui );
dur:=( ur-Uer[j] )/( Rg[i]*incVal );
dui:=( ui-Uei[j] )/( Rg[i]*incVal );
Fh[i,j]:=2*(dur*(Uer[j]-Umr[j])+dui*(Uei[j]-Umi[j]));
end;
setApproximationItem( Rg[i], i ); {restore si[i] value}
end;
{Central:U'(i)= (U(i+1)-U(i-1))/2h}
2 : for i:=1 to mCur do {get dF/dSI[i] value}
Ñòðàíèöû: 1, 2, 3, 4, 5, 6, 7, 8
