Logo Search packages:      
Sourcecode: necpp version File versions  Download package

ground.c

/******* Translated to the C language by N. Kyriazis  20 Aug 2003 *******/
/*                                                    */
/* Program NEC(input,tape5=input,output,tape11,tape12,tape13,tape14,    */
/* tape15,tape16,tape20,tape21)                                   */
/*                                                    */
/* Numerical Electromagnetics Code (NEC2)  developed at Lawrence  */
/* Livermore lab., Livermore, CA.  (contact G. Burke at 415-422-8414    */
/* for problems with the NEC code. For problems with the vax implem-    */
/* entation, contact J. Breakall at 415-422-8196 or E. Domning at 415   */
/* 422-5936)                                                */
/* file created 4/11/80.                                    */
/*                                                    */
/*                ***********Notice**********                     */
/* This computer code material was prepared as an account of work       */
/* sponsored by the United States government.  Neither the United       */
/* States nor the United States Department Of Energy, nor any of  */
/* their employees, nor any of their contractors, subcontractors,       */
/* or their employees, makes any warranty, express or implied, or */
/* assumes any legal liability or responsibility for the accuracy,      */
/* completeness or usefulness of any information, apparatus, product    */
/* or process disclosed, or represents that its use would not infringe  */
/* privately-owned rights.                                  */
/*                                                    */
/************************************************************************/

#include "nec2c.h"

/* pointers to input/output files */
extern FILE *input_fp, *output_fp, *plot_fp;

/* common  /data/ */
extern int n, np, m, mp, ipsym, npm, np2m, np3m; /* n+m,n+2m,n+3m */
extern int *icon1, *icon2, *itag;
extern long double *x, *y, *z, *si, *bi;
extern long double *x2, *y2, *z2, *cab, *sab, *salp;
extern long double *t1x, *t1y, *t1z, *t2x, *t2y, *t2z;
extern long double *px, *py, *pz, *pbi, *psalp;
extern long double wlam;

/* common  /dataj/ */
extern int iexk, ind1, indd1, ind2, indd2, ipgnd;
extern long double s, b, xj, yj, zj, cabj, sabj, salpj, rkh;
extern long double t1xj, t1yj, t1zj, t2xj, t2yj, t2zj;
extern complex long double  exk, eyk, ezk, exs, eys, ezs, exc, eyc, ezc;

/* for the sommerfeld ground option, rom2 integrates over the source */
/* common  /incom/ */
extern int isnor;
extern long double xo, yo, zo, sn, xsn, ysn;

/* common  /gwav/ */
extern long double r1, r2, zmh, zph;
extern complex long double u, u2, xx1, xx2;

/* common  /gnd/ */
extern int ksymp, ifar, iperf, nradl;
extern long double t2, cl, ch, scrwl, scrwr;
extern complex long double zrati, zrati2, t1, frati;

/*-------------------------------------------------------------------*/

/* segment to obtain the total field due to ground.  the method of */
/* variable interval width romberg integration is used.  there are 9 */
/* field components - the x, y, and z components due to constant, */
/* sine, and cosine current distributions. */
void rom2( long double a, long double b, complex long double *sum, long double dmin )

{
      int i, ns, nt, flag=TRUE;
      int nts = 4, nx = 1, n = 9;
      long double ze, ep, zend, dz=0., dzot=0., tmag1, tmag2, tr, ti;
      long double z, s; /***also global***/
      long double rx = 1.0e-4;
      complex long double g1[9], g2[9], g3[9], g4[9], g5[9];
      complex long double t00, t01[9], t10[9], t02, t11, t20[9];

      z= a;
      ze= b;
      s= b- a;

      if( s < 0.)
      {
            fprintf( output_fp, "\n  ERROR - B LESS THAN A IN ROM2" );
            stop(-1);
      }

      ep= s/(1.e4* npm);
      zend= ze- ep;

      for( i = 0; i < n; i++ )
            sum[i]=CPLX_00;

      ns= nx;
      nt=0;
      sflds( z, g1);

      while( TRUE )
      {
            if( flag == TRUE)
            {
                  dz= s/ ns;
                  if( z+ dz > ze)
                  {
                        dz= ze- z;
                        if( dz <= ep)
                              return;
                  }

                  dzot= dz*.5;
                  sflds( z+ dzot, g3);
                  sflds( z+ dz, g5);

            } /* if( flag ) */

            tmag1=0.;
            tmag2=0.;

            /* evaluate 3 point romberg result and test convergence. */
            for( i = 0; i < n; i++ )
            {
                  t00=( g1[i]+ g5[i])* dzot;
                  t01[i]=( t00+ dz* g3[i])*.5;
                  t10[i]=(4.* t01[i]- t00)/3.;
                  if( i > 2)
                        continue;

                  tr= creal( t01[i]);
                  ti= cimag( t01[i]);
                  tmag1= tmag1+ tr* tr+ ti* ti;
                  tr= creal( t10[i]);
                  ti= cimag( t10[i]);
                  tmag2= tmag2+ tr* tr+ ti* ti;

            } /* for( i = 0; i < n; i++ ) */

            tmag1= sqrtl( tmag1);
            tmag2= sqrtl( tmag2);
            test( tmag1, tmag2, &tr, 0., 0., &ti, dmin);

            if( tr <= rx)
            {
                  for( i = 0; i < n; i++ )
                        sum[i] += t10[i];
                  nt += 2;

                  z += dz;
                  if( z > zend)
                        return;

                  for( i = 0; i < n; i++ )
                        g1[i]= g5[i];

                  if( (nt >= nts) && (ns > nx) )
                  {
                        ns= ns/2;
                        nt=1;
                  }
                  flag = TRUE;
                  continue;

            } /* if( tr <= rx) */

            sflds( z+ dz*.25, g2);
            sflds( z+ dz*.75, g4);
            tmag1=0.;
            tmag2=0.;

            /* evaluate 5 point romberg result and test convergence. */
            for( i = 0; i < n; i++ )
            {
                  t02=( t01[i]+ dzot*( g2[i]+ g4[i]))*.5;
                  t11=( 4.0 * t02- t01[i] )/3.;
                  t20[i]=(16.* t11- t10[i])/15.;
                  if( i > 2)
                        continue;

                  tr= creal( t11);
                  ti= cimag( t11);
                  tmag1= tmag1+ tr* tr+ ti* ti;
                  tr= creal( t20[i]);
                  ti= cimag( t20[i]);
                  tmag2= tmag2+ tr* tr+ ti* ti;

            } /* for( i = 0; i < n; i++ ) */

            tmag1= sqrtl( tmag1);
            tmag2= sqrtl( tmag2);
            test( tmag1, tmag2, &tr, 0.,0., &ti, dmin);

            if( tr > rx)
            {
                  nt=0;
                  if( ns < npm )
                  {
                        ns= ns*2;
                        dz= s/ ns;
                        dzot= dz*.5;

                        for( i = 0; i < n; i++ )
                        {
                              g5[i]= g3[i];
                              g3[i]= g2[i];
                        }

                        flag=FALSE;
                        continue;

                  } /* if( ns < npm) */

                  fprintf( output_fp,
                           "\n  ROM2 -- STEP SIZE LIMITED AT Z = %12.5LE", z );

            } /* if( tr > rx) */

            for( i = 0; i < n; i++ )
                  sum[i]= sum[i]+ t20[i];
            nt= nt+1;

            z= z+ dz;
            if( z > zend)
                  return;

            for( i = 0; i < n; i++ )
                  g1[i]= g5[i];

            flag = TRUE;
            if( (nt < nts) || (ns <= nx) )
                  continue;

            ns= ns/2;
            nt=1;

      } /* while( TRUE ) */

}
/*-----------------------------------------------------------------------*/

/* sfldx returns the field due to ground for a current element on */
/* the source segment at t relative to the segment center. */
void sflds( long double t, complex long double *e )

{
      long double xt, yt, zt, rhx, rhy, rhs, rho, phx, phy;
      long double cph, sph, zphs, r2s, rk, sfac, thet;
      complex long double  erv, ezv, erh, ezh, eph, er, et, hrv, hzv, hrh;

      xt= xj+ t* cabj;
      yt= yj+ t* sabj;
      zt= zj+ t* salpj;
      rhx= xo- xt;
      rhy= yo- yt;
      rhs= rhx* rhx+ rhy* rhy;
      rho= sqrtl( rhs);

      if( rho <= 0.)
      {
            rhx=1.;
            rhy=0.;
            phx=0.;
            phy=1.;
      }
      else
      {
            rhx= rhx/ rho;
            rhy= rhy/ rho;
            phx=- rhy;
            phy= rhx;
      }

      cph= rhx* xsn+ rhy* ysn;
      sph= rhy* xsn- rhx* ysn;

      if( fabsl( cph) < 1.0e-10)
            cph=0.;
      if( fabsl( sph) < 1.0e-10)
            sph=0.;

      zph= zo+ zt;
      zphs= zph* zph;
      r2s= rhs+ zphs;
      r2= sqrtl( r2s);
      rk= r2* TP;
      xx2= cmplx( cosl( rk),- sinl( rk));

      /* use norton approximation for field due to ground.  current is */
      /* lumped at segment center with current moment for constant, sine, */
      /* or cosine distribution. */
      if( isnor != 1)
      {
            zmh=1.;
            r1=1.;
            xx1=0.;
            gwave( &erv, &ezv, &erh, &ezh, &eph);

            et=-CONST1* frati* xx2/( r2s* r2);
            er=2.* et* cmplx(1.0, rk);
            et= et* cmplx(1.0 - rk* rk, rk);
            hrv=( er+ et)* rho* zph/ r2s;
            hzv=( zphs* er- rhs* et)/ r2s;
            hrh=( rhs* er- zphs* et)/ r2s;
            erv= erv- hrv;
            ezv= ezv- hzv;
            erh= erh+ hrh;
            ezh= ezh+ hrv;
            eph= eph+ et;
            erv= erv* salpj;
            ezv= ezv* salpj;
            erh= erh* sn* cph;
            ezh= ezh* sn* cph;
            eph= eph* sn* sph;
            erh= erv+ erh;
            e[0]=( erh* rhx+ eph* phx)* s;
            e[1]=( erh* rhy+ eph* phy)* s;
            e[2]=( ezv+ ezh)* s;
            e[3]=0.;
            e[4]=0.;
            e[5]=0.;
            sfac= PI* s;
            sfac= sinl( sfac)/ sfac;
            e[6]= e[0]* sfac;
            e[7]= e[1]* sfac;
            e[8]= e[2]* sfac;

            return;
      } /* if( isnor != 1) */

      /* interpolate in sommerfeld field tables */
      if( rho >= 1.0e-12)
            thet= atanl( zph/ rho);
      else
            thet= POT;

      /* combine vertical and horizontal components and convert */
      /* to x,y,z components. multiply by exp(-jkr)/r. */
      intrp( r2, thet, &erv, &ezv, &erh, &eph );
      xx2= xx2/ r2;
      sfac= sn* cph;
      erh= xx2*( salpj* erv+ sfac* erh);
      ezh= xx2*( salpj* ezv- sfac* erv);
      /* x,y,z fields for constant current */
      eph= sn* sph* xx2* eph;
      e[0]= erh* rhx+ eph* phx;
      e[1]= erh* rhy+ eph* phy;
      e[2]= ezh;
      /* x,y,z fields for sine current */
      rk= TP* t;
      sfac= sinl( rk);
      e[3]= e[0]* sfac;
      e[4]= e[1]* sfac;
      /* x,y,z fields for cosine current */
      e[5]= e[2]* sfac;
      sfac= cosl( rk);
      e[6]= e[0]* sfac;
      e[7]= e[1]* sfac;
      e[8]= e[2]* sfac;

      return;
}
/*-----------------------------------------------------------------------*/


Generated by  Doxygen 1.6.0   Back to index