si0[ 3 ]:=par0[ 1 ]; {Alfa- ratio of approx.}
mCur:=3;
end;
setApproximationType( apDT ); {approx. type for direct problem}
setApproximationData( si0, mCur ); {approx. data for direct problem}
nApprox := ( nApprox AND $0F ); {XXXXYYYY->0000YYYY}
fHypTg := (( nApprox AND apHypTg ) = apHypTg );
fMulti := (( nApprox AND apExp ) = 0 ) AND NOT fHypTg; {It's not an EXP approx.}
if fMulti then
begin
for i:=1 to nPoints do
begin
Gr[ 1,i ]:=SiMax[ i ];
Gr[ 2,i ]:=SiMin[ i ];
Rg[ i ]:=( Gr[ 1,i ] + Gr[ 2,i ] )/2; {zero estimate of SI}
Rgs[ i ]:=1E33; {biggest integer}
end;
mLast:=nPoints; {loop for every node of approx.}
mCur :=1; {to begin from the only node of approx}
end
else
if fHypTg then
begin
Gr[ 1,1 ]:= siMax[ 1 ]; Gr[ 2,1 ]:= siMin[ 1 ]; Rgs[ 1 ]:=1E33;
Gr[ 1,2 ]:=parMax[ 2 ]; Gr[ 2,2 ]:=parMin[ 2 ]; Rgs[ 2 ]:=1E33;
Gr[ 1,3 ]:=parMax[ 3 ]; Gr[ 2,3 ]:=parMin[ 3 ]; Rgs[ 3 ]:=1E33;
Gr[ 1,4 ]:=parMax[ 4 ]; Gr[ 2,4 ]:=parMin[ 4 ]; Rgs[ 4 ]:=1E33;
for i:=1 to 4 do Rg[ i ]:=( Gr[ 1,i ] + Gr[ 2,i ] )/2;
mLast:=1;
mCur:=4;
end
else
begin
Gr[ 1,1 ]:= siMax[1]; Gr[2,1]:= siMin[1]; Rgs[ 1 ]:=1E33;
Gr[ 1,2 ]:= siMax[nPoints]; Gr[2,2]:= siMin[nPoints]; Rgs[ 2 ]:=1E33;
Gr[ 1,3 ]:= parMax[1]; Gr[2,3]:= parMin[1]; Rgs[ 3 ]:=1E33;
for i:=1 to 3 do Rg[ i ]:=( Gr[ 1,i ] + Gr[ 2,i ] )/2;
mLast:=1;
mCur :=3;
end;
initConst( nLayers, parMaxH, parMaxX , parEps, parEqlB );{set probe params}
end;
procedure directTask; {emulate voltage measurements [with error]}
begin
for i:=1 to nFreqs do
begin
getVoltage( freqs[i], Umr[ i ], Umi[ i ] ); {"measured" Uvn*}
if ( epsU > 0 ) then {add measurement error}
begin
randomize; Umr[ i ]:=Umr[ i ]*( 1 + epsU*( random-0.5 ) );
randomize; Umi[ i ]:=Umi[ i ]*( 1 + epsU*( random-0.5 ) );
end;
end;
writeln('* Voltage measurements have been emulated');
setApproximationType( nApprox ); {approx. type for inverse problem}
setApproximationData( Rg, mCur ); {approx. data for inverse problem}
end;
procedure reduceSILimits; {evaluate SI for m+1 points of approx. using aG}
var
x0, x1, xL, dx, Gr1, Gr2 : real;
j, k : byte;
begin
{----------------------------- get SI min/max for m+1 points of approximation}
dx:=1/( nPoints-1 );
for i:=1 to m+1 do
begin
k:=1;
x1:=0;
x0:=( i-1 )/m;
for j:=1 to nPoints-1 do
begin
xL:=( j-1 )/( nPoints-1 );
if( ( xL < x0 ) AND ( x0 <= xL+dx ) )then
begin
k:=j;
x1:=xL;
end;
end;
Gr[ 1,i ]:=siMax[ k ] + ( siMax[ k+1 ]-siMax[ k ] )*( x0-x1 )/dx;
Gr[ 2,i ]:=siMin[ k ] + ( siMin[ k+1 ]-siMin[ k ] )*( x0-x1 )/dx;
end;
{------------------------------------- get SI for m+1 points of approximation}
for i:=1 to m+1 do
begin
Rg[i]:=getSiFunction( (i-1)/m );
if ( Rg[i] > Gr[1,i] )then Rg[i]:=Gr[1,i];
if ( Rg[i] < Gr[2,i] )then Rg[i]:=Gr[2,i];
if m > 1 then {There're more than 1 point of approx.}
begin
Gr1:= Rg[i]+( Gr[1,i]-Rg[i] )*aG; {reduce upper bound}
Gr2:= Rg[i]-( Rg[i]-Gr[2,i] )*aG; {reduce lower bound}
if ( Gr1 < Gr[1,i] )then Gr[1,i]:=Gr1; {test overflow}
if ( Gr2 > Gr[2,i] )then Gr[2,i]:=Gr2;
end;
end;
setApproximationData( Rg , m+1 );
end;
procedure resultMessage; {to announce new results}
begin
if fMulti then
begin
writeln(' current nodal values of conductivity');
write(' si : ');for i:=1 to m do write(Rg[i] :6:3,' ');writeln;
write(' max: ');for i:=1 to m do write(Gr[1,i]:6:3,' ');writeln;
write(' min: ');for i:=1 to m do write(Gr[2,i]:6:3,' ');writeln;
end
else
begin
for i:=1 to nPoints do si[i]:=getSiFunction( ( i-1 )/( nPoints-1 ) );
if fHypTg then
saveHypTgResults
else
saveExpResults;
end;
end;
procedure clockMessage; {user-friendly message}
begin
writeln('***********************************************************');
write( '* approximation points number :',m:3,' * Time '); clock;
writeln('***********************************************************');
end;
procedure done; {final message}
begin
Sound(222); Delay(111); Sound(444); Delay(111); NoSound; {beep}
write('* Task processing time '); clock; saveTime;
writeln('* Status: Inverse problem has been successfully evaluated.');
end;
Begin
about;
loadData;
initParameters;
directTask;
for m:=1 to mLast do
begin
if fMulti then
begin
mCur:=m;
clockMessage;
end;
doMinimization; {main part of work}
setApproximationData( Rg, mCur ); {set new approx. data}
resultMessage;
if(( fMulti )AND( m < nPoints ))then reduceSILimits;
end;
done;
End.
Ï1.5 Ìîäóëü ãëîáàëüíûõ äàííûõ EData
Unit EData;
Interface
Uses DOS;
Const
maxPAR = 40; {nodes of approximation max number}
maxFUN = 20; {excitation frequencies max number}
maxSPC = 4; {support approximation values number}
iterImax = 50; {max number of internal iterations}
Const
apSpline = 1; {approximation type identifiers}
apHypTg = 3;
apExp = 2;
apPWCon = 4;
apPWLin = 8;
Type
Parameters = array[ 1..maxPAR ] of real; {Si,Mu data}
Functionals = array[ 1..maxFUN ] of real; {Voltage data}
SpecialPar = array[ 1..maxSPC ] of real; {Special data}
Var
hThick : real; {thickness of slab}
nPoints : integer; {nodes of approximation number within [ 0,H ]}
nLayers : integer; {number of piecewise constant SI layers within[ 0,H ]}
nFreqs : integer; {number of excitation frequencies}
nStab : integer; {required number of true digits in SI estimation}
epsU : real; {relative error of measurement Uvn*}
nApprox : byte; {approximation type identifier}
incVal : real; {step for numerical differ.}
parMaxH : integer; {max number of integration steps}
parMaxX : real; {upper bound of integration}
parEps : real; {error of integration}
derivType: byte; {1 for right; 2 for central}
Var
freqs : Functionals; {frequencies of excitment Uvn*}
Umr, Umi : Functionals; { Re(Uvn*),Im(Uvn*) for "measured" data}
Uer, Uei : Functionals; { Re(Uvn*),Im(Uvn*) for estimated data}
mu : Parameters; {relative permeability nodal values}
si, si0 : Parameters; {conductivity approximation nodal values}
siMin, siMax : Parameters; {conductivity nodal values restrictions}
par0 : SpecialPar; {alfa,si2,beta,gamma - for exp&HypTg}
parMin,parMax: SpecialPar; - min/max
zLayer : Parameters; {relative borders of slab layers [0,1]}
Var
aG : real; {scale factor for SImin/max}
Ft : real; {current discrepancy functional value}
fMulti : boolean; {TRUE if it isn't an EXP-approximation}
fHypTg : boolean; {TRUE for Hyperbolic tg approximation}
parEqlB : boolean; {TRUE if b[i]=const}
mCur : integer; {current number of approximation nodes}
inFileName : string; {data file name}
outFileName : string; {results file name}
Var
Rg : Parameters; {current SI estimation}
RgS : Parameters; {previous SI estimation}
Gr : array [ 1..2 ,1..maxPAR ] of real; {SI max/min}
Fh : array [ 1..maxPAR , 1..maxFUN ] of real; {current discrepancies}
Type
TTime=record
H, M, S, S100 : word; {hour,min,sec,sec/100}
end;
Var
clk1, clk2 : TTime; {start&finish time}
procedure loadData; {load all user-defined data from file}
Implementation
procedure loadData;
var
i,eqlB : integer;
FF : text;
begin
assign( FF, outFileName ); {clear output file}
rewrite( FF );
close( FF );
assign( FF, inFileName ); {read input file}
reset( FF );
readln( FF );
Ñòðàíèöû: 1, 2, 3, 4, 5, 6, 7, 8
