Решение обратной задачи вихретокового контроля

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

     readln( FF );

     readln( FF, hThick, nPoints, nLayers, nFreqs, nStab, epsU, aG, nApprox );

     readln( FF );

     for i:=1 to nFreqs do read( FF, freqs[i] );

     readln( FF );

     readln( FF );

     readln( FF );

     for i:=1 to nPoints do readln( FF, si[i], siMin[i], siMax[i] );

     readln( FF );

     readln( FF );

     readln( FF , incVal, parMaxH, parMaxX, parEps, derivType, eqlB );

     readln( FF );

     readln( FF );

     for i:=1 to maxSPC do readln( FF, par0[i] , parMin[i] , parMax[i] );

     readln( FF );

     if ( eqlB=0 )then

     begin

          for i:=1 to nLayers+1 do read( FF, zLayer[i] );

          parEqlB:=false;

     end

     else parEqlB:=true;

     close( FF );

     for i:=1 to maxPAR do mu[i]:=1;

end;

Var

     str : string;

Begin

     if( ParamCount = 1 )then str:=ParamStr(1)

     else

     begin

          write('Enter I/O file name, please: ');

          readln( str );

     end;

     inFileName :=str+'.txt';

     outFileName:=str+'.lst';

End.

П1.6 Модуль работы с файлами EFile

Unit EFile;

Interface

Uses

    DOS, EData;

function  isStable( ns : integer; var RG1,RG2 ) : boolean;

function  saveResults( ns,iter : integer ) : boolean;

procedure saveExpResults;

procedure saveHypTgResults;

procedure clock;

procedure saveTime;

Implementation

Var

   FF : text;

   i  : byte;

function decimalDegree( n:integer ) : real;{10^n}

var

     s:real; i:byte;

begin

     s:=1;

     for i:=1 to n do s:=s*10;

     decimalDegree:=s;

end;

function isStable( ns:integer ; var RG1,RG2 ) : boolean;

var

     m  : real;

     R1 : Parameters absolute RG1;

     R2 : Parameters absolute RG2;

begin

     isStable:=TRUE;

     m:=decimalDegree( ns-1 );

     for i:=1 to mCur do

     begin

          if NOT(( ABS( R2[i]-R1[i] )*m )<=ABS( R2[i]) ) then isStable:=FALSE;

          RgS[i]:=Rg[i];

     end;

end;

function saveResults( ns , iter : integer ) : boolean;

var

     sum : real;

begin

     sum:=0;

     for i:=1 to nFreqs do sum:=sum + Fh[1,i];

     sum:=SQRT( sum/nFreqs );

     assign( FF , outFileName );

     append( FF );

     write(      iter:2, ' <”>', sum:10:7, ' Rg=' );

     write( FF , iter:2, ' <”>', sum:10:7, ' Rg=');

     for i:=1 to mCur do

     begin

          write(      Rg[i]:6:3, ' ');

          write( FF , Rg[i]:6:3, ' ');

     end;

     writeln;

     writeln( FF );

     close( FF );

     saveResults:=isStable( ns , Rgs , Rg );

end;

procedure saveExpResults;

begin

     assign( FF , outFileName );

     append( FF );

     writeln(      ' siE=',Rg[2]:6:3,' siI=',Rg[1]:6:3,' alfa=',Rg[3]:6:3);

     writeln( FF , ' siE=',Rg[2]:6:3,' siI=',Rg[1]:6:3,' alfa=',Rg[3]:6:3);

     write(      ' SI: ');

     write( FF , ' SI: ');

     for i:=1 to nPoints do

     begin

          write(      si[i]:6:3,' ');

          write( FF , si[i]:6:3,' ');

     end;

     writeln;

     writeln( FF );

     close( FF );

end;

procedure saveHypTgResults;

begin

     assign( FF , outFileName );

     append( FF );

     writeln(      ' si1=',Rg[2]:6:3,' si2=',Rg[1]:6:3,' beta=',Rg[3]:6:3,' gamma=',Rg[4]:6:3);

     writeln( FF , ' si1=',Rg[2]:6:3,' si2=',Rg[1]:6:3,' beta=',Rg[3]:6:3,' gamma=',Rg[4]:6:3);

     write(      ' SI: ');

     write( FF , ' SI: ');

     for i:=1 to nPoints do

     begin

          write(      si[i]:6:3,' ');

          write( FF , si[i]:6:3,' ');

     end;

     writeln;

     writeln( FF );

     close( FF );

end;

procedure clock;                                                  {t2 = t2-t1}

var

     H1,M1,S1,H2,M2,S2,sec1,sec2 : longint;

begin

     GetTime( clk2.H, clk2.M, clk2.S, clk2.S100 );              {current time}

     H2:=clk2.H; M2:=clk2.M; S2:=clk2.S; H1:=clk1.H; M1:=clk1.M; S1:=clk1.S;

     sec2:= ( H2*60 + M2 )*60 + S2;

     sec1:= ( H1*60 + M1 )*60 + S1;

     if( sec2 < sec1 )then sec2:=sec2 + 85020;                     {+23.59.59}

     sec2:=sec2 - sec1;

     clk2.H := sec2 div 3600; sec2:=sec2 - clk2.H*3600;

     clk2.M := sec2 div 60;   sec2:=sec2 - clk2.M*60;

     clk2.S := sec2;

     writeln( clk2.H:2, ':', clk2.M:2, ':', clk2.S:2 );

end;

procedure saveTime;

begin

     assign( FF , outFileName );

     append( FF );

     write( FF ,'* Processing time ',clk2.H:2, ':', clk2.M:2, ':', clk2.S:2 );

     close( FF );

end;

End.

П1.7 Модуль решения прямой задачи ВТК для НВТП EDirect

{****************************************************************************}

{ ERIN submodule : EDirect , 15.02.99,    (C) 1999 by Nikita U.Dolgov        }

{****************************************************************************}

{  Estimates Uvn* for Eddy current testing of inhomogeneous multilayer slab  }

{ with  surface( flat ) probe.                                               }

{  It can do it using one of five types of conductivity approximation :      }

{Spline, Exponential, Piecewise constant, Piecewise linear,Hyperbolic tangent}

{****************************************************************************}

{$F+}

Unit EDirect;

Interface

Uses EData, EMath;

Type

     siFunc   = function( x:real ) : real;

Var

     getSiFunction : siFunc; {for external getting SI estimate}

procedure initConst( par1,par2:integer; par3,par4:real; par5:boolean );

procedure getVoltage( freq : real ; var ur,ui : real );   { Uvn* = ur + j*ui }

procedure setApproximationType( approx : byte );          { type of approx.  }

procedure setApproximationItem( SIG:real ; N : byte );    { set SIGMA[ N ]}

procedure setApproximationData( var SIG; nVal : byte );   { SIGMA[1..nVal]   }

procedure getApproximationData( var SIG ; var N : byte ); { get SIGMA[ N ]}

Implementation

Const

     PI23 = 2000*pi;          {2*pi*KHz}

     mu0  = 4*pi*1E-7;        {magnetic const}

Var

     appSigma : Parameters;   {conductivity approximation data buffer}

     appCount : byte;         {size of conductivity approximation data buffer}

     appType  : byte;         {conductivity approximation type identifier}

Type

     commonInfo=record

                      w       : real;    {cyclical excitation frecuency}

                      R       : real;    {equivalent radius of probe}

                      H       : real;    {generalized lift-off of probe}

                      Kr      : real;    {parameter of probe}

                      eps     : real;    {error of integration}

                      xMax    : real;    {upper bound of integration}

                      steps   : integer; {current number of integration steps}

                      maxsteps: integer; {max number of integration steps}

                      Nlay    : integer;    {number of layers in slab}

                      sigma   : Parameters; {conductivity of layers}

                      m       : Parameters; {relative permeability of layers}

                      b       : Parameters; {thickness of layers}

                      zCentre : Parameters; {centre of layer}

                   end;

     procFunc = procedure( x:real; var result:complex);

Var

     siB, siC, siD : Parameters;          {support for Spline approx.}

     cInfo         : commonInfo;          {one-way access low level info}

function siSpline( x:real ) : real;{Spline approximation}

begin

     if( appCount = 1 )then

         siSpline := appSigma[ 1 ]

     else

         siSpline:=Seval( appCount, x, appSigma, siB, siC, siD);

end;

function siExp( x:real ) : real;{Exponential approximation}

begin

     siExp:=(appSigma[2]-appSigma[1])*EXP( -appSigma[3]*(1-x) ) + appSigma[1];

end;

function siPWConst( x:real ) : real;{Piecewise constant approximation}

var

     dx, dh : real; i : byte;

begin

     if( appCount = 1 )then siPWConst := appSigma[ 1 ]

     else

     begin

          dh:=1/( appCount-1 );

          dx:=dh/2;

          i:=1;

          while( x > dx ) do

          begin

               i:=i + 1;

               dx:=dx + dh;

          end;

          siPWConst:=appSigma[ i ];

     end;

end;

function siPWLinear( x:real ) : real;{Piecewise linear approximation}

var

     dx, dh : real;

     i      : byte;

begin

     if( appCount = 1 )then siPWLinear := appSigma[ 1 ]

     else

     begin

          dh:=1/( appCount-1 );

          dx:=0;

          i:=1;

          repeat

                i:=i + 1;

                dx:=dx + dh;

          until( x <= dx );

          siPWLinear:=appSigma[i-1]+( appSigma[i]-appSigma[i-1] )*( x/dh+2-i);

     end;

end;

function siHyperTg( x:real ) : real;{Hyperbolic tangent approximation}

begin

     siHyperTg:=appSigma[2]+(appSigma[1]-appSigma[2])*(1+th((appSigma[3]-x)/appSigma[4]))/2;

end;

procedure setApproximationType( approx : byte );

begin

     appType := approx;

     write('*  conductivity approximation type : ');

     case approx of

                   apSpline : begin

                                   writeln('SPLINE');

                                   getSiFunction := siSpline;

                              end;

                   apExp    : begin

                                   writeln('EXP');

                                   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}

              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 getVoltage( freqs[j],ur2[j],ui2[j] );

                   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:=( ur2[j]-ur )/( 2*Rg[i]*incVal );

                        dui:=( ui2[j]-ui )/( 2*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;

     end;

     setApproximationData( si , k );

end;

procedure doMinimization;

const

     mp1Max = maxPAR + 1;

     mp2Max = maxPAR + 2;

     m2Max  = 2*( maxPAR + maxFUN );

     m21Max = m2Max + 1;

     n2Max  = 2*maxFUN;

     m1Max  = maxPAR + n2Max;

     n1Max  = n2Max + 1;

     mm1Max = maxPAR + n1Max;

     minDh  : real = 0.001;  {criterion of an exit from golden section method}

var

     A : array [ 1 .. m1Max , 1 .. m21Max ] of real;

     B : array [ 1 .. m1Max]  of real;             

     Sx: array [ 1 .. m21Max] of real;             

     Zt: array [ 1 .. maxPAR] of real;

     Nb: array [ 1 .. m1Max]  of integer;

     N0: array [ 1 .. m21Max] of integer;          

     a1, a2, dh, r, tt, tp, tl, cv, cv1, cl, cp : real;

     n2, n1, mp1, mp2, mm1, m1, m2, m21 : integer;    

     ain   : real;   

     apn   : real;   

     iq    : integer;

     k1    : integer;

     n11   : integer;

     ip    : integer;

     iterI : integer;

     i,j,k : integer;

label

     102 ,103 ,104 ,105 ,106 ,107 ,108;

begin

     n2:=2*nFreqs; n1:=n2+1; m1:=mCur+n2;

     mp1:=mCur+1; mp2:=mCur+2; mm1:=mCur+n1;

     m2:=2*( mCur+nFreqs ); m21:=m2+1;

     for k:=1 to m1Max do

         for i:=1 to m21Max do

             A[k,i]:=0;                   

     iterI:=0;

102:

     iterI:=iterI+1;

     getFunctional( 0 );                  

     for i:=1 to nFreqs do b[i]:= -Fh[1,i];

     getFunctional( derivType );          

     for k:=1 to mCur do

     begin

          Zt[k]:=Rg[k];

          for i:=1 to nFreqs do

          begin

               A[i,k+1]:=Fh[k,i];         

               A[i+nFreqs,k+1]:=-A[i,k+1];

          end;

          for i:=1 to nFreqs do B[i]:=B[i]+Rg[k]*A[i,k+1];

     end;

     for i:=1 to nFreqs do B[i+nFreqs]:=-B[i];

     for i:=n1 to m1 do B[i]:=Gr[1,i-n2]-Gr[2,i-n2];   

     for i:=1 to m1 do

     begin

          if i<=n2 then

             for k:=2 to mp1 do B[i]:=B[i]-A[i,k]*Gr[2,k-1];

          A[i,1]:=-1;

          if i>n2 then

          begin

               A[i,1]:=0;

               for k:=2 to mp1 do           

                   if i-n2=k-1 then A[i,k]:=1

                   else A[i,k]:=0;

          end;

          for k:=mp2 to m21 do              

              if k-mp1=i then A[i,k]:=1

              else A[i,k]:=0;

     end;

     k1:=1;

     for k:=1 to n2 do

         if B[k1]>B[k] then k1:=k;        

     for k:=1 to mp1 do A[k1,k]:=-A[k1,k];

     A[k1,mCur+1+k1]:=0;

     B[k1]:=-B[k1];

     for i:=1 to n2 do

         if i<>k1 then

         begin

              B[i]:=B[i]+B[k1];

              for k:=1 to mm1 do A[i,k]:=A[i,k]+A[k1,k];

         end;

     for i:=mp2 to m21 do

     begin

          Sx[i]:=B[i-mp1];

          Nb[i-mp1]:=i;

     end;

     for i:=1 to mp1 do Sx[i]:=0;

     Sx[1]:=B[k1];

     Sx[mp1+k1]:=0;

     Nb[k1]:=1;

103:

     for i:=2 to m21 do N0[i]:=0;       

104:

     for i:=m21 downto 2 do

         if N0[i]=0 then n11:=i;

     for k:=2 to m21 do

         if ((A[k1,n11]<A[k1,k]) AND (k<>N0[k])) then n11:=k;

     if A[k1,n11]<=0 then goto 105;       

     iq:=0;

     for i:=1 to m1 do

         if i<>k1 then

         begin

              if A[i,n11]>0 then

              begin

                   iq:=iq+1;

                   if iq=1 then

                   begin

                        Sx[n11]:=B[i]/A[i,n11]; ip:=i;

                   end

                   else

                   begin

                        if Sx[n11]>B[i]/A[i,n11] then

                        begin

                             Sx[n11]:=B[i]/A[i,n11]; ip:=i;

                        end;

                   end;

              end

              else

                  if iq=0 then

                  begin

                       N0[n11]:=n11;

                       goto 104;

                  end;

         end;

     Sx[Nb[ip]]:=0;

     Nb[ip]:=n11;

     B[ip]:=B[ip]/A[ip,n11];

     apn:=A[ip,n11];

     for k:=2 to m21 do A[ip,k]:=A[ip,k]/apn;

     for i:=1 to m1 do

         if i<>ip then

         begin

              ain:=A[i,n11];

              B[i]:=-B[ip]*ain+B[i];

              for j:=1 to m21 do A[i,j]:=-ain*A[ip,j]+A[i,j];

         end;

     for i:=1 to m1 do Sx[Nb[i]]:=B[i];

     goto 103;

105:

     for k:=1 to mCur do Sx[k+1]:=Sx[k+1]+Gr[2,k];

     a1:=0;

     a2:=1.;

     dh:=a2-a1;

     r:=0.618033;

     tl:=a1+r*r*dh;

     tp:=a1+r*dh;

     j:=1;

108:

     if j=1 then tt:=tl else tt:=tp;

106:

     for i:=1 to mCur do Rg[i]:=Zt[i]+tt*(Sx[i+1]-Zt[i]);

     getFunctional( 0 );                                

     cv:=abs(Fh[1,1]);

     if nFreqs>1 then

        for k:=2 to nFreqs do

        begin

             cv1:=abs(Fh[1,k]);

             if cv<cv1 then cv:=cv1;

        end;

     if (j=1) or (j=3) then cl:=cv

     else cp:=cv;

     if j=1 then

     begin

          j:=2;

          goto 108;

     end;

     if dh<MinDh then goto 107;

     if cl>cp then

     begin

          a1:=tl; dh:=a2-a1; tl:=tp; tp:=a1+r*dh  ; tt:=tp; cl:=cp; j:=4;

     end

     else

     begin

          a2:=tp; dh:=tp-a1; tp:=tl; tl:=a1+r*r*dh; tt:=tl; cp:=cl; j:=3;

     end;

     goto 106;

107:

     if (iterI < iterImax)AND(NOT saveResults( nStab,iterI )) then goto 102;

end;

End.

Приложение 2 - Удельная электрическая проводимость материалов

Приведем сводку справочных данных согласно[7-9].

Материал

smin ,[МСм/м]

smax ,[МСм/м]

Немагнитные стали

0.4

1.8

Бронзы (БрБ, Бр2, Бр9)

6.8

17

Латуни (ЛС59, ЛС62)

13.5

17.8

Магниевые сплавы (МЛ5-МЛ15)

5.8

18.5

Титановые сплавы (ОТ4, ВТ3-ВТ16)

0.48

2.15

Алюминиевые сплавы (В95, Д16, Д19)

15.1

26.9