Subversion Repositories Nec2c

Rev

Blame | Last modification | View Log | Download | RSS feed

  1. /*** Translated to the C language by N. Kyriazis  20 Aug 2003 ***
  2.  
  3.   Program NEC(input,tape5=input,output,tape11,tape12,tape13,tape14,
  4.   tape15,tape16,tape20,tape21)
  5.  
  6.   Numerical Electromagnetics Code (NEC2)  developed at Lawrence
  7.   Livermore lab., Livermore, CA.  (contact G. Burke at 415-422-8414
  8.   for problems with the NEC code. For problems with the vax implem-
  9.   entation, contact J. Breakall at 415-422-8196 or E. Domning at 415
  10.   422-5936)
  11.   file created 4/11/80.
  12.  
  13.                                 ***********Notice**********
  14.  This computer code material was prepared as an account of work
  15.  sponsored by the United States government.  Neither the United
  16.  States nor the United States Department Of Energy, nor any of
  17.  their employees, nor any of their contractors, subcontractors,
  18.  or their employees, makes any warranty, express or implied, or
  19.  assumes any legal liability or responsibility for the accuracy,
  20.  completeness or usefulness of any information, apparatus, product
  21.  or process disclosed, or represents that its use would not infringe
  22.  privately-owned rights.
  23.  
  24.  *******************************************************************/
  25.  
  26. #include "nec2c.h"
  27.  
  28. /* pointers to input/output files */
  29. extern FILE *input_fp, *output_fp, *plot_fp;
  30.  
  31. /* common  /data/ */
  32. extern data_t data;
  33.  
  34. /* common  /dataj/ */
  35. extern dataj_t dataj;
  36.  
  37. /* common  /incom/ */
  38. incom_t incom;
  39. ;
  40.  
  41. /* common  /gwav/ */
  42. extern gwav_t gwav;
  43.  
  44. /* common  /gnd/ */
  45. extern gnd_t gnd;
  46.  
  47. /*-------------------------------------------------------------------*/
  48.  
  49. /* segment to obtain the total field due to ground.  the method of */
  50. /* variable interval width romberg integration is used.  there are 9 */
  51. /* field components - the x, y, and z components due to constant, */
  52. /* sine, and cosine current distributions. */
  53. void rom2 (double a, double b, complex double *sum, double dmin)
  54. {
  55.         int i, ns, nt, flag = TRUE;
  56.         int nts = 4, nx = 1, n = 9;
  57.         double ze, ep, zend, dz = 0., dzot = 0., tmag1, tmag2, tr, ti;
  58.         double z, s; /***also global***/
  59.         double rx = 1.0e-4;
  60.         complex double g1[9], g2[9], g3[9], g4[9], g5[9];
  61.         complex double t00, t01[9], t10[9], t02, t11, t20[9];
  62.  
  63.         z = a;
  64.         ze = b;
  65.         s = b - a;
  66.  
  67.         if (s < 0.)
  68.         {
  69.                 fprintf (output_fp, "\n  ERROR - B LESS THAN A IN ROM2");
  70.                 stop (-1);
  71.         }
  72.  
  73.         ep = s / (1.e4 * data.npm);
  74.         zend = ze - ep;
  75.  
  76.         for (i = 0; i < n; i++)
  77.                 sum[i] = CPLX_00;
  78.  
  79.         ns = nx;
  80.         nt = 0;
  81.         sflds (z, g1);
  82.  
  83.         while (TRUE)
  84.         {
  85.                 if (flag)
  86.                 {
  87.                         dz = s / ns;
  88.                         if (z + dz > ze)
  89.                         {
  90.                                 dz = ze - z;
  91.                                 if (dz <= ep)
  92.                                         return;
  93.                         }
  94.  
  95.                         dzot = dz * .5;
  96.                         sflds (z + dzot, g3);
  97.                         sflds (z + dz, g5);
  98.  
  99.                 } /* if( flag ) */
  100.  
  101.                 tmag1 = 0.;
  102.                 tmag2 = 0.;
  103.  
  104.                 /* evaluate 3 point romberg result and test convergence. */
  105.                 for (i = 0; i < n; i++)
  106.                 {
  107.                         t00 = (g1[i] + g5[i]) * dzot;
  108.                         t01[i] = (t00 + dz * g3[i]) * .5;
  109.                         t10[i] = (4. * t01[i] - t00) / 3.;
  110.                         if (i > 2)
  111.                                 continue;
  112.  
  113.                         tr = creal (t01[i]);
  114.                         ti = cimag (t01[i]);
  115.                         tmag1 = tmag1 + tr * tr + ti * ti;
  116.                         tr = creal (t10[i]);
  117.                         ti = cimag (t10[i]);
  118.                         tmag2 = tmag2 + tr * tr + ti * ti;
  119.  
  120.                 } /* for( i = 0; i < n; i++ ) */
  121.  
  122.                 tmag1 = sqrtl (tmag1);
  123.                 tmag2 = sqrtl (tmag2);
  124.                 test (tmag1, tmag2, &tr, 0., 0., &ti, dmin);
  125.  
  126.                 if (tr <= rx)
  127.                 {
  128.                         for (i = 0; i < n; i++)
  129.                                 sum[i] += t10[i];
  130.                         nt += 2;
  131.  
  132.                         z += dz;
  133.                         if (z > zend)
  134.                                 return;
  135.  
  136.                         for (i = 0; i < n; i++)
  137.                                 g1[i] = g5[i];
  138.  
  139.                         if ((nt >= nts) && (ns > nx))
  140.                         {
  141.                                 ns = ns / 2;
  142.                                 nt = 1;
  143.                         }
  144.                         flag = TRUE;
  145.                         continue;
  146.  
  147.                 } /* if( tr <= rx) */
  148.  
  149.                 sflds (z + dz * .25, g2);
  150.                 sflds (z + dz * .75, g4);
  151.                 tmag1 = 0.;
  152.                 tmag2 = 0.;
  153.  
  154.                 /* evaluate 5 point romberg result and test convergence. */
  155.                 for (i = 0; i < n; i++)
  156.                 {
  157.                         t02 = (t01[i] + dzot * (g2[i] + g4[i])) * .5;
  158.                         t11 = (4.0 * t02 - t01[i]) / 3.;
  159.                         t20[i] = (16. * t11 - t10[i]) / 15.;
  160.                         if (i > 2)
  161.                                 continue;
  162.  
  163.                         tr = creal (t11);
  164.                         ti = cimag (t11);
  165.                         tmag1 = tmag1 + tr * tr + ti * ti;
  166.                         tr = creal (t20[i]);
  167.                         ti = cimag (t20[i]);
  168.                         tmag2 = tmag2 + tr * tr + ti * ti;
  169.  
  170.                 } /* for( i = 0; i < n; i++ ) */
  171.  
  172.                 tmag1 = sqrtl (tmag1);
  173.                 tmag2 = sqrtl (tmag2);
  174.                 test (tmag1, tmag2, &tr, 0., 0., &ti, dmin);
  175.  
  176.                 if (tr > rx)
  177.                 {
  178.                         nt = 0;
  179.                         if (ns < data.npm)
  180.                         {
  181.                                 ns = ns * 2;
  182.                                 dz = s / ns;
  183.                                 dzot = dz * .5;
  184.  
  185.                                 for (i = 0; i < n; i++)
  186.                                 {
  187.                                         g5[i] = g3[i];
  188.                                         g3[i] = g2[i];
  189.                                 }
  190.  
  191.                                 flag = FALSE;
  192.                                 continue;
  193.  
  194.                         } /* if( ns < npm) */
  195.  
  196.                         fprintf (output_fp, "\n  ROM2 -- STEP SIZE LIMITED AT Z = %12.5le", z);
  197.  
  198.                 } /* if( tr > rx) */
  199.  
  200.                 for (i = 0; i < n; i++)
  201.                         sum[i] = sum[i] + t20[i];
  202.                 nt = nt + 1;
  203.  
  204.                 z = z + dz;
  205.                 if (z > zend)
  206.                         return;
  207.  
  208.                 for (i = 0; i < n; i++)
  209.                         g1[i] = g5[i];
  210.  
  211.                 flag = TRUE;
  212.                 if ((nt < nts) || (ns <= nx))
  213.                         continue;
  214.  
  215.                 ns = ns / 2;
  216.                 nt = 1;
  217.  
  218.         } /* while( TRUE ) */
  219. }
  220.  
  221. /*-----------------------------------------------------------------------*/
  222.  
  223. /* sfldx returns the field due to ground for a current element on */
  224. /* the source segment at t relative to the segment center. */
  225. void sflds (double t, complex double *e)
  226. {
  227.         double xt, yt, zt, rhx, rhy, rhs, rho, phx, phy;
  228.         double cph, sph, zphs, r2s, rk, sfac, thet;
  229.         complex double erv, ezv, erh, ezh, eph, er, et, hrv, hzv, hrh;
  230.  
  231.         xt = dataj.xj + t * dataj.cabj;
  232.         yt = dataj.yj + t * dataj.sabj;
  233.         zt = dataj.zj + t * dataj.salpj;
  234.         rhx = incom.xo - xt;
  235.         rhy = incom.yo - yt;
  236.         rhs = rhx * rhx + rhy * rhy;
  237.         rho = sqrtl (rhs);
  238.  
  239.         if (rho <= 0.)
  240.         {
  241.                 rhx = 1.;
  242.                 rhy = 0.;
  243.                 phx = 0.;
  244.                 phy = 1.;
  245.         }
  246.         else
  247.         {
  248.                 rhx = rhx / rho;
  249.                 rhy = rhy / rho;
  250.                 phx = -rhy;
  251.                 phy = rhx;
  252.         }
  253.  
  254.         cph = rhx * incom.xsn + rhy * incom.ysn;
  255.         sph = rhy * incom.xsn - rhx * incom.ysn;
  256.  
  257.         if (fabsl (cph) < 1.0e-10)
  258.                 cph = 0.;
  259.         if (fabsl (sph) < 1.0e-10)
  260.                 sph = 0.;
  261.  
  262.         gwav.zph = incom.zo + zt;
  263.         zphs = gwav.zph * gwav.zph;
  264.         r2s = rhs + zphs;
  265.         gwav.r2 = sqrtl (r2s);
  266.         rk = gwav.r2 * TP;
  267.         gwav.xx2 = cmplx (cosl (rk), -sinl (rk));
  268.  
  269.         /* use norton approximation for field due to ground.  current is */
  270.         /* lumped at segment center with current moment for constant, sine, */
  271.         /* or cosine distribution. */
  272.         if (incom.isnor != 1)
  273.         {
  274.                 gwav.zmh = 1.;
  275.                 gwav.r1 = 1.;
  276.                 gwav.xx1 = 0.;
  277.                 gwave (&erv, &ezv, &erh, &ezh, &eph);
  278.  
  279.                 et = -CONST1 * gnd.frati * gwav.xx2 / (r2s * gwav.r2);
  280.                 er = 2. * et * cmplx (1.0, rk);
  281.                 et = et * cmplx (1.0 - rk * rk, rk);
  282.                 hrv = (er + et) * rho * gwav.zph / r2s;
  283.                 hzv = (zphs * er - rhs * et) / r2s;
  284.                 hrh = (rhs * er - zphs * et) / r2s;
  285.                 erv = erv - hrv;
  286.                 ezv = ezv - hzv;
  287.                 erh = erh + hrh;
  288.                 ezh = ezh + hrv;
  289.                 eph = eph + et;
  290.                 erv = erv * dataj.salpj;
  291.                 ezv = ezv * dataj.salpj;
  292.                 erh = erh * incom.sn * cph;
  293.                 ezh = ezh * incom.sn * cph;
  294.                 eph = eph * incom.sn * sph;
  295.                 erh = erv + erh;
  296.                 e[0] = (erh * rhx + eph * phx) * dataj.s;
  297.                 e[1] = (erh * rhy + eph * phy) * dataj.s;
  298.                 e[2] = (ezv + ezh) * dataj.s;
  299.                 e[3] = 0.;
  300.                 e[4] = 0.;
  301.                 e[5] = 0.;
  302.                 sfac = PI * dataj.s;
  303.                 sfac = sinl (sfac) / sfac;
  304.                 e[6] = e[0] * sfac;
  305.                 e[7] = e[1] * sfac;
  306.                 e[8] = e[2] * sfac;
  307.  
  308.                 return;
  309.         } /* if( smat.isnor != 1) */
  310.  
  311.         /* interpolate in sommerfeld field tables */
  312.         if (rho >= 1.0e-12)
  313.                 thet = atanl (gwav.zph / rho);
  314.         else
  315.                 thet = POT;
  316.  
  317.         /* combine vertical and horizontal components and convert */
  318.         /* to x,y,z components. multiply by exp(-jkr)/r. */
  319.         intrp (gwav.r2, thet, &erv, &ezv, &erh, &eph);
  320.         gwav.xx2 = gwav.xx2 / gwav.r2;
  321.         sfac = incom.sn * cph;
  322.         erh = gwav.xx2 * (dataj.salpj * erv + sfac * erh);
  323.         ezh = gwav.xx2 * (dataj.salpj * ezv - sfac * erv);
  324.         /* x,y,z fields for constant current */
  325.         eph = incom.sn * sph * gwav.xx2 * eph;
  326.         e[0] = erh * rhx + eph * phx;
  327.         e[1] = erh * rhy + eph * phy;
  328.         e[2] = ezh;
  329.         /* x,y,z fields for sine current */
  330.         rk = TP * t;
  331.         sfac = sinl (rk);
  332.         e[3] = e[0] * sfac;
  333.         e[4] = e[1] * sfac;
  334.         /* x,y,z fields for cosine current */
  335.         e[5] = e[2] * sfac;
  336.         sfac = cosl (rk);
  337.         e[6] = e[0] * sfac;
  338.         e[7] = e[1] * sfac;
  339.         e[8] = e[2] * sfac;
  340.  
  341.         return;
  342. }
  343.  
  344. /*-----------------------------------------------------------------------*/
  345.