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. /* common  /data/ */
  29. extern data_t data;
  30.  
  31. /* common  /gnd/ */
  32. extern gnd_t gnd;
  33.  
  34. /* common  /crnt/ */
  35. extern crnt_t crnt;
  36.  
  37. /* common  /gwav/ */
  38. extern gwav_t gwav;
  39.  
  40. /* common  /fpat/ */
  41. extern fpat_t fpat;
  42.  
  43. /* pointers to input/output files */
  44. extern FILE *input_fp, *output_fp, *plot_fp;
  45.  
  46. /* common  /save/ */
  47. extern save_t save;
  48.  
  49. /* common  /plot/ */
  50. extern plot_t plot;
  51.  
  52. /*-----------------------------------------------------------------------*/
  53.  
  54. /* ffld calculates the far zone radiated electric fields, */
  55. /* the factor exp(j*k*r)/(r/lamda) not included */
  56. void ffld (double thet, double phi, complex double *eth, complex double *eph)
  57. {
  58.         int k, i, ip, jump;
  59.         double phx, phy, roz, rozs, thx, thy, thz, rox, roy;
  60.         double tthet = 0., darg = 0., omega, el, sill, top, bot, a;
  61.         double too, boo, b, c, d, rr, ri, arg, dr, rfl, rrz;
  62.         complex double cix = CPLX_00, ciy = CPLX_00, ciz = CPLX_00;
  63.         complex double exa, ccx = CPLX_00, ccy = CPLX_00, ccz = CPLX_00, cdp;
  64.         complex double zrsin, rrv = CPLX_00, rrh = CPLX_00, rrv1 = CPLX_00;
  65.         complex double rrh1 = CPLX_00, rrv2 = CPLX_00, rrh2 = CPLX_00;
  66.         complex double tix, tiy, tiz, zscrn, ex = CPLX_00, ey = CPLX_00, ez = CPLX_00, gx, gy,
  67.                                              gz;
  68.  
  69.         phx = -sinl (phi);
  70.         phy = cosl (phi);
  71.         roz = cosl (thet);
  72.         rozs = roz;
  73.         thx = roz * phy;
  74.         thy = -roz * phx;
  75.         thz = -sinl (thet);
  76.         rox = -thz * phy;
  77.         roy = thz * phx;
  78.  
  79.         jump = FALSE;
  80.         if (data.n != 0)
  81.         {
  82.                 /* loop for structure image if any */
  83.                 /* calculation of reflection coeffecients */
  84.                 for (k = 0; k < gnd.ksymp; k++)
  85.                 {
  86.                         if (k != 0)
  87.                         {
  88.                                 /* for perfect ground */
  89.                                 if (gnd.iperf == 1)
  90.                                 {
  91.                                         rrv = -CPLX_10;
  92.                                         rrh = -CPLX_10;
  93.                                 }
  94.                                 else
  95.                                 {
  96.                                         /* for infinite planar ground */
  97.                                         zrsin =
  98.                                             csqrtl (1. - gnd.zrati * gnd.zrati * thz * thz);
  99.                                         rrv = -(roz - gnd.zrati * zrsin) /
  100.                                               (roz + gnd.zrati * zrsin);
  101.                                         rrh = (gnd.zrati * roz - zrsin) /
  102.                                               (gnd.zrati * roz + zrsin);
  103.  
  104.                                 } /* if( gnd.iperf == 1) */
  105.  
  106.                                 /* for the cliff problem, two reflction coefficients calculated
  107.                                  */
  108.                                 if (gnd.ifar > 1)
  109.                                 {
  110.                                         rrv1 = rrv;
  111.                                         rrh1 = rrh;
  112.                                         tthet = tanl (thet);
  113.  
  114.                                         if (gnd.ifar != 4)
  115.                                         {
  116.                                                 zrsin = csqrtl (
  117.                                                     1. - gnd.zrati2 * gnd.zrati2 * thz * thz);
  118.                                                 rrv2 = -(roz - gnd.zrati2 * zrsin) /
  119.                                                        (roz + gnd.zrati2 * zrsin);
  120.                                                 rrh2 = (gnd.zrati2 * roz - zrsin) /
  121.                                                        (gnd.zrati2 * roz + zrsin);
  122.                                                 darg = -TP * 2. * gnd.ch * roz;
  123.                                         }
  124.                                 } /* if( gnd.ifar > 1) */
  125.  
  126.                                 roz = -roz;
  127.                                 ccx = cix;
  128.                                 ccy = ciy;
  129.                                 ccz = ciz;
  130.  
  131.                         } /* if( k != 0 ) */
  132.  
  133.                         cix = CPLX_00;
  134.                         ciy = CPLX_00;
  135.                         ciz = CPLX_00;
  136.  
  137.                         /* loop over structure segments */
  138.                         for (i = 0; i < data.n; i++)
  139.                         {
  140.                                 omega =
  141.                                     -(rox * data.cab[i] + roy * data.sab[i] +
  142.                                       roz * data.salp[i]);
  143.                                 el = PI * data.si[i];
  144.                                 sill = omega * el;
  145.                                 top = el + sill;
  146.                                 bot = el - sill;
  147.  
  148.                                 if (fabsl (omega) >= 1.0e-7)
  149.                                         a = 2. * sinl (sill) / omega;
  150.                                 else
  151.                                         a = (2. - omega * omega * el * el / 3.) * el;
  152.  
  153.                                 if (fabsl (top) >= 1.0e-7)
  154.                                         too = sinl (top) / top;
  155.                                 else
  156.                                         too = 1. - top * top / 6.;
  157.  
  158.                                 if (fabsl (bot) >= 1.0e-7)
  159.                                         boo = sinl (bot) / bot;
  160.                                 else
  161.                                         boo = 1. - bot * bot / 6.;
  162.  
  163.                                 b = el * (boo - too);
  164.                                 c = el * (boo + too);
  165.                                 rr = a * crnt.air[i] + b * crnt.bii[i] + c * crnt.cir[i];
  166.                                 ri = a * crnt.aii[i] - b * crnt.bir[i] + c * crnt.cii[i];
  167.                                 arg =
  168.                                     TP * (data.x[i] * rox + data.y[i] * roy + data.z[i] * roz);
  169.  
  170.                                 if ((k != 1) || (gnd.ifar < 2))
  171.                                 {
  172.                                         /* summation for far field integral */
  173.                                         exa = cmplx (cosl (arg), sinl (arg)) * cmplx (rr, ri);
  174.                                         cix = cix + exa * data.cab[i];
  175.                                         ciy = ciy + exa * data.sab[i];
  176.                                         ciz = ciz + exa * data.salp[i];
  177.                                         continue;
  178.                                 }
  179.  
  180.                                 /* calculation of image contribution */
  181.                                 /* in cliff and ground screen problems */
  182.  
  183.                                 /* specular point distance */
  184.                                 dr = data.z[i] * tthet;
  185.  
  186.                                 d = dr * phy + data.x[i];
  187.                                 if (gnd.ifar == 2)
  188.                                 {
  189.                                         if ((gnd.cl - d) > 0.)
  190.                                         {
  191.                                                 rrv = rrv1;
  192.                                                 rrh = rrh1;
  193.                                         }
  194.                                         else
  195.                                         {
  196.                                                 rrv = rrv2;
  197.                                                 rrh = rrh2;
  198.                                                 arg = arg + darg;
  199.                                         }
  200.                                 } /* if( gnd.ifar == 2) */
  201.                                 else
  202.                                 {
  203.                                         d = sqrtl (
  204.                                             d * d +
  205.                                             (data.y[i] - dr * phx) * (data.y[i] - dr * phx));
  206.                                         if (gnd.ifar == 3)
  207.                                         {
  208.                                                 if ((gnd.cl - d) > 0.)
  209.                                                 {
  210.                                                         rrv = rrv1;
  211.                                                         rrh = rrh1;
  212.                                                 }
  213.                                                 else
  214.                                                 {
  215.                                                         rrv = rrv2;
  216.                                                         rrh = rrh2;
  217.                                                         arg = arg + darg;
  218.                                                 }
  219.                                         } /* if( gnd.ifar == 3) */
  220.                                         else
  221.                                         {
  222.                                                 if ((gnd.scrwl - d) >= 0.)
  223.                                                 {
  224.                                                         /* radial wire ground screen reflection
  225.                                                          * coefficient */
  226.                                                         d = d + gnd.t2;
  227.                                                         zscrn = gnd.t1 * d * logl (d / gnd.t2);
  228.                                                         zscrn = (zscrn * gnd.zrati) /
  229.                                                                 (ETA * gnd.zrati + zscrn);
  230.                                                         zrsin = csqrtl (
  231.                                                             1. - zscrn * zscrn * thz * thz);
  232.                                                         rrv = (roz + zscrn * zrsin) /
  233.                                                               (-roz + zscrn * zrsin);
  234.                                                         rrh = (zscrn * roz + zrsin) /
  235.                                                               (zscrn * roz - zrsin);
  236.                                                 } /* if(( gnd.scrwl- d) < 0.) */
  237.                                                 else
  238.                                                 {
  239.                                                         if (gnd.ifar == 4)
  240.                                                         {
  241.                                                                 rrv = rrv1;
  242.                                                                 rrh = rrh1;
  243.                                                         } /* if( gnd.ifar == 4) */
  244.                                                         else
  245.                                                         {
  246.                                                                 if (gnd.ifar == 5)
  247.                                                                         d = dr * phy +
  248.                                                                             data.x[i];
  249.  
  250.                                                                 if ((gnd.cl - d) > 0.)
  251.                                                                 {
  252.                                                                         rrv = rrv1;
  253.                                                                         rrh = rrh1;
  254.                                                                 }
  255.                                                                 else
  256.                                                                 {
  257.                                                                         rrv = rrv2;
  258.                                                                         rrh = rrh2;
  259.                                                                         arg = arg + darg;
  260.                                                                 } /* if(( gnd.cl- d) > 0.) */
  261.  
  262.                                                         } /* if( gnd.ifar == 4) */
  263.  
  264.                                                 } /* if(( gnd.scrwl- d) < 0.) */
  265.  
  266.                                         } /* if( gnd.ifar == 3) */
  267.  
  268.                                 } /* if( gnd.ifar == 2) */
  269.  
  270.                                 /* contribution of each image segment modified by */
  271.                                 /* reflection coef, for cliff and ground screen problems */
  272.                                 exa = cmplx (cosl (arg), sinl (arg)) * cmplx (rr, ri);
  273.                                 tix = exa * data.cab[i];
  274.                                 tiy = exa * data.sab[i];
  275.                                 tiz = exa * data.salp[i];
  276.                                 cdp = (tix * phx + tiy * phy) * (rrh - rrv);
  277.                                 cix = cix + tix * rrv + cdp * phx;
  278.                                 ciy = ciy + tiy * rrv + cdp * phy;
  279.                                 ciz = ciz - tiz * rrv;
  280.  
  281.                         } /* for( i = 0; i < n; i++ ) */
  282.  
  283.                         if (k == 0)
  284.                                 continue;
  285.  
  286.                         /* calculation of contribution of structure image for infinite ground
  287.                          */
  288.                         if (gnd.ifar < 2)
  289.                         {
  290.                                 cdp = (cix * phx + ciy * phy) * (rrh - rrv);
  291.                                 cix = ccx + cix * rrv + cdp * phx;
  292.                                 ciy = ccy + ciy * rrv + cdp * phy;
  293.                                 ciz = ccz - ciz * rrv;
  294.                         }
  295.                         else
  296.                         {
  297.                                 cix = cix + ccx;
  298.                                 ciy = ciy + ccy;
  299.                                 ciz = ciz + ccz;
  300.                         }
  301.  
  302.                 } /* for( k=0; k < gnd.ksymp; k++ ) */
  303.  
  304.                 if (data.m > 0)
  305.                         jump = TRUE;
  306.                 else
  307.                 {
  308.                         *eth = (cix * thx + ciy * thy + ciz * thz) * CONST3;
  309.                         *eph = (cix * phx + ciy * phy) * CONST3;
  310.                         return;
  311.                 }
  312.  
  313.         } /* if( n != 0) */
  314.  
  315.         if (!jump)
  316.         {
  317.                 cix = CPLX_00;
  318.                 ciy = CPLX_00;
  319.                 ciz = CPLX_00;
  320.         }
  321.  
  322.         /* electric field components */
  323.         roz = rozs;
  324.         rfl = -1.;
  325.         for (ip = 0; ip < gnd.ksymp; ip++)
  326.         {
  327.                 rfl = -rfl;
  328.                 rrz = roz * rfl;
  329.                 fflds (rox, roy, rrz, &crnt.cur[data.n], &gx, &gy, &gz);
  330.  
  331.                 if (ip != 1)
  332.                 {
  333.                         ex = gx;
  334.                         ey = gy;
  335.                         ez = gz;
  336.                         continue;
  337.                 }
  338.  
  339.                 if (gnd.iperf == 1)
  340.                 {
  341.                         gx = -gx;
  342.                         gy = -gy;
  343.                         gz = -gz;
  344.                 }
  345.                 else
  346.                 {
  347.                         rrv = csqrtl (1. - gnd.zrati * gnd.zrati * thz * thz);
  348.                         rrh = gnd.zrati * roz;
  349.                         rrh = (rrh - rrv) / (rrh + rrv);
  350.                         rrv = gnd.zrati * rrv;
  351.                         rrv = -(roz - rrv) / (roz + rrv);
  352.                         *eth = (gx * phx + gy * phy) * (rrh - rrv);
  353.                         gx = gx * rrv + *eth * phx;
  354.                         gy = gy * rrv + *eth * phy;
  355.                         gz = gz * rrv;
  356.  
  357.                 } /* if( gnd.iperf == 1) */
  358.  
  359.                 ex = ex + gx;
  360.                 ey = ey + gy;
  361.                 ez = ez - gz;
  362.  
  363.         } /* for( ip = 0; ip < gnd.ksymp; ip++ ) */
  364.  
  365.         ex = ex + cix * CONST3;
  366.         ey = ey + ciy * CONST3;
  367.         ez = ez + ciz * CONST3;
  368.         *eth = ex * thx + ey * thy + ez * thz;
  369.         *eph = ex * phx + ey * phy;
  370.  
  371.         return;
  372. }
  373.  
  374. /*-----------------------------------------------------------------------*/
  375.  
  376. /* calculates the xyz components of the electric */
  377. /* field due to surface currents */
  378. void fflds (
  379.     double rox,
  380.     double roy,
  381.     double roz,
  382.     complex double *scur,
  383.     complex double *ex,
  384.     complex double *ey,
  385.     complex double *ez)
  386. {
  387.         double *xs, *ys, *zs, *s;
  388.         int j, i, k;
  389.         double arg;
  390.         complex double ct;
  391.  
  392.         xs = data.px;
  393.         ys = data.py;
  394.         zs = data.pz;
  395.         s = data.pbi;
  396.  
  397.         *ex = CPLX_00;
  398.         *ey = CPLX_00;
  399.         *ez = CPLX_00;
  400.  
  401.         i = -1;
  402.         for (j = 0; j < data.m; j++)
  403.         {
  404.                 i++;
  405.                 arg = TP * (rox * xs[i] + roy * ys[i] + roz * zs[i]);
  406.                 ct = cmplx (cosl (arg) * s[i], sinl (arg) * s[i]);
  407.                 k = 3 * j;
  408.                 *ex += scur[k] * ct;
  409.                 *ey += scur[k + 1] * ct;
  410.                 *ez += scur[k + 2] * ct;
  411.         }
  412.  
  413.         ct = rox * *ex + roy * *ey + roz * *ez;
  414.         *ex = CONST4 * (ct * rox - *ex);
  415.         *ey = CONST4 * (ct * roy - *ey);
  416.         *ez = CONST4 * (ct * roz - *ez);
  417.  
  418.         return;
  419. }
  420.  
  421. /*-----------------------------------------------------------------------*/
  422.  
  423. /* gfld computes the radiated field including ground wave. */
  424. void gfld (
  425.     double rho,
  426.     double phi,
  427.     double rz,
  428.     complex double *eth,
  429.     complex double *epi,
  430.     complex double *erd,
  431.     complex double ux,
  432.     int ksymp)
  433. {
  434.         int i, k;
  435.         double b, r, thet, arg, phx, phy, rx, ry, dx, dy, dz, rix, riy, rhs, rhp;
  436.         double rhx, rhy, calp, cbet, sbet, cph, sph, el, rfl, riz, thx, thy, thz;
  437.         double rxyz, rnx, rny, rnz, omega, sill, top, bot, a, too, boo, c, rr, ri;
  438.         complex double cix, ciy, ciz, exa, erv;
  439.         complex double ezv, erh, eph, ezh, ex, ey;
  440.  
  441.         r = sqrtl (rho * rho + rz * rz);
  442.         if ((ksymp == 1) || (cabs (ux) > .5) || (r > 1.e5))
  443.         {
  444.                 /* computation of space wave only */
  445.                 if (rz >= 1.0e-20)
  446.                         thet = atanl (rho / rz);
  447.                 else
  448.                         thet = PI * .5;
  449.  
  450.                 ffld (thet, phi, eth, epi);
  451.                 arg = -TP * r;
  452.                 exa = cmplx (cosl (arg), sinl (arg)) / r;
  453.                 *eth = *eth * exa;
  454.                 *epi = *epi * exa;
  455.                 *erd = CPLX_00;
  456.                 return;
  457.         } /* if( (ksymp == 1) && (cabs(ux) > .5) && (r > 1.e5) ) */
  458.  
  459.         /* computation of space and ground waves. */
  460.         gwav.u = ux;
  461.         gwav.u2 = gwav.u * gwav.u;
  462.         phx = -sinl (phi);
  463.         phy = cosl (phi);
  464.         rx = rho * phy;
  465.         ry = -rho * phx;
  466.         cix = CPLX_00;
  467.         ciy = CPLX_00;
  468.         ciz = CPLX_00;
  469.  
  470.         /* summation of field from individual segments */
  471.         for (i = 0; i < data.n; i++)
  472.         {
  473.                 dx = data.cab[i];
  474.                 dy = data.sab[i];
  475.                 dz = data.salp[i];
  476.                 rix = rx - data.x[i];
  477.                 riy = ry - data.y[i];
  478.                 rhs = rix * rix + riy * riy;
  479.                 rhp = sqrtl (rhs);
  480.  
  481.                 if (rhp >= 1.0e-6)
  482.                 {
  483.                         rhx = rix / rhp;
  484.                         rhy = riy / rhp;
  485.                 }
  486.                 else
  487.                 {
  488.                         rhx = 1.;
  489.                         rhy = 0.;
  490.                 }
  491.  
  492.                 calp = 1. - dz * dz;
  493.                 if (calp >= 1.0e-6)
  494.                 {
  495.                         calp = sqrtl (calp);
  496.                         cbet = dx / calp;
  497.                         sbet = dy / calp;
  498.                         cph = rhx * cbet + rhy * sbet;
  499.                         sph = rhy * cbet - rhx * sbet;
  500.                 }
  501.                 else
  502.                 {
  503.                         cph = rhx;
  504.                         sph = rhy;
  505.                 }
  506.  
  507.                 el = PI * data.si[i];
  508.                 rfl = -1.;
  509.  
  510.                 /* integration of (current)*(phase factor) over segment and image for */
  511.                 /* constant, sine, and cosine current distributions */
  512.                 for (k = 0; k < 2; k++)
  513.                 {
  514.                         rfl = -rfl;
  515.                         riz = rz - data.z[i] * rfl;
  516.                         rxyz = sqrtl (rix * rix + riy * riy + riz * riz);
  517.                         rnx = rix / rxyz;
  518.                         rny = riy / rxyz;
  519.                         rnz = riz / rxyz;
  520.                         omega = -(rnx * dx + rny * dy + rnz * dz * rfl);
  521.                         sill = omega * el;
  522.                         top = el + sill;
  523.                         bot = el - sill;
  524.  
  525.                         if (fabsl (omega) >= 1.0e-7)
  526.                                 a = 2. * sinl (sill) / omega;
  527.                         else
  528.                                 a = (2. - omega * omega * el * el / 3.) * el;
  529.  
  530.                         if (fabsl (top) >= 1.0e-7)
  531.                                 too = sinl (top) / top;
  532.                         else
  533.                                 too = 1. - top * top / 6.;
  534.  
  535.                         if (fabsl (bot) >= 1.0e-7)
  536.                                 boo = sinl (bot) / bot;
  537.                         else
  538.                                 boo = 1. - bot * bot / 6.;
  539.  
  540.                         b = el * (boo - too);
  541.                         c = el * (boo + too);
  542.                         rr = a * crnt.air[i] + b * crnt.bii[i] + c * crnt.cir[i];
  543.                         ri = a * crnt.aii[i] - b * crnt.bir[i] + c * crnt.cii[i];
  544.                         arg = TP * (data.x[i] * rnx + data.y[i] * rny + data.z[i] * rnz * rfl);
  545.                         exa = cmplx (cosl (arg), sinl (arg)) * cmplx (rr, ri) / TP;
  546.  
  547.                         if (k != 1)
  548.                         {
  549.                                 gwav.xx1 = exa;
  550.                                 gwav.r1 = rxyz;
  551.                                 gwav.zmh = riz;
  552.                                 continue;
  553.                         }
  554.  
  555.                         gwav.xx2 = exa;
  556.                         gwav.r2 = rxyz;
  557.                         gwav.zph = riz;
  558.  
  559.                 } /* for( k = 0; k < 2; k++ ) */
  560.  
  561.                 /* call subroutine to compute the field */
  562.                 /* of segment including ground wave. */
  563.                 gwave (&erv, &ezv, &erh, &ezh, &eph);
  564.                 erh = erh * cph * calp + erv * dz;
  565.                 eph = eph * sph * calp;
  566.                 ezh = ezh * cph * calp + ezv * dz;
  567.                 ex = erh * rhx - eph * rhy;
  568.                 ey = erh * rhy + eph * rhx;
  569.                 cix = cix + ex;
  570.                 ciy = ciy + ey;
  571.                 ciz = ciz + ezh;
  572.  
  573.         } /* for( i = 0; i < n; i++ ) */
  574.  
  575.         arg = -TP * r;
  576.         exa = cmplx (cosl (arg), sinl (arg));
  577.         cix = cix * exa;
  578.         ciy = ciy * exa;
  579.         ciz = ciz * exa;
  580.         rnx = rx / r;
  581.         rny = ry / r;
  582.         rnz = rz / r;
  583.         thx = rnz * phy;
  584.         thy = -rnz * phx;
  585.         thz = -rho / r;
  586.         *eth = cix * thx + ciy * thy + ciz * thz;
  587.         *epi = cix * phx + ciy * phy;
  588.         *erd = cix * rnx + ciy * rny + ciz * rnz;
  589.  
  590.         return;
  591. }
  592.  
  593. /*-----------------------------------------------------------------------*/
  594.  
  595. /* compute radiation pattern, gain, normalized gain */
  596. void rdpat (void)
  597. {
  598.         char *hpol[3] = {"LINEAR", "RIGHT ", "LEFT  "};
  599.         char *igtp[2] = {"- POWER GAINS -", "- DIRECTIVE GAINS -"};
  600.         char *igax[4] = {" MAJOR", " MINOR", "VERT. ", "HOR.  "};
  601.         char *igntp[5] = {
  602.             " MAJOR AXIS", "  MINOR AXIS", "    VERTICAL", "  HORIZONTAL", "       TOTAL "};
  603.  
  604.         char *hclif = NULL, *isens;
  605.         int i, j, jump, itmp1, itmp2, kth, kph, itmp3, itmp4;
  606.         double exrm = 0., exra = 0., prad, gcon, gcop, gmax, pint, tmp1, tmp2;
  607.         double phi, pha, thet, tha, erdm = 0., erda = 0., ethm2, ethm, *gain = NULL;
  608.         double etha, ephm2, ephm, epha, tilta, emajr2, eminr2, axrat;
  609.         double dfaz, dfaz2, cdfaz, tstor1 = 0., tstor2, stilta, gnmj;
  610.         double gnmn, gnv, gnh, gtot, tmp3, tmp4, da, tmp5, tmp6;
  611.         complex double eth, eph, erd;
  612.  
  613.         /* Allocate memory to gain buffer */
  614.         if (fpat.inor > 0)
  615.                 mem_alloc ((void *) &gain, fpat.nth * fpat.nph * sizeof (double));
  616.  
  617.         if (gnd.ifar > 1)
  618.         {
  619.                 fprintf (
  620.                     output_fp,
  621.                     "\n\n\n\n"
  622.                     "                               "
  623.                     "- - - FAR FIELD GROUND PARAMETERS - - -\n");
  624.  
  625.                 jump = FALSE;
  626.                 if (gnd.ifar > 3)
  627.                 {
  628.                         fprintf (
  629.                             output_fp,
  630.                             "\n"
  631.                             "                               "
  632.                             "--- RADIAL WIRE GROUND SCREEN ---\n"
  633.                             "                               "
  634.                             "NUM OF WIRES= %d\n"
  635.                             "                               "
  636.                             "WIRE LENGTH= %8.2f METERS\n"
  637.                             "                               "
  638.                             "WIRE RADIUS= %10.3E METERS",
  639.                             gnd.nradl,
  640.                             save.scrwlt,
  641.                             save.scrwrt);
  642.  
  643.                         if (gnd.ifar == 4)
  644.                                 jump = TRUE;
  645.  
  646.                 } /* if( gnd.ifar > 3) */
  647.  
  648.                 if (!jump)
  649.                 {
  650.                         if ((gnd.ifar == 2) || (gnd.ifar == 5))
  651.                                 hclif = "LINEAR";
  652.                         if ((gnd.ifar == 3) || (gnd.ifar == 6))
  653.                                 hclif = "CIRCULAR";
  654.  
  655.                         gnd.cl = fpat.clt / data.wlam;
  656.                         gnd.ch = fpat.cht / data.wlam;
  657.                         gnd.zrati2 =
  658.                             csqrtl (1. / cmplx (fpat.epsr2, -fpat.sig2 * data.wlam * 59.96));
  659.  
  660.                         fprintf (
  661.                             output_fp,
  662.                             "\n"
  663.                             "                               "
  664.                             "- - %s CLIFF - -\n"
  665.                             "                               "
  666.                             "EDGE DISTANCE= %9.2f METERS\n"
  667.                             "                               "
  668.                             "       HEIGHT= %9.2f METERS\n"
  669.                             "                               "
  670.                             "--- SECOND MEDIUM ---\n"
  671.                             "                               "
  672.                             "RELATIVE DIELECTRIC CONST= %10.3f\n"
  673.                             "                               "
  674.                             "      GROUND CONDUCTIVITY= %10.3f MHOS",
  675.                             hclif,
  676.                             fpat.clt,
  677.                             fpat.cht,
  678.                             fpat.epsr2,
  679.                             fpat.sig2);
  680.  
  681.                 } /* if( ! jump ) */
  682.  
  683.         } /* if( gnd.ifar > 1) */
  684.  
  685.         if (gnd.ifar == 1)
  686.         {
  687.                 fprintf (
  688.                     output_fp,
  689.                     "\n\n\n"
  690.                     "                             "
  691.                     "------- RADIATED FIELDS NEAR GROUND --------\n\n"
  692.                     "    ------- LOCATION -------     --- E(THETA) ---    "
  693.                     " ---- E(PHI) ----    --- E(RADIAL) ---\n"
  694.                     "      RHO    PHI        Z           MAG    PHASE     "
  695.                     "    MAG    PHASE        MAG     PHASE \n"
  696.                     "    METERS   DEG.     METERS      VOLTS/M DEGREES   "
  697.                     "   VOLTS/M   DEG.     VOLTS/M    DEG. ");
  698.         }
  699.         else
  700.         {
  701.                 itmp1 = 2 * fpat.iax;
  702.                 itmp2 = itmp1 + 1;
  703.  
  704.                 fprintf (
  705.                     output_fp,
  706.                     "\n\n\n\n"
  707.                     "                                                "
  708.                     "- - - RADIATION PATTERNS - - -\n");
  709.  
  710.                 if (fpat.rfld >= 1.0e-20)
  711.                 {
  712.                         exrm = 1. / fpat.rfld;
  713.                         exra = fpat.rfld / data.wlam;
  714.                         exra = -360. * (exra - floorl (exra));
  715.  
  716.                         fprintf (
  717.                             output_fp,
  718.                             "\n"
  719.                             "                             "
  720.                             "RANGE: %13.6E METERS\n"
  721.                             "                             "
  722.                             "EXP(-JKR)/R: %12.5E AT PHASE: %7.2f DEGREES\n",
  723.                             fpat.rfld,
  724.                             exrm,
  725.                             exra);
  726.                 }
  727.  
  728.                 fprintf (
  729.                     output_fp,
  730.                     "\n"
  731.                     "  - - ANGLES - -   %23s       - - - POLARIZATION - - -  "
  732.                     "  - - - E(THETA) - - -    - - - E(PHI) - - -\n"
  733.                     "  THETA     PHI        %6s  %6s  TOTAL      AXIAL     TILT  "
  734.                     " SENSE     MAGNITUDE    PHASE      MAGNITUDE    PHASE \n"
  735.                     " DEGREES  DEGREES       DB      DB      DB        RATIO     "
  736.                     "DEG.              VOLTS/M    DEGREES      VOLTS/M    DEGREES",
  737.                     igtp[fpat.ipd],
  738.                     igax[itmp1],
  739.                     igax[itmp2]);
  740.  
  741.         } /* if( gnd.ifar == 1) */
  742.  
  743.         if ((fpat.ixtyp == 0) || (fpat.ixtyp == 5))
  744.         {
  745.                 gcop = data.wlam * data.wlam * 2. * PI / (376.73 * fpat.pinr);
  746.                 prad = fpat.pinr - fpat.ploss - fpat.pnlr;
  747.                 gcon = gcop;
  748.                 if (fpat.ipd != 0)
  749.                         gcon = gcon * fpat.pinr / prad;
  750.         }
  751.         else if (fpat.ixtyp == 4)
  752.         {
  753.                 fpat.pinr = 394.51 * fpat.xpr6 * fpat.xpr6 * data.wlam * data.wlam;
  754.                 gcop = data.wlam * data.wlam * 2. * PI / (376.73 * fpat.pinr);
  755.                 prad = fpat.pinr - fpat.ploss - fpat.pnlr;
  756.                 gcon = gcop;
  757.                 if (fpat.ipd != 0)
  758.                         gcon = gcon * fpat.pinr / prad;
  759.         }
  760.         else
  761.         {
  762.                 prad = 0.;
  763.                 gcon = 4. * PI / (1. + fpat.xpr6 * fpat.xpr6);
  764.                 gcop = gcon;
  765.         }
  766.  
  767.         i = 0;
  768.         gmax = -1.e+10;
  769.         pint = 0.;
  770.         tmp1 = fpat.dph * TA;
  771.         tmp2 = .5 * fpat.dth * TA;
  772.         phi = fpat.phis - fpat.dph;
  773.  
  774.         for (kph = 1; kph <= fpat.nph; kph++)
  775.         {
  776.                 phi += fpat.dph;
  777.                 pha = phi * TA;
  778.                 thet = fpat.thets - fpat.dth;
  779.  
  780.                 for (kth = 1; kth <= fpat.nth; kth++)
  781.                 {
  782.                         thet += fpat.dth;
  783.                         if ((gnd.ksymp == 2) && (thet > 90.01) && (gnd.ifar != 1))
  784.                                 continue;
  785.  
  786.                         tha = thet * TA;
  787.                         if (gnd.ifar != 1)
  788.                                 ffld (tha, pha, &eth, &eph);
  789.                         else
  790.                         {
  791.                                 gfld (
  792.                                     fpat.rfld / data.wlam,
  793.                                     pha,
  794.                                     thet / data.wlam,
  795.                                     &eth,
  796.                                     &eph,
  797.                                     &erd,
  798.                                     gnd.zrati,
  799.                                     gnd.ksymp);
  800.                                 erdm = cabs (erd);
  801.                                 erda = cang (erd);
  802.                         }
  803.  
  804.                         ethm2 = creal (eth * conjl (eth));
  805.                         ethm = sqrtl (ethm2);
  806.                         etha = cang (eth);
  807.                         ephm2 = creal (eph * conjl (eph));
  808.                         ephm = sqrtl (ephm2);
  809.                         epha = cang (eph);
  810.  
  811.                         /* elliptical polarization calc. */
  812.                         if (gnd.ifar != 1)
  813.                         {
  814.                                 if ((ethm2 <= 1.0e-20) && (ephm2 <= 1.0e-20))
  815.                                 {
  816.                                         tilta = 0.;
  817.                                         emajr2 = 0.;
  818.                                         eminr2 = 0.;
  819.                                         axrat = 0.;
  820.                                         isens = " ";
  821.                                 }
  822.                                 else
  823.                                 {
  824.                                         dfaz = epha - etha;
  825.                                         if (epha >= 0.)
  826.                                                 dfaz2 = dfaz - 360.;
  827.                                         else
  828.                                                 dfaz2 = dfaz + 360.;
  829.  
  830.                                         if (fabsl (dfaz) > fabsl (dfaz2))
  831.                                                 dfaz = dfaz2;
  832.  
  833.                                         cdfaz = cosl (dfaz * TA);
  834.                                         tstor1 = ethm2 - ephm2;
  835.                                         tstor2 = 2. * ephm * ethm * cdfaz;
  836.                                         tilta = .5 * atan2l (tstor2, tstor1);
  837.                                         stilta = sinl (tilta);
  838.                                         tstor1 = tstor1 * stilta * stilta;
  839.                                         tstor2 = tstor2 * stilta * cosl (tilta);
  840.                                         emajr2 = -tstor1 + tstor2 + ethm2;
  841.                                         eminr2 = tstor1 - tstor2 + ephm2;
  842.                                         if (eminr2 < 0.)
  843.                                                 eminr2 = 0.;
  844.  
  845.                                         axrat = sqrtl (eminr2 / emajr2);
  846.                                         tilta = tilta * TD;
  847.                                         if (axrat <= 1.0e-5)
  848.                                                 isens = hpol[0];
  849.                                         else if (dfaz <= 0.)
  850.                                                 isens = hpol[1];
  851.                                         else
  852.                                                 isens = hpol[2];
  853.  
  854.                                 } /* if( (ethm2 <= 1.0e-20) && (ephm2 <= 1.0e-20) ) */
  855.  
  856.                                 gnmj = db10 (gcon * emajr2);
  857.                                 gnmn = db10 (gcon * eminr2);
  858.                                 gnv = db10 (gcon * ethm2);
  859.                                 gnh = db10 (gcon * ephm2);
  860.                                 gtot = db10 (gcon * (ethm2 + ephm2));
  861.  
  862.                                 if (fpat.inor > 0)
  863.                                 {
  864.                                         i++;
  865.                                         switch (fpat.inor)
  866.                                         {
  867.                                         case 1:
  868.                                                 tstor1 = gnmj;
  869.                                                 break;
  870.  
  871.                                         case 2:
  872.                                                 tstor1 = gnmn;
  873.                                                 break;
  874.  
  875.                                         case 3:
  876.                                                 tstor1 = gnv;
  877.                                                 break;
  878.  
  879.                                         case 4:
  880.                                                 tstor1 = gnh;
  881.                                                 break;
  882.  
  883.                                         case 5:
  884.                                                 tstor1 = gtot;
  885.                                         }
  886.  
  887.                                         gain[i - 1] = tstor1;
  888.                                         if (tstor1 > gmax)
  889.                                                 gmax = tstor1;
  890.  
  891.                                 } /* if( fpat.inor > 0) */
  892.  
  893.                                 if (fpat.iavp != 0)
  894.                                 {
  895.                                         tstor1 = gcop * (ethm2 + ephm2);
  896.                                         tmp3 = tha - tmp2;
  897.                                         tmp4 = tha + tmp2;
  898.  
  899.                                         if (kth == 1)
  900.                                                 tmp3 = tha;
  901.                                         else if (kth == fpat.nth)
  902.                                                 tmp4 = tha;
  903.  
  904.                                         da = fabsl (tmp1 * (cosl (tmp3) - cosl (tmp4)));
  905.                                         if ((kph == 1) || (kph == fpat.nph))
  906.                                                 da *= .5;
  907.                                         pint += tstor1 * da;
  908.  
  909.                                         if (fpat.iavp == 2)
  910.                                                 continue;
  911.                                 }
  912.  
  913.                                 if (fpat.iax != 1)
  914.                                 {
  915.                                         tmp5 = gnmj;
  916.                                         tmp6 = gnmn;
  917.                                 }
  918.                                 else
  919.                                 {
  920.                                         tmp5 = gnv;
  921.                                         tmp6 = gnh;
  922.                                 }
  923.  
  924.                                 ethm = ethm * data.wlam;
  925.                                 ephm = ephm * data.wlam;
  926.  
  927.                                 if (fpat.rfld >= 1.0e-20)
  928.                                 {
  929.                                         ethm = ethm * exrm;
  930.                                         etha = etha + exra;
  931.                                         ephm = ephm * exrm;
  932.                                         epha = epha + exra;
  933.                                 }
  934.  
  935.                                 fprintf (
  936.                                     output_fp,
  937.                                     "\n"
  938.                                     "%8.2f%9.2f   %8.2f%8.2f%8.2f%11.5f"
  939.                                     " %8.2f  %5s    %11.5E%9.2f    %11.5E%9.2f",
  940.                                     thet,
  941.                                     phi,
  942.                                     tmp5,
  943.                                     tmp6,
  944.                                     gtot,
  945.                                     axrat,
  946.                                     tilta,
  947.                                     isens,
  948.                                     ethm,
  949.                                     etha,
  950.                                     ephm,
  951.                                     epha);
  952.  
  953.                                 if (plot.iplp1 != 3)
  954.                                         continue;
  955.  
  956.                                 if (plot.iplp3 != 0)
  957.                                 {
  958.                                         if (plot.iplp2 == 1)
  959.                                         {
  960.                                                 if (plot.iplp3 == 1)
  961.                                                         fprintf (
  962.                                                             plot_fp,
  963.                                                             "%12.5E %12.5E %12.5E\n",
  964.                                                             (float) thet,
  965.                                                             (float) ethm,
  966.                                                             (float) etha);
  967.                                                 else if (plot.iplp3 == 2)
  968.                                                         fprintf (
  969.                                                             plot_fp,
  970.                                                             "%12.5E %12.5E %12.5E\n",
  971.                                                             (float) thet,
  972.                                                             (float) ephm,
  973.                                                             (float) epha);
  974.                                         }
  975.  
  976.                                         if (plot.iplp2 == 2)
  977.                                         {
  978.                                                 if (plot.iplp3 == 1)
  979.                                                         fprintf (
  980.                                                             plot_fp,
  981.                                                             "%12.4E %12.4E %12.4E\n",
  982.                                                             phi,
  983.                                                             ethm,
  984.                                                             etha);
  985.                                                 else if (plot.iplp3 == 2)
  986.                                                         fprintf (
  987.                                                             plot_fp,
  988.                                                             "%12.4E %12.4E %12.4E\n",
  989.                                                             phi,
  990.                                                             ephm,
  991.                                                             epha);
  992.                                         }
  993.                                 }
  994.  
  995.                                 if (plot.iplp4 == 0)
  996.                                         continue;
  997.  
  998.                                 if (plot.iplp2 == 1)
  999.                                 {
  1000.                                         switch (plot.iplp4)
  1001.                                         {
  1002.                                         case 1:
  1003.                                                 fprintf (
  1004.                                                     plot_fp, "%12.4E %12.4E\n", thet, tmp5);
  1005.                                                 break;
  1006.                                         case 2:
  1007.                                                 fprintf (
  1008.                                                     plot_fp, "%12.4E %12.4E\n", thet, tmp6);
  1009.                                                 break;
  1010.                                         case 3:
  1011.                                                 fprintf (
  1012.                                                     plot_fp, "%12.4E %12.4E\n", thet, gtot);
  1013.                                         }
  1014.                                 }
  1015.  
  1016.                                 if (plot.iplp2 == 2)
  1017.                                 {
  1018.                                         switch (plot.iplp4)
  1019.                                         {
  1020.                                         case 1:
  1021.                                                 fprintf (
  1022.                                                     plot_fp, "%12.4E %12.4E\n", phi, tmp5);
  1023.                                                 break;
  1024.                                         case 2:
  1025.                                                 fprintf (
  1026.                                                     plot_fp, "%12.4E %12.4E\n", phi, tmp6);
  1027.                                                 break;
  1028.                                         case 3:
  1029.                                                 fprintf (
  1030.                                                     plot_fp, "%12.4E %12.4E\n", phi, gtot);
  1031.                                         }
  1032.                                 }
  1033.  
  1034.                                 continue;
  1035.                         } /* if( gnd.ifar != 1) */
  1036.  
  1037.                         fprintf (
  1038.                             output_fp,
  1039.                             "\n"
  1040.                             " %9.2f %7.2f %9.2f  %11.4E %7.2f  %11.4E %7.2f  %11.4E %7.2f",
  1041.                             fpat.rfld,
  1042.                             phi,
  1043.                             thet,
  1044.                             ethm,
  1045.                             etha,
  1046.                             ephm,
  1047.                             epha,
  1048.                             erdm,
  1049.                             erda);
  1050.  
  1051.                 } /* for( kth = 1; kth <= fpat.nth; kth++ ) */
  1052.  
  1053.         } /* for( kph = 1; kph <= fpat.nph; kph++ ) */
  1054.  
  1055.         if (fpat.iavp != 0)
  1056.         {
  1057.                 tmp3 = fpat.thets * TA;
  1058.                 tmp4 = tmp3 + fpat.dth * TA * (double) (fpat.nth - 1);
  1059.                 tmp3 = fabsl (
  1060.                     fpat.dph * TA * (double) (fpat.nph - 1) * (cosl (tmp3) - cosl (tmp4)));
  1061.                 pint /= tmp3;
  1062.                 tmp3 /= PI;
  1063.  
  1064.                 fprintf (
  1065.                     output_fp,
  1066.                     "\n\n\n"
  1067.                     "   AVERAGE POWER GAIN= %11.5E       SOLID ANGLE"
  1068.                     " USED IN AVERAGING=( %6.4f)*PI STERADIANS.\n",
  1069.                     pint,
  1070.                     tmp3);
  1071.         }
  1072.  
  1073.         if (fpat.inor == 0)
  1074.                 return;
  1075.  
  1076.         if (fabsl (fpat.gnor) > 1.0e-20)
  1077.                 gmax = fpat.gnor;
  1078.         itmp1 = (fpat.inor - 1);
  1079.  
  1080.         fprintf (
  1081.             output_fp,
  1082.             "\n\n\n"
  1083.             "                             "
  1084.             " ---------- NORMALIZED GAIN ----------\n"
  1085.             "                                      %6s GAIN\n"
  1086.             "                                  "
  1087.             " NORMALIZATION FACTOR: %.2lf db\n\n"
  1088.             "    ---- ANGLES ----                ---- ANGLES ----"
  1089.             "                ---- ANGLES ----\n"
  1090.             "    THETA      PHI        GAIN      THETA      PHI  "
  1091.             "      GAIN      THETA      PHI       GAIN\n"
  1092.             "   DEGREES   DEGREES        DB     DEGREES   DEGREES "
  1093.             "       DB     DEGREES   DEGREES       DB",
  1094.             igntp[itmp1],
  1095.             gmax);
  1096.  
  1097.         itmp2 = fpat.nph * fpat.nth;
  1098.         itmp1 = (itmp2 + 2) / 3;
  1099.         itmp2 = itmp1 * 3 - itmp2;
  1100.         itmp3 = itmp1;
  1101.         itmp4 = 2 * itmp1;
  1102.  
  1103.         if (itmp2 == 2)
  1104.                 itmp4--;
  1105.  
  1106.         for (i = 0; i < itmp1; i++)
  1107.         {
  1108.                 itmp3++;
  1109.                 itmp4++;
  1110.                 j = i / fpat.nth;
  1111.                 tmp1 = fpat.thets + (double) (i - j * fpat.nth) * fpat.dth;
  1112.                 tmp2 = fpat.phis + (double) (j) *fpat.dph;
  1113.                 j = (itmp3 - 1) / fpat.nth;
  1114.                 tmp3 = fpat.thets + (double) (itmp3 - j * fpat.nth - 1) * fpat.dth;
  1115.                 tmp4 = fpat.phis + (double) (j) *fpat.dph;
  1116.                 j = (itmp4 - 1) / fpat.nth;
  1117.                 tmp5 = fpat.thets + (double) (itmp4 - j * fpat.nth - 1) * fpat.dth;
  1118.                 tmp6 = fpat.phis + (double) (j) *fpat.dph;
  1119.                 tstor1 = gain[i] - gmax;
  1120.  
  1121.                 if (((i + 1) == itmp1) && (itmp2 != 0))
  1122.                 {
  1123.                         if (itmp2 != 2)
  1124.                         {
  1125.                                 tstor2 = gain[itmp3 - 1] - gmax;
  1126.                                 fprintf (
  1127.                                     output_fp,
  1128.                                     "\n"
  1129.                                     " %9.2f %9.2f %9.2f   %9.2f %9.2f %9.2f   ",
  1130.                                     tmp1,
  1131.                                     tmp2,
  1132.                                     tstor1,
  1133.                                     tmp3,
  1134.                                     tmp4,
  1135.                                     tstor2);
  1136.                                 return;
  1137.                         }
  1138.  
  1139.                         fprintf (
  1140.                             output_fp,
  1141.                             "\n"
  1142.                             " %9.2f %9.2f %9.2f   ",
  1143.                             tmp1,
  1144.                             tmp2,
  1145.                             tstor1);
  1146.                         return;
  1147.  
  1148.                 } /* if( ((i+1) == itmp1) && (itmp2 != 0) ) */
  1149.  
  1150.                 tstor2 = gain[itmp3 - 1] - gmax;
  1151.                 pint = gain[itmp4 - 1] - gmax;
  1152.  
  1153.                 fprintf (
  1154.                     output_fp,
  1155.                     "\n"
  1156.                     " %9.2f %9.2f %9.2f   %9.2f %9.2f %9.2f   %9.2f %9.2f %9.2f",
  1157.                     tmp1,
  1158.                     tmp2,
  1159.                     tstor1,
  1160.                     tmp3,
  1161.                     tmp4,
  1162.                     tstor2,
  1163.                     tmp5,
  1164.                     tmp6,
  1165.                     pint);
  1166.  
  1167.         } /* for( i = 0; i < itmp1; i++ ) */
  1168.  
  1169.         free_ptr ((void *) &gain);
  1170.  
  1171.         return;
  1172. }
  1173.  
  1174. /*-----------------------------------------------------------------------*/
  1175.