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].