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  /dataj/ */
  29. dataj_t dataj;
  30.  
  31. /* common  /gnd/ */
  32. extern gnd_t gnd;
  33.  
  34. /* common  /incom/ */
  35. extern incom_t incom;
  36.  
  37. /* common  /tmi/ */
  38. tmi_t tmi;
  39.  
  40. /*common  /tmh/ */
  41. static tmh_t tmh;
  42.  
  43. /* common  /gwav/ */
  44. extern gwav_t gwav;
  45.  
  46. /* common  /data/ */
  47. extern data_t data;
  48.  
  49. /* common  /crnt/ */
  50. extern crnt_t crnt;
  51.  
  52. /* common  /fpat/ */
  53. extern fpat_t fpat;
  54.  
  55. /* common  /plot/ */
  56. extern plot_t plot;
  57.  
  58. /* pointers to input/output files */
  59. extern FILE *input_fp, *output_fp, *plot_fp;
  60.  
  61. /*-------------------------------------------------------------------*/
  62.  
  63. /* compute near e fields of a segment with sine, cosine, and */
  64. /* constant currents.  ground effect included. */
  65. void efld (double xi, double yi, double zi, double ai, int ij)
  66. {
  67. #define txk egnd[0]
  68. #define tyk egnd[1]
  69. #define tzk egnd[2]
  70. #define txs egnd[3]
  71. #define tys egnd[4]
  72. #define tzs egnd[5]
  73. #define txc egnd[6]
  74. #define tyc egnd[7]
  75. #define tzc egnd[8]
  76.  
  77.         int ip;
  78.         double xij, yij, ijx, rfl, salpr, zij, zp, rhox;
  79.         double rhoy, rhoz, rh, r, rmag, cth, px, py;
  80.         double xymag, xspec, yspec, rhospc, dmin, shaf;
  81.         complex double epx, epy, refs, refps, zrsin, zratx, zscrn;
  82.         complex double tezs, ters, tezc, terc, tezk, terk, egnd[9];
  83.  
  84.         xij = xi - dataj.xj;
  85.         yij = yi - dataj.yj;
  86.         ijx = ij;
  87.         rfl = -1.;
  88.  
  89.         for (ip = 0; ip < gnd.ksymp; ip++)
  90.         {
  91.                 if (ip == 1)
  92.                         ijx = 1;
  93.                 rfl = -rfl;
  94.                 salpr = dataj.salpj * rfl;
  95.                 zij = zi - rfl * dataj.zj;
  96.                 zp = xij * dataj.cabj + yij * dataj.sabj + zij * salpr;
  97.                 rhox = xij - dataj.cabj * zp;
  98.                 rhoy = yij - dataj.sabj * zp;
  99.                 rhoz = zij - salpr * zp;
  100.  
  101.                 rh = sqrtl (rhox * rhox + rhoy * rhoy + rhoz * rhoz + ai * ai);
  102.                 if (rh <= 1.e-10)
  103.                 {
  104.                         rhox = 0.;
  105.                         rhoy = 0.;
  106.                         rhoz = 0.;
  107.                 }
  108.                 else
  109.                 {
  110.                         rhox = rhox / rh;
  111.                         rhoy = rhoy / rh;
  112.                         rhoz = rhoz / rh;
  113.                 }
  114.  
  115.                 /* lumped current element approx. for large separations */
  116.                 r = sqrtl (zp * zp + rh * rh);
  117.                 if (r >= dataj.rkh)
  118.                 {
  119.                         rmag = TP * r;
  120.                         cth = zp / r;
  121.                         px = rh / r;
  122.                         txk = cmplx (cosl (rmag), -sinl (rmag));
  123.                         py = TP * r * r;
  124.                         tyk = ETA * cth * txk * cmplx (1.0, -1.0 / rmag) / py;
  125.                         tzk = ETA * px * txk * cmplx (1.0, rmag - 1.0 / rmag) / (2. * py);
  126.                         tezk = tyk * cth - tzk * px;
  127.                         terk = tyk * px + tzk * cth;
  128.                         rmag = sinl (PI * dataj.s) / PI;
  129.                         tezc = tezk * rmag;
  130.                         terc = terk * rmag;
  131.                         tezk = tezk * dataj.s;
  132.                         terk = terk * dataj.s;
  133.                         txs = CPLX_00;
  134.                         tys = CPLX_00;
  135.                         tzs = CPLX_00;
  136.  
  137.                 } /* if( r >= dataj.rkh) */
  138.  
  139.                 if (r < dataj.rkh)
  140.                 {
  141.                         /* eksc for thin wire approx. or ekscx for extended t.w. approx. */
  142.                         if (dataj.iexk != 1)
  143.                                 eksc (
  144.                                     dataj.s,
  145.                                     zp,
  146.                                     rh,
  147.                                     TP,
  148.                                     ijx,
  149.                                     &tezs,
  150.                                     &ters,
  151.                                     &tezc,
  152.                                     &terc,
  153.                                     &tezk,
  154.                                     &terk);
  155.                         else
  156.                                 ekscx (
  157.                                     dataj.b,
  158.                                     dataj.s,
  159.                                     zp,
  160.                                     rh,
  161.                                     TP,
  162.                                     ijx,
  163.                                     dataj.ind1,
  164.                                     dataj.ind2,
  165.                                     &tezs,
  166.                                     &ters,
  167.                                     &tezc,
  168.                                     &terc,
  169.                                     &tezk,
  170.                                     &terk);
  171.  
  172.                         txs = tezs * dataj.cabj + ters * rhox;
  173.                         tys = tezs * dataj.sabj + ters * rhoy;
  174.                         tzs = tezs * salpr + ters * rhoz;
  175.  
  176.                 } /* if( r < dataj.rkh) */
  177.  
  178.                 txk = tezk * dataj.cabj + terk * rhox;
  179.                 tyk = tezk * dataj.sabj + terk * rhoy;
  180.                 tzk = tezk * salpr + terk * rhoz;
  181.                 txc = tezc * dataj.cabj + terc * rhox;
  182.                 tyc = tezc * dataj.sabj + terc * rhoy;
  183.                 tzc = tezc * salpr + terc * rhoz;
  184.  
  185.                 if (ip == 1)
  186.                 {
  187.                         if (gnd.iperf <= 0)
  188.                         {
  189.                                 zratx = gnd.zrati;
  190.                                 rmag = r;
  191.                                 xymag = sqrtl (xij * xij + yij * yij);
  192.  
  193.                                 /* set parameters for radial wire ground screen. */
  194.                                 if (gnd.nradl != 0)
  195.                                 {
  196.                                         xspec =
  197.                                             (xi * dataj.zj + zi * dataj.xj) / (zi + dataj.zj);
  198.                                         yspec =
  199.                                             (yi * dataj.zj + zi * dataj.yj) / (zi + dataj.zj);
  200.                                         rhospc = sqrtl (
  201.                                             xspec * xspec + yspec * yspec + gnd.t2 * gnd.t2);
  202.  
  203.                                         if (rhospc <= gnd.scrwl)
  204.                                         {
  205.                                                 zscrn =
  206.                                                     gnd.t1 * rhospc * logl (rhospc / gnd.t2);
  207.                                                 zratx = (zscrn * gnd.zrati) /
  208.                                                         (ETA * gnd.zrati + zscrn);
  209.                                         }
  210.                                 } /* if( gnd.nradl != 0) */
  211.  
  212.                                 /* calculation of reflection coefficients when ground is
  213.                                  * specified. */
  214.                                 if (xymag <= 1.0e-6)
  215.                                 {
  216.                                         px = 0.;
  217.                                         py = 0.;
  218.                                         cth = 1.;
  219.                                         zrsin = CPLX_10;
  220.                                 }
  221.                                 else
  222.                                 {
  223.                                         px = -yij / xymag;
  224.                                         py = xij / xymag;
  225.                                         cth = zij / rmag;
  226.                                         zrsin =
  227.                                             csqrtl (1.0 - zratx * zratx * (1.0 - cth * cth));
  228.  
  229.                                 } /* if( xymag <= 1.0e-6) */
  230.  
  231.                                 refs = (cth - zratx * zrsin) / (cth + zratx * zrsin);
  232.                                 refps = -(zratx * cth - zrsin) / (zratx * cth + zrsin);
  233.                                 refps = refps - refs;
  234.                                 epy = px * txk + py * tyk;
  235.                                 epx = px * epy;
  236.                                 epy = py * epy;
  237.                                 txk = refs * txk + refps * epx;
  238.                                 tyk = refs * tyk + refps * epy;
  239.                                 tzk = refs * tzk;
  240.                                 epy = px * txs + py * tys;
  241.                                 epx = px * epy;
  242.                                 epy = py * epy;
  243.                                 txs = refs * txs + refps * epx;
  244.                                 tys = refs * tys + refps * epy;
  245.                                 tzs = refs * tzs;
  246.                                 epy = px * txc + py * tyc;
  247.                                 epx = px * epy;
  248.                                 epy = py * epy;
  249.                                 txc = refs * txc + refps * epx;
  250.                                 tyc = refs * tyc + refps * epy;
  251.                                 tzc = refs * tzc;
  252.  
  253.                         } /* if( gnd.iperf <= 0) */
  254.  
  255.                         dataj.exk = dataj.exk - txk * gnd.frati;
  256.                         dataj.eyk = dataj.eyk - tyk * gnd.frati;
  257.                         dataj.ezk = dataj.ezk - tzk * gnd.frati;
  258.                         dataj.exs = dataj.exs - txs * gnd.frati;
  259.                         dataj.eys = dataj.eys - tys * gnd.frati;
  260.                         dataj.ezs = dataj.ezs - tzs * gnd.frati;
  261.                         dataj.exc = dataj.exc - txc * gnd.frati;
  262.                         dataj.eyc = dataj.eyc - tyc * gnd.frati;
  263.                         dataj.ezc = dataj.ezc - tzc * gnd.frati;
  264.                         continue;
  265.  
  266.                 } /* if( ip == 1) */
  267.  
  268.                 dataj.exk = txk;
  269.                 dataj.eyk = tyk;
  270.                 dataj.ezk = tzk;
  271.                 dataj.exs = txs;
  272.                 dataj.eys = tys;
  273.                 dataj.ezs = tzs;
  274.                 dataj.exc = txc;
  275.                 dataj.eyc = tyc;
  276.                 dataj.ezc = tzc;
  277.  
  278.         } /* for( ip = 0; ip < gnd.ksymp; ip++ ) */
  279.  
  280.         if (gnd.iperf != 2)
  281.                 return;
  282.  
  283.         /* field due to ground using sommerfeld/norton */
  284.         incom.sn = sqrtl (dataj.cabj * dataj.cabj + dataj.sabj * dataj.sabj);
  285.         if (incom.sn >= 1.0e-5)
  286.         {
  287.                 incom.xsn = dataj.cabj / incom.sn;
  288.                 incom.ysn = dataj.sabj / incom.sn;
  289.         }
  290.         else
  291.         {
  292.                 incom.sn = 0.;
  293.                 incom.xsn = 1.;
  294.                 incom.ysn = 0.;
  295.         }
  296.  
  297.         /* displace observation point for thin wire approximation */
  298.         zij = zi + dataj.zj;
  299.         salpr = -dataj.salpj;
  300.         rhox = dataj.sabj * zij - salpr * yij;
  301.         rhoy = salpr * xij - dataj.cabj * zij;
  302.         rhoz = dataj.cabj * yij - dataj.sabj * xij;
  303.         rh = rhox * rhox + rhoy * rhoy + rhoz * rhoz;
  304.  
  305.         if (rh <= 1.e-10)
  306.         {
  307.                 incom.xo = xi - ai * incom.ysn;
  308.                 incom.yo = yi + ai * incom.xsn;
  309.                 incom.zo = zi;
  310.         }
  311.         else
  312.         {
  313.                 rh = ai / sqrtl (rh);
  314.                 if (rhoz < 0.)
  315.                         rh = -rh;
  316.                 incom.xo = xi + rh * rhox;
  317.                 incom.yo = yi + rh * rhoy;
  318.                 incom.zo = zi + rh * rhoz;
  319.  
  320.         } /* if( rh <= 1.e-10) */
  321.  
  322.         r = xij * xij + yij * yij + zij * zij;
  323.         if (r <= .95)
  324.         {
  325.                 /* field from interpolation is integrated over segment */
  326.                 incom.isnor = 1;
  327.                 dmin = dataj.exk * conjl (dataj.exk) + dataj.eyk * conjl (dataj.eyk) +
  328.                        dataj.ezk * conjl (dataj.ezk);
  329.                 dmin = .01 * sqrtl (dmin);
  330.                 shaf = .5 * dataj.s;
  331.                 rom2 (-shaf, shaf, egnd, dmin);
  332.         }
  333.         else
  334.         {
  335.                 /* norton field equations and lumped current element approximation */
  336.                 incom.isnor = 2;
  337.                 sflds (0., egnd);
  338.         } /* if( r <= .95) */
  339.  
  340.         if (r > .95)
  341.         {
  342.                 zp = xij * dataj.cabj + yij * dataj.sabj + zij * salpr;
  343.                 rh = r - zp * zp;
  344.                 if (rh <= 1.e-10)
  345.                         dmin = 0.;
  346.                 else
  347.                         dmin = sqrtl (rh / (rh + ai * ai));
  348.  
  349.                 if (dmin <= .95)
  350.                 {
  351.                         px = 1. - dmin;
  352.                         terk = (txk * dataj.cabj + tyk * dataj.sabj + tzk * salpr) * px;
  353.                         txk = dmin * txk + terk * dataj.cabj;
  354.                         tyk = dmin * tyk + terk * dataj.sabj;
  355.                         tzk = dmin * tzk + terk * salpr;
  356.                         ters = (txs * dataj.cabj + tys * dataj.sabj + tzs * salpr) * px;
  357.                         txs = dmin * txs + ters * dataj.cabj;
  358.                         tys = dmin * tys + ters * dataj.sabj;
  359.                         tzs = dmin * tzs + ters * salpr;
  360.                         terc = (txc * dataj.cabj + tyc * dataj.sabj + tzc * salpr) * px;
  361.                         txc = dmin * txc + terc * dataj.cabj;
  362.                         tyc = dmin * tyc + terc * dataj.sabj;
  363.                         tzc = dmin * tzc + terc * salpr;
  364.  
  365.                 } /* if( dmin <= .95) */
  366.  
  367.         } /* if( r > .95) */
  368.  
  369.         dataj.exk = dataj.exk + txk;
  370.         dataj.eyk = dataj.eyk + tyk;
  371.         dataj.ezk = dataj.ezk + tzk;
  372.         dataj.exs = dataj.exs + txs;
  373.         dataj.eys = dataj.eys + tys;
  374.         dataj.ezs = dataj.ezs + tzs;
  375.         dataj.exc = dataj.exc + txc;
  376.         dataj.eyc = dataj.eyc + tyc;
  377.         dataj.ezc = dataj.ezc + tzc;
  378.  
  379.         return;
  380. }
  381.  
  382. /*-----------------------------------------------------------------------*/
  383.  
  384. /* compute e field of sine, cosine, and constant */
  385. /* current filaments by thin wire approximation. */
  386. void eksc (
  387.     double s,
  388.     double z,
  389.     double rh,
  390.     double xk,
  391.     int ij,
  392.     complex double *ezs,
  393.     complex double *ers,
  394.     complex double *ezc,
  395.     complex double *erc,
  396.     complex double *ezk,
  397.     complex double *erk)
  398. {
  399.         double rhk, sh, shk, ss, cs, z1a, z2a, cint, sint;
  400.         complex double gz1, gz2, gp1, gp2, gzp1, gzp2;
  401.  
  402.         tmi.ij = ij;
  403.         tmi.zpk = xk * z;
  404.         rhk = xk * rh;
  405.         tmi.rkb2 = rhk * rhk;
  406.         sh = .5 * s;
  407.         shk = xk * sh;
  408.         ss = sinl (shk);
  409.         cs = cosl (shk);
  410.         z2a = sh - z;
  411.         z1a = -(sh + z);
  412.         gx (z1a, rh, xk, &gz1, &gp1);
  413.         gx (z2a, rh, xk, &gz2, &gp2);
  414.         gzp1 = gp1 * z1a;
  415.         gzp2 = gp2 * z2a;
  416.         *ezs = CONST1 * ((gz2 - gz1) * cs * xk - (gzp2 + gzp1) * ss);
  417.         *ezc = -CONST1 * ((gz2 + gz1) * ss * xk + (gzp2 - gzp1) * cs);
  418.         *erk = CONST1 * (gp2 - gp1) * rh;
  419.         intx (-shk, shk, rhk, ij, &cint, &sint);
  420.         *ezk = -CONST1 * (gzp2 - gzp1 + xk * xk * cmplx (cint, -sint));
  421.         gzp1 = gzp1 * z1a;
  422.         gzp2 = gzp2 * z2a;
  423.  
  424.         if (rh >= 1.0e-10)
  425.         {
  426.                 *ers = -CONST1 *
  427.                        ((gzp2 + gzp1 + gz2 + gz1) * ss - (z2a * gz2 - z1a * gz1) * cs * xk) /
  428.                        rh;
  429.                 *erc = -CONST1 *
  430.                        ((gzp2 - gzp1 + gz2 - gz1) * cs + (z2a * gz2 + z1a * gz1) * ss * xk) /
  431.                        rh;
  432.                 return;
  433.         }
  434.  
  435.         *ers = CPLX_00;
  436.         *erc = CPLX_00;
  437.  
  438.         return;
  439. }
  440.  
  441. /*-----------------------------------------------------------------------*/
  442.  
  443. /* compute e field of sine, cosine, and constant current */
  444. /* filaments by extended thin wire approximation. */
  445. void ekscx (
  446.     double bx,
  447.     double s,
  448.     double z,
  449.     double rhx,
  450.     double xk,
  451.     int ij,
  452.     int inx1,
  453.     int inx2,
  454.     complex double *ezs,
  455.     complex double *ers,
  456.     complex double *ezc,
  457.     complex double *erc,
  458.     complex double *ezk,
  459.     complex double *erk)
  460. {
  461.         int ira;
  462.         double b, rh, sh, rhk, shk, ss, cs, z1a;
  463.         double z2a, a2, bk, bk2, cint, sint;
  464.         complex double gz1, gz2, gzp1, gzp2, gr1, gr2;
  465.         complex double grp1, grp2, grk1, grk2, gzz1, gzz2;
  466.  
  467.         if (rhx >= bx)
  468.         {
  469.                 rh = rhx;
  470.                 b = bx;
  471.                 ira = 0;
  472.         }
  473.         else
  474.         {
  475.                 rh = bx;
  476.                 b = rhx;
  477.                 ira = 1;
  478.         }
  479.  
  480.         sh = .5 * s;
  481.         tmi.ij = ij;
  482.         tmi.zpk = xk * z;
  483.         rhk = xk * rh;
  484.         tmi.rkb2 = rhk * rhk;
  485.         shk = xk * sh;
  486.         ss = sinl (shk);
  487.         cs = cosl (shk);
  488.         z2a = sh - z;
  489.         z1a = -(sh + z);
  490.         a2 = b * b;
  491.  
  492.         if (inx1 != 2)
  493.                 gxx (z1a, rh, b, a2, xk, ira, &gz1, &gzp1, &gr1, &grp1, &grk1, &gzz1);
  494.         else
  495.         {
  496.                 gx (z1a, rhx, xk, &gz1, &grk1);
  497.                 gzp1 = grk1 * z1a;
  498.                 gr1 = gz1 / rhx;
  499.                 grp1 = gzp1 / rhx;
  500.                 grk1 = grk1 * rhx;
  501.                 gzz1 = CPLX_00;
  502.         }
  503.  
  504.         if (inx2 != 2)
  505.                 gxx (z2a, rh, b, a2, xk, ira, &gz2, &gzp2, &gr2, &grp2, &grk2, &gzz2);
  506.         else
  507.         {
  508.                 gx (z2a, rhx, xk, &gz2, &grk2);
  509.                 gzp2 = grk2 * z2a;
  510.                 gr2 = gz2 / rhx;
  511.                 grp2 = gzp2 / rhx;
  512.                 grk2 = grk2 * rhx;
  513.                 gzz2 = CPLX_00;
  514.         }
  515.  
  516.         *ezs = CONST1 * ((gz2 - gz1) * cs * xk - (gzp2 + gzp1) * ss);
  517.         *ezc = -CONST1 * ((gz2 + gz1) * ss * xk + (gzp2 - gzp1) * cs);
  518.         *ers = -CONST1 * ((z2a * grp2 + z1a * grp1 + gr2 + gr1) * ss -
  519.                           (z2a * gr2 - z1a * gr1) * cs * xk);
  520.         *erc = -CONST1 * ((z2a * grp2 - z1a * grp1 + gr2 - gr1) * cs +
  521.                           (z2a * gr2 + z1a * gr1) * ss * xk);
  522.         *erk = CONST1 * (grk2 - grk1);
  523.         intx (-shk, shk, rhk, ij, &cint, &sint);
  524.         bk = b * xk;
  525.         bk2 = bk * bk * .25;
  526.         *ezk = -CONST1 * (gzp2 - gzp1 + xk * xk * (1. - bk2) * cmplx (cint, -sint) -
  527.                           bk2 * (gzz2 - gzz1));
  528.  
  529.         return;
  530. }
  531.  
  532. /*-----------------------------------------------------------------------*/
  533.  
  534. /* integrand for h field of a wire */
  535. void gh (double zk, double *hr, double *hi)
  536. {
  537.         double rs, r, ckr, skr, rr2, rr3;
  538.  
  539.         rs = zk - tmh.zpka;
  540.         rs = tmh.rhks + rs * rs;
  541.         r = sqrtl (rs);
  542.         ckr = cosl (r);
  543.         skr = sinl (r);
  544.         rr2 = 1. / rs;
  545.         rr3 = rr2 / r;
  546.         *hr = skr * rr2 + ckr * rr3;
  547.         *hi = ckr * rr2 - skr * rr3;
  548.  
  549.         return;
  550. }
  551.  
  552. /*-----------------------------------------------------------------------*/
  553.  
  554. /* gwave computes the electric field, including ground wave, of a */
  555. /* current element over a ground plane using formulas of k.a. norton */
  556. /* (proc. ire, sept., 1937, pp.1203,1236.) */
  557.  
  558. void gwave (
  559.     complex double *erv,
  560.     complex double *ezv,
  561.     complex double *erh,
  562.     complex double *ezh,
  563.     complex double *eph)
  564. {
  565.         double sppp, sppp2, cppp2, cppp, spp, spp2, cpp2, cpp;
  566.         complex double rk1, rk2, t1, t2, t3, t4, p1, rv;
  567.         complex double omr, w, f, q1, rh, v, g, xr1, xr2;
  568.         complex double x1, x2, x3, x4, x5, x6, x7;
  569.  
  570.         sppp = gwav.zmh / gwav.r1;
  571.         sppp2 = sppp * sppp;
  572.         cppp2 = 1. - sppp2;
  573.  
  574.         if (cppp2 < 1.0e-20)
  575.                 cppp2 = 1.0e-20;
  576.  
  577.         cppp = sqrtl (cppp2);
  578.         spp = gwav.zph / gwav.r2;
  579.         spp2 = spp * spp;
  580.         cpp2 = 1. - spp2;
  581.  
  582.         if (cpp2 < 1.0e-20)
  583.                 cpp2 = 1.0e-20;
  584.  
  585.         cpp = sqrtl (cpp2);
  586.         rk1 = -TPJ * gwav.r1;
  587.         rk2 = -TPJ * gwav.r2;
  588.         t1 = 1. - gwav.u2 * cpp2;
  589.         t2 = csqrtl (t1);
  590.         t3 = (1. - 1. / rk1) / rk1;
  591.         t4 = (1. - 1. / rk2) / rk2;
  592.         p1 = rk2 * gwav.u2 * t1 / (2. * cpp2);
  593.         rv = (spp - gwav.u * t2) / (spp + gwav.u * t2);
  594.         omr = 1. - rv;
  595.         w = 1. / omr;
  596.         w = (4.0 + 0.0fj) * p1 * w * w;
  597.         fbar (w, &f);
  598.         q1 = rk2 * t1 / (2. * gwav.u2 * cpp2);
  599.         rh = (t2 - gwav.u * spp) / (t2 + gwav.u * spp);
  600.         v = 1. / (1. + rh);
  601.         v = (4.0 + 0.0fj) * q1 * v * v;
  602.         fbar (v, &g);
  603.         xr1 = gwav.xx1 / gwav.r1;
  604.         xr2 = gwav.xx2 / gwav.r2;
  605.         x1 = cppp2 * xr1;
  606.         x2 = rv * cpp2 * xr2;
  607.         x3 = omr * cpp2 * f * xr2;
  608.         x4 = gwav.u * t2 * spp * 2. * xr2 / rk2;
  609.         x5 = xr1 * t3 * (1. - 3. * sppp2);
  610.         x6 = xr2 * t4 * (1. - 3. * spp2);
  611.         *ezv = (x1 + x2 + x3 - x4 - x5 - x6) * (-CONST4);
  612.         x1 = sppp * cppp * xr1;
  613.         x2 = rv * spp * cpp * xr2;
  614.         x3 = cpp * omr * gwav.u * t2 * f * xr2;
  615.         x4 = spp * cpp * omr * xr2 / rk2;
  616.         x5 = 3. * sppp * cppp * t3 * xr1;
  617.         x6 = cpp * gwav.u * t2 * omr * xr2 / rk2 * .5;
  618.         x7 = 3. * spp * cpp * t4 * xr2;
  619.         *erv = -(x1 + x2 - x3 + x4 - x5 + x6 - x7) * (-CONST4);
  620.         *ezh = -(x1 - x2 + x3 - x4 - x5 - x6 + x7) * (-CONST4);
  621.         x1 = sppp2 * xr1;
  622.         x2 = rv * spp2 * xr2;
  623.         x4 = gwav.u2 * t1 * omr * f * xr2;
  624.         x5 = t3 * (1. - 3. * cppp2) * xr1;
  625.         x6 = t4 * (1. - 3. * cpp2) * (1. - gwav.u2 * (1. + rv) - gwav.u2 * omr * f) * xr2;
  626.         x7 = gwav.u2 * cpp2 * omr * (1. - 1. / rk2) *
  627.              (f * (gwav.u2 * t1 - spp2 - 1. / rk2) + 1. / rk2) * xr2;
  628.         *erh = (x1 - x2 - x4 - x5 + x6 + x7) * (-CONST4);
  629.         x1 = xr1;
  630.         x2 = rh * xr2;
  631.         x3 = (rh + 1.) * g * xr2;
  632.         x4 = t3 * xr1;
  633.         x5 = t4 * (1. - gwav.u2 * (1. + rv) - gwav.u2 * omr * f) * xr2;
  634.         x6 =
  635.             .5 * gwav.u2 * omr * (f * (gwav.u2 * t1 - spp2 - 1. / rk2) + 1. / rk2) * xr2 / rk2;
  636.         *eph = -(x1 - x2 + x3 - x4 + x5 + x6) * (-CONST4);
  637.  
  638.         return;
  639. }
  640.  
  641. /*-----------------------------------------------------------------------*/
  642.  
  643. /* segment end contributions for thin wire approx. */
  644. void gx (double zz, double rh, double xk, complex double *gz, complex double *gzp)
  645. {
  646.         double r, r2, rkz;
  647.  
  648.         r2 = zz * zz + rh * rh;
  649.         r = sqrtl (r2);
  650.         rkz = xk * r;
  651.         *gz = cmplx (cosl (rkz), -sinl (rkz)) / r;
  652.         *gzp = -cmplx (1.0, rkz) * *gz / r2;
  653.  
  654.         return;
  655. }
  656.  
  657. /*-----------------------------------------------------------------------*/
  658.  
  659. /* segment end contributions for ext. thin wire approx. */
  660. void gxx (
  661.     double zz,
  662.     double rh,
  663.     double a,
  664.     double a2,
  665.     double xk,
  666.     int ira,
  667.     complex double *g1,
  668.     complex double *g1p,
  669.     complex double *g2,
  670.     complex double *g2p,
  671.     complex double *g3,
  672.     complex double *gzp)
  673. {
  674.         double r, r2, r4, rk, rk2, rh2, t1, t2;
  675.         complex double gz, c1, c2, c3;
  676.  
  677.         r2 = zz * zz + rh * rh;
  678.         r = sqrtl (r2);
  679.         r4 = r2 * r2;
  680.         rk = xk * r;
  681.         rk2 = rk * rk;
  682.         rh2 = rh * rh;
  683.         t1 = .25 * a2 * rh2 / r4;
  684.         t2 = .5 * a2 / r2;
  685.         c1 = cmplx (1.0, rk);
  686.         c2 = 3. * c1 - rk2;
  687.         c3 = cmplx (6.0, rk) * rk2 - 15. * c1;
  688.         gz = cmplx (cosl (rk), -sinl (rk)) / r;
  689.         *g2 = gz * (1. + t1 * c2);
  690.         *g1 = *g2 - t2 * c1 * gz;
  691.         gz = gz / r2;
  692.         *g2p = gz * (t1 * c3 - c1);
  693.         *gzp = t2 * c2 * gz;
  694.         *g3 = *g2p + *gzp;
  695.         *g1p = *g3 * zz;
  696.  
  697.         if (ira != 1)
  698.         {
  699.                 *g3 = (*g3 + *gzp) * rh;
  700.                 *gzp = -zz * c1 * gz;
  701.  
  702.                 if (rh <= 1.0e-10)
  703.                 {
  704.                         *g2 = 0.;
  705.                         *g2p = 0.;
  706.                         return;
  707.                 }
  708.  
  709.                 *g2 = *g2 / rh;
  710.                 *g2p = *g2p * zz / rh;
  711.                 return;
  712.  
  713.         } /* if( ira != 1) */
  714.  
  715.         t2 = .5 * a;
  716.         *g2 = -t2 * c1 * gz;
  717.         *g2p = t2 * gz * c2 / r2;
  718.         *g3 = rh2 * *g2p - a * gz * c1;
  719.         *g2p = *g2p * zz;
  720.         *gzp = -zz * c1 * gz;
  721.  
  722.         return;
  723. }
  724.  
  725. /*-----------------------------------------------------------------------*/
  726.  
  727. /* hfk computes the h field of a uniform current */
  728. /* filament by numerical integration */
  729. void hfk (double el1, double el2, double rhk, double zpkx, double *sgr, double *sgi)
  730. {
  731.         int nx = 1, nma = 65536, nts = 4;
  732.         int ns, nt;
  733.         int flag = TRUE;
  734.         double rx = 1.0e-4;
  735.         double z, ze, s, ep, zend, dz = 0., zp, dzot = 0., t00r, g1r, g5r = 0, t00i;
  736.         double g1i, g5i = 0., t01r, g3r = 0, t01i, g3i = 0, t10r, t10i, te1i, te1r, t02r;
  737.         double g2r, g4r, t02i, g2i, g4i, t11r, t11i, t20r, t20i, te2i, te2r;
  738.  
  739.         tmh.zpka = zpkx;
  740.         tmh.rhks = rhk * rhk;
  741.         z = el1;
  742.         ze = el2;
  743.         s = ze - z;
  744.         ep = s / (10. * nma);
  745.         zend = ze - ep;
  746.         *sgr = 0.0;
  747.         *sgi = 0.0;
  748.         ns = nx;
  749.         nt = 0;
  750.         gh (z, &g1r, &g1i);
  751.  
  752.         while (TRUE)
  753.         {
  754.                 if (flag)
  755.                 {
  756.                         dz = s / ns;
  757.                         zp = z + dz;
  758.  
  759.                         if (zp > ze)
  760.                         {
  761.                                 dz = ze - z;
  762.                                 if (fabsl (dz) <= ep)
  763.                                 {
  764.                                         *sgr = *sgr * rhk * .5;
  765.                                         *sgi = *sgi * rhk * .5;
  766.                                         return;
  767.                                 }
  768.                         }
  769.  
  770.                         dzot = dz * .5;
  771.                         zp = z + dzot;
  772.                         gh (zp, &g3r, &g3i);
  773.                         zp = z + dz;
  774.                         gh (zp, &g5r, &g5i);
  775.  
  776.                 } /* if( flag ) */
  777.  
  778.                 t00r = (g1r + g5r) * dzot;
  779.                 t00i = (g1i + g5i) * dzot;
  780.                 t01r = (t00r + dz * g3r) * 0.5;
  781.                 t01i = (t00i + dz * g3i) * 0.5;
  782.                 t10r = (4.0 * t01r - t00r) / 3.0;
  783.                 t10i = (4.0 * t01i - t00i) / 3.0;
  784.  
  785.                 test (t01r, t10r, &te1r, t01i, t10i, &te1i, 0.);
  786.                 if ((te1i <= rx) && (te1r <= rx))
  787.                 {
  788.                         *sgr = *sgr + t10r;
  789.                         *sgi = *sgi + t10i;
  790.                         nt += 2;
  791.  
  792.                         z += dz;
  793.                         if (z >= zend)
  794.                         {
  795.                                 *sgr = *sgr * rhk * .5;
  796.                                 *sgi = *sgi * rhk * .5;
  797.                                 return;
  798.                         }
  799.  
  800.                         g1r = g5r;
  801.                         g1i = g5i;
  802.                         if (nt >= nts)
  803.                                 if (ns > nx)
  804.                                 {
  805.                                         ns = ns / 2;
  806.                                         nt = 1;
  807.                                 }
  808.                         flag = TRUE;
  809.                         continue;
  810.  
  811.                 } /* if( (te1i <= rx) && (te1r <= rx) ) */
  812.  
  813.                 zp = z + dz * 0.25;
  814.                 gh (zp, &g2r, &g2i);
  815.                 zp = z + dz * 0.75;
  816.                 gh (zp, &g4r, &g4i);
  817.                 t02r = (t01r + dzot * (g2r + g4r)) * 0.5;
  818.                 t02i = (t01i + dzot * (g2i + g4i)) * 0.5;
  819.                 t11r = (4.0 * t02r - t01r) / 3.0;
  820.                 t11i = (4.0 * t02i - t01i) / 3.0;
  821.                 t20r = (16.0 * t11r - t10r) / 15.0;
  822.                 t20i = (16.0 * t11i - t10i) / 15.0;
  823.  
  824.                 test (t11r, t20r, &te2r, t11i, t20i, &te2i, 0.);
  825.                 if ((te2i > rx) || (te2r > rx))
  826.                 {
  827.                         nt = 0;
  828.                         if (ns >= nma)
  829.                                 fprintf (output_fp, "\n  STEP SIZE LIMITED AT Z= %10.5f", z);
  830.                         else
  831.                         {
  832.                                 ns = ns * 2;
  833.                                 dz = s / ns;
  834.                                 dzot = dz * 0.5;
  835.                                 g5r = g3r;
  836.                                 g5i = g3i;
  837.                                 g3r = g2r;
  838.                                 g3i = g2i;
  839.  
  840.                                 flag = FALSE;
  841.                                 continue;
  842.                         }
  843.  
  844.                 } /* if( (te2i > rx) || (te2r > rx) ) */
  845.  
  846.                 *sgr = *sgr + t20r;
  847.                 *sgi = *sgi + t20i;
  848.                 nt++;
  849.  
  850.                 z += dz;
  851.                 if (z >= zend)
  852.                 {
  853.                         *sgr = *sgr * rhk * .5;
  854.                         *sgi = *sgi * rhk * .5;
  855.                         return;
  856.                 }
  857.  
  858.                 g1r = g5r;
  859.                 g1i = g5i;
  860.                 if (nt >= nts)
  861.                         if (ns > nx)
  862.                         {
  863.                                 ns = ns / 2;
  864.                                 nt = 1;
  865.                         }
  866.                 flag = TRUE;
  867.  
  868.         } /* while( TRUE ) */
  869. }
  870.  
  871. /*-----------------------------------------------------------------------*/
  872.  
  873. /* hintg computes the h field of a patch current */
  874. void hintg (double xi, double yi, double zi)
  875. {
  876.         int ip;
  877.         double rx, ry, rfl, xymag, pxx, pyy, cth;
  878.         double rz, rsq, r, rk, cr, sr, t1zr, t2zr;
  879.         complex double gam, f1x, f1y, f1z, f2x, f2y, f2z, rrv, rrh;
  880.  
  881.         rx = xi - dataj.xj;
  882.         ry = yi - dataj.yj;
  883.         rfl = -1.;
  884.         dataj.exk = CPLX_00;
  885.         dataj.eyk = CPLX_00;
  886.         dataj.ezk = CPLX_00;
  887.         dataj.exs = CPLX_00;
  888.         dataj.eys = CPLX_00;
  889.         dataj.ezs = CPLX_00;
  890.  
  891.         for (ip = 1; ip <= gnd.ksymp; ip++)
  892.         {
  893.                 rfl = -rfl;
  894.                 rz = zi - dataj.zj * rfl;
  895.                 rsq = rx * rx + ry * ry + rz * rz;
  896.  
  897.                 if (rsq < 1.0e-20)
  898.                         continue;
  899.  
  900.                 r = sqrtl (rsq);
  901.                 rk = TP * r;
  902.                 cr = cosl (rk);
  903.                 sr = sinl (rk);
  904.                 gam = -(cmplx (cr, -sr) + rk * cmplx (sr, cr)) / (FPI * rsq * r) * dataj.s;
  905.                 dataj.exc = gam * rx;
  906.                 dataj.eyc = gam * ry;
  907.                 dataj.ezc = gam * rz;
  908.                 t1zr = dataj.t1zj * rfl;
  909.                 t2zr = dataj.t2zj * rfl;
  910.                 f1x = dataj.eyc * t1zr - dataj.ezc * dataj.t1yj;
  911.                 f1y = dataj.ezc * dataj.t1xj - dataj.exc * t1zr;
  912.                 f1z = dataj.exc * dataj.t1yj - dataj.eyc * dataj.t1xj;
  913.                 f2x = dataj.eyc * t2zr - dataj.ezc * dataj.t2yj;
  914.                 f2y = dataj.ezc * dataj.t2xj - dataj.exc * t2zr;
  915.                 f2z = dataj.exc * dataj.t2yj - dataj.eyc * dataj.t2xj;
  916.  
  917.                 if (ip != 1)
  918.                 {
  919.                         if (gnd.iperf == 1)
  920.                         {
  921.                                 f1x = -f1x;
  922.                                 f1y = -f1y;
  923.                                 f1z = -f1z;
  924.                                 f2x = -f2x;
  925.                                 f2y = -f2y;
  926.                                 f2z = -f2z;
  927.                         }
  928.                         else
  929.                         {
  930.                                 xymag = sqrtl (rx * rx + ry * ry);
  931.                                 if (xymag <= 1.0e-6)
  932.                                 {
  933.                                         pxx = 0.;
  934.                                         pyy = 0.;
  935.                                         cth = 1.;
  936.                                         rrv = CPLX_10;
  937.                                 }
  938.                                 else
  939.                                 {
  940.                                         pxx = -ry / xymag;
  941.                                         pyy = rx / xymag;
  942.                                         cth = rz / r;
  943.                                         rrv = csqrtl (
  944.                                             1. - gnd.zrati * gnd.zrati * (1. - cth * cth));
  945.  
  946.                                 } /* if( xymag <= 1.0e-6) */
  947.  
  948.                                 rrh = gnd.zrati * cth;
  949.                                 rrh = (rrh - rrv) / (rrh + rrv);
  950.                                 rrv = gnd.zrati * rrv;
  951.                                 rrv = -(cth - rrv) / (cth + rrv);
  952.                                 gam = (f1x * pxx + f1y * pyy) * (rrv - rrh);
  953.                                 f1x = f1x * rrh + gam * pxx;
  954.                                 f1y = f1y * rrh + gam * pyy;
  955.                                 f1z = f1z * rrh;
  956.                                 gam = (f2x * pxx + f2y * pyy) * (rrv - rrh);
  957.                                 f2x = f2x * rrh + gam * pxx;
  958.                                 f2y = f2y * rrh + gam * pyy;
  959.                                 f2z = f2z * rrh;
  960.  
  961.                         } /* if( gnd.iperf == 1) */
  962.  
  963.                 } /* if( ip != 1) */
  964.  
  965.                 dataj.exk += f1x;
  966.                 dataj.eyk += f1y;
  967.                 dataj.ezk += f1z;
  968.                 dataj.exs += f2x;
  969.                 dataj.eys += f2y;
  970.                 dataj.ezs += f2z;
  971.  
  972.         } /* for( ip = 1; ip <= gnd.ksymp; ip++ ) */
  973.  
  974.         return;
  975. }
  976.  
  977. /*-----------------------------------------------------------------------*/
  978.  
  979. /* hsfld computes the h field for constant, sine, and */
  980. /* cosine current on a segment including ground effects. */
  981. void hsfld (double xi, double yi, double zi, double ai)
  982. {
  983.         int ip;
  984.         double xij, yij, rfl, salpr, zij, zp, rhox, rhoy, rhoz, rh, phx;
  985.         double phy, phz, rmag, xymag, xspec, yspec, rhospc, px, py, cth;
  986.         complex double hpk, hps, hpc, qx, qy, qz, rrv, rrh, zratx;
  987.  
  988.         xij = xi - dataj.xj;
  989.         yij = yi - dataj.yj;
  990.         rfl = -1.;
  991.  
  992.         for (ip = 0; ip < gnd.ksymp; ip++)
  993.         {
  994.                 rfl = -rfl;
  995.                 salpr = dataj.salpj * rfl;
  996.                 zij = zi - rfl * dataj.zj;
  997.                 zp = xij * dataj.cabj + yij * dataj.sabj + zij * salpr;
  998.                 rhox = xij - dataj.cabj * zp;
  999.                 rhoy = yij - dataj.sabj * zp;
  1000.                 rhoz = zij - salpr * zp;
  1001.                 rh = sqrtl (rhox * rhox + rhoy * rhoy + rhoz * rhoz + ai * ai);
  1002.  
  1003.                 if (rh <= 1.0e-10)
  1004.                 {
  1005.                         dataj.exk = 0.;
  1006.                         dataj.eyk = 0.;
  1007.                         dataj.ezk = 0.;
  1008.                         dataj.exs = 0.;
  1009.                         dataj.eys = 0.;
  1010.                         dataj.ezs = 0.;
  1011.                         dataj.exc = 0.;
  1012.                         dataj.eyc = 0.;
  1013.                         dataj.ezc = 0.;
  1014.                         continue;
  1015.                 }
  1016.  
  1017.                 rhox = rhox / rh;
  1018.                 rhoy = rhoy / rh;
  1019.                 rhoz = rhoz / rh;
  1020.                 phx = dataj.sabj * rhoz - salpr * rhoy;
  1021.                 phy = salpr * rhox - dataj.cabj * rhoz;
  1022.                 phz = dataj.cabj * rhoy - dataj.sabj * rhox;
  1023.  
  1024.                 hsflx (dataj.s, rh, zp, &hpk, &hps, &hpc);
  1025.  
  1026.                 if (ip == 1)
  1027.                 {
  1028.                         if (gnd.iperf != 1)
  1029.                         {
  1030.                                 zratx = gnd.zrati;
  1031.                                 rmag = sqrtl (zp * zp + rh * rh);
  1032.                                 xymag = sqrtl (xij * xij + yij * yij);
  1033.  
  1034.                                 /* set parameters for radial wire ground screen. */
  1035.                                 if (gnd.nradl != 0)
  1036.                                 {
  1037.                                         xspec =
  1038.                                             (xi * dataj.zj + zi * dataj.xj) / (zi + dataj.zj);
  1039.                                         yspec =
  1040.                                             (yi * dataj.zj + zi * dataj.yj) / (zi + dataj.zj);
  1041.                                         rhospc = sqrtl (
  1042.                                             xspec * xspec + yspec * yspec + gnd.t2 * gnd.t2);
  1043.  
  1044.                                         if (rhospc <= gnd.scrwl)
  1045.                                         {
  1046.                                                 rrv = gnd.t1 * rhospc * logl (rhospc / gnd.t2);
  1047.                                                 zratx = (rrv * gnd.zrati) /
  1048.                                                         (ETA * gnd.zrati + rrv);
  1049.                                         }
  1050.                                 }
  1051.  
  1052.                                 /* calculation of reflection coefficients when ground is
  1053.                                  * specified. */
  1054.                                 if (xymag <= 1.0e-6)
  1055.                                 {
  1056.                                         px = 0.;
  1057.                                         py = 0.;
  1058.                                         cth = 1.;
  1059.                                         rrv = CPLX_10;
  1060.                                 }
  1061.                                 else
  1062.                                 {
  1063.                                         px = -yij / xymag;
  1064.                                         py = xij / xymag;
  1065.                                         cth = zij / rmag;
  1066.                                         rrv = csqrtl (1. - zratx * zratx * (1. - cth * cth));
  1067.                                 }
  1068.  
  1069.                                 rrh = zratx * cth;
  1070.                                 rrh = -(rrh - rrv) / (rrh + rrv);
  1071.                                 rrv = zratx * rrv;
  1072.                                 rrv = (cth - rrv) / (cth + rrv);
  1073.                                 qy = (phx * px + phy * py) * (rrv - rrh);
  1074.                                 qx = qy * px + phx * rrh;
  1075.                                 qy = qy * py + phy * rrh;
  1076.                                 qz = phz * rrh;
  1077.                                 dataj.exk = dataj.exk - hpk * qx;
  1078.                                 dataj.eyk = dataj.eyk - hpk * qy;
  1079.                                 dataj.ezk = dataj.ezk - hpk * qz;
  1080.                                 dataj.exs = dataj.exs - hps * qx;
  1081.                                 dataj.eys = dataj.eys - hps * qy;
  1082.                                 dataj.ezs = dataj.ezs - hps * qz;
  1083.                                 dataj.exc = dataj.exc - hpc * qx;
  1084.                                 dataj.eyc = dataj.eyc - hpc * qy;
  1085.                                 dataj.ezc = dataj.ezc - hpc * qz;
  1086.                                 continue;
  1087.  
  1088.                         } /* if( gnd.iperf != 1 ) */
  1089.  
  1090.                         dataj.exk = dataj.exk - hpk * phx;
  1091.                         dataj.eyk = dataj.eyk - hpk * phy;
  1092.                         dataj.ezk = dataj.ezk - hpk * phz;
  1093.                         dataj.exs = dataj.exs - hps * phx;
  1094.                         dataj.eys = dataj.eys - hps * phy;
  1095.                         dataj.ezs = dataj.ezs - hps * phz;
  1096.                         dataj.exc = dataj.exc - hpc * phx;
  1097.                         dataj.eyc = dataj.eyc - hpc * phy;
  1098.                         dataj.ezc = dataj.ezc - hpc * phz;
  1099.                         continue;
  1100.  
  1101.                 } /* if( ip == 1 ) */
  1102.  
  1103.                 dataj.exk = hpk * phx;
  1104.                 dataj.eyk = hpk * phy;
  1105.                 dataj.ezk = hpk * phz;
  1106.                 dataj.exs = hps * phx;
  1107.                 dataj.eys = hps * phy;
  1108.                 dataj.ezs = hps * phz;
  1109.                 dataj.exc = hpc * phx;
  1110.                 dataj.eyc = hpc * phy;
  1111.                 dataj.ezc = hpc * phz;
  1112.  
  1113.         } /* for( ip = 0; ip < gnd.ksymp; ip++ ) */
  1114.  
  1115.         return;
  1116. }
  1117.  
  1118. /*-----------------------------------------------------------------------*/
  1119.  
  1120. /* calculates h field of sine cosine, and constant current of segment */
  1121. void hsflx (
  1122.     double s,
  1123.     double rh,
  1124.     double zpx,
  1125.     complex double *hpk,
  1126.     complex double *hps,
  1127.     complex double *hpc)
  1128. {
  1129.         double r1, r2, zp, z2a, hss, dh, z1;
  1130.         double rhz, dk, cdk, sdk, hkr, hki, rh2;
  1131.         complex double fjk, ekr1, ekr2, t1, t2, cons;
  1132.  
  1133.         fjk = -TPJ;
  1134.         if (rh >= 1.0e-10)
  1135.         {
  1136.                 if (zpx >= 0.)
  1137.                 {
  1138.                         zp = zpx;
  1139.                         hss = 1.;
  1140.                 }
  1141.                 else
  1142.                 {
  1143.                         zp = -zpx;
  1144.                         hss = -1.;
  1145.                 }
  1146.  
  1147.                 dh = .5 * s;
  1148.                 z1 = zp + dh;
  1149.                 z2a = zp - dh;
  1150.                 if (z2a >= 1.0e-7)
  1151.                         rhz = rh / z2a;
  1152.                 else
  1153.                         rhz = 1.;
  1154.  
  1155.                 dk = TP * dh;
  1156.                 cdk = cosl (dk);
  1157.                 sdk = sinl (dk);
  1158.                 hfk (-dk, dk, rh * TP, zp * TP, &hkr, &hki);
  1159.                 *hpk = cmplx (hkr, hki);
  1160.  
  1161.                 if (rhz >= 1.0e-3)
  1162.                 {
  1163.                         rh2 = rh * rh;
  1164.                         r1 = sqrtl (rh2 + z1 * z1);
  1165.                         r2 = sqrtl (rh2 + z2a * z2a);
  1166.                         ekr1 = cexp (fjk * r1);
  1167.                         ekr2 = cexp (fjk * r2);
  1168.                         t1 = z1 * ekr1 / r1;
  1169.                         t2 = z2a * ekr2 / r2;
  1170.                         *hps = (cdk * (ekr2 - ekr1) - CPLX_01 * sdk * (t2 + t1)) * hss;
  1171.                         *hpc = -sdk * (ekr2 + ekr1) - CPLX_01 * cdk * (t2 - t1);
  1172.                         cons = -CPLX_01 / (2. * TP * rh);
  1173.                         *hps = cons * *hps;
  1174.                         *hpc = cons * *hpc;
  1175.                         return;
  1176.  
  1177.                 } /* if( rhz >= 1.0e-3) */
  1178.  
  1179.                 ekr1 = cmplx (cdk, sdk) / (z2a * z2a);
  1180.                 ekr2 = cmplx (cdk, -sdk) / (z1 * z1);
  1181.                 t1 = TP * (1. / z1 - 1. / z2a);
  1182.                 t2 = cexp (fjk * zp) * rh / PI8;
  1183.                 *hps = t2 * (t1 + (ekr1 + ekr2) * sdk) * hss;
  1184.                 *hpc = t2 * (-CPLX_01 * t1 + (ekr1 - ekr2) * cdk);
  1185.                 return;
  1186.  
  1187.         } /* if( rh >= 1.0e-10) */
  1188.  
  1189.         *hps = CPLX_00;
  1190.         *hpc = CPLX_00;
  1191.         *hpk = CPLX_00;
  1192.  
  1193.         return;
  1194. }
  1195.  
  1196. /*-----------------------------------------------------------------------*/
  1197.  
  1198. /* nefld computes the near field at specified points in space after */
  1199. /* the structure currents have been computed. */
  1200. void nefld (
  1201.     double xob,
  1202.     double yob,
  1203.     double zob,
  1204.     complex double *ex,
  1205.     complex double *ey,
  1206.     complex double *ez)
  1207. {
  1208.         int i, ix, ipr, iprx, jc, ipa;
  1209.         double zp, xi, ax;
  1210.         complex double acx, bcx, ccx;
  1211.  
  1212.         *ex = CPLX_00;
  1213.         *ey = CPLX_00;
  1214.         *ez = CPLX_00;
  1215.         ax = 0.;
  1216.  
  1217.         if (data.n != 0)
  1218.         {
  1219.                 for (i = 0; i < data.n; i++)
  1220.                 {
  1221.                         dataj.xj = xob - data.x[i];
  1222.                         dataj.yj = yob - data.y[i];
  1223.                         dataj.zj = zob - data.z[i];
  1224.                         zp = data.cab[i] * dataj.xj + data.sab[i] * dataj.yj +
  1225.                              data.salp[i] * dataj.zj;
  1226.  
  1227.                         if (fabsl (zp) > 0.5001 * data.si[i])
  1228.                                 continue;
  1229.  
  1230.                         zp = dataj.xj * dataj.xj + dataj.yj * dataj.yj + dataj.zj * dataj.zj -
  1231.                              zp * zp;
  1232.                         dataj.xj = data.bi[i];
  1233.  
  1234.                         if (zp > 0.9 * dataj.xj * dataj.xj)
  1235.                                 continue;
  1236.  
  1237.                         ax = dataj.xj;
  1238.                         break;
  1239.  
  1240.                 } /* for( i = 0; i < n; i++ ) */
  1241.  
  1242.                 for (i = 0; i < data.n; i++)
  1243.                 {
  1244.                         ix = i + 1;
  1245.                         dataj.s = data.si[i];
  1246.                         dataj.b = data.bi[i];
  1247.                         dataj.xj = data.x[i];
  1248.                         dataj.yj = data.y[i];
  1249.                         dataj.zj = data.z[i];
  1250.                         dataj.cabj = data.cab[i];
  1251.                         dataj.sabj = data.sab[i];
  1252.                         dataj.salpj = data.salp[i];
  1253.  
  1254.                         if (dataj.iexk != 0)
  1255.                         {
  1256.                                 ipr = data.icon1[i];
  1257.  
  1258.                                 if (ipr > PCHCON)
  1259.                                         dataj.ind1 = 2;
  1260.                                 else if (ipr < 0)
  1261.                                 {
  1262.                                         ipr = -ipr;
  1263.                                         iprx = ipr - 1;
  1264.  
  1265.                                         if (-data.icon1[iprx] != ix)
  1266.                                                 dataj.ind1 = 2;
  1267.                                         else
  1268.                                         {
  1269.                                                 xi = fabsl (
  1270.                                                     dataj.cabj * data.cab[iprx] +
  1271.                                                     dataj.sabj * data.sab[iprx] +
  1272.                                                     dataj.salpj * data.salp[iprx]);
  1273.                                                 if ((xi < 0.999999) ||
  1274.                                                     (fabsl (data.bi[iprx] / dataj.b - 1.) >
  1275.                                                      1.0e-6))
  1276.                                                         dataj.ind1 = 2;
  1277.                                                 else
  1278.                                                         dataj.ind1 = 0;
  1279.                                         }
  1280.                                 } /* if( ipr < 0 ) */
  1281.                                 else if (ipr == 0)
  1282.                                         dataj.ind1 = 1;
  1283.                                 else
  1284.                                 {
  1285.                                         iprx = ipr - 1;
  1286.  
  1287.                                         if (ipr != ix)
  1288.                                         {
  1289.                                                 if (data.icon2[iprx] != ix)
  1290.                                                         dataj.ind1 = 2;
  1291.                                                 else
  1292.                                                 {
  1293.                                                         xi = fabsl (
  1294.                                                             dataj.cabj * data.cab[iprx] +
  1295.                                                             dataj.sabj * data.sab[iprx] +
  1296.                                                             dataj.salpj * data.salp[iprx]);
  1297.                                                         if ((xi < 0.999999) ||
  1298.                                                             (fabsl (
  1299.                                                                  data.bi[iprx] / dataj.b -
  1300.                                                                  1.) > 1.0e-6))
  1301.                                                                 dataj.ind1 = 2;
  1302.                                                         else
  1303.                                                                 dataj.ind1 = 0;
  1304.                                                 }
  1305.                                         } /* if( ipr != ix ) */
  1306.                                         else
  1307.                                         {
  1308.                                                 if (dataj.cabj * dataj.cabj +
  1309.                                                         dataj.sabj * dataj.sabj >
  1310.                                                     1.0e-8)
  1311.                                                         dataj.ind1 = 2;
  1312.                                                 else
  1313.                                                         dataj.ind1 = 0;
  1314.                                         }
  1315.                                 } /* else */
  1316.  
  1317.                                 ipr = data.icon2[i];
  1318.  
  1319.                                 if (ipr > PCHCON)
  1320.                                         dataj.ind2 = 2;
  1321.                                 else if (ipr < 0)
  1322.                                 {
  1323.                                         ipr = -ipr;
  1324.                                         iprx = ipr - 1;
  1325.  
  1326.                                         if (-data.icon2[iprx] != ix)
  1327.                                                 dataj.ind1 = 2;
  1328.                                         else
  1329.                                         {
  1330.                                                 xi = fabsl (
  1331.                                                     dataj.cabj * data.cab[iprx] +
  1332.                                                     dataj.sabj * data.sab[iprx] +
  1333.                                                     dataj.salpj * data.salp[iprx]);
  1334.                                                 if ((xi < 0.999999) ||
  1335.                                                     (fabsl (data.bi[iprx] / dataj.b - 1.) >
  1336.                                                      1.0e-6))
  1337.                                                         dataj.ind1 = 2;
  1338.                                                 else
  1339.                                                         dataj.ind1 = 0;
  1340.                                         }
  1341.                                 } /* if( ipr < 0 ) */
  1342.                                 else if (ipr == 0)
  1343.                                         dataj.ind2 = 1;
  1344.                                 else
  1345.                                 {
  1346.                                         iprx = ipr - 1;
  1347.  
  1348.                                         if (ipr != ix)
  1349.                                         {
  1350.                                                 if (data.icon1[iprx] != ix)
  1351.                                                         dataj.ind2 = 2;
  1352.                                                 else
  1353.                                                 {
  1354.                                                         xi = fabsl (
  1355.                                                             dataj.cabj * data.cab[iprx] +
  1356.                                                             dataj.sabj * data.sab[iprx] +
  1357.                                                             dataj.salpj * data.salp[iprx]);
  1358.                                                         if ((xi < 0.999999) ||
  1359.                                                             (fabsl (
  1360.                                                                  data.bi[iprx] / dataj.b -
  1361.                                                                  1.) > 1.0e-6))
  1362.                                                                 dataj.ind2 = 2;
  1363.                                                         else
  1364.                                                                 dataj.ind2 = 0;
  1365.                                                 }
  1366.                                         } /* if( ipr != (i+1) ) */
  1367.                                         else
  1368.                                         {
  1369.                                                 if (dataj.cabj * dataj.cabj +
  1370.                                                         dataj.sabj * dataj.sabj >
  1371.                                                     1.0e-8)
  1372.                                                         dataj.ind1 = 2;
  1373.                                                 else
  1374.                                                         dataj.ind1 = 0;
  1375.                                         }
  1376.  
  1377.                                 } /* else */
  1378.  
  1379.                         } /* if( dataj.iexk != 0) */
  1380.  
  1381.                         efld (xob, yob, zob, ax, 1);
  1382.                         acx = cmplx (crnt.air[i], crnt.aii[i]);
  1383.                         bcx = cmplx (crnt.bir[i], crnt.bii[i]);
  1384.                         ccx = cmplx (crnt.cir[i], crnt.cii[i]);
  1385.                         *ex += dataj.exk * acx + dataj.exs * bcx + dataj.exc * ccx;
  1386.                         *ey += dataj.eyk * acx + dataj.eys * bcx + dataj.eyc * ccx;
  1387.                         *ez += dataj.ezk * acx + dataj.ezs * bcx + dataj.ezc * ccx;
  1388.  
  1389.                 } /* for( i = 0; i < n; i++ ) */
  1390.  
  1391.                 if (data.m == 0)
  1392.                         return;
  1393.  
  1394.         } /* if( n != 0) */
  1395.  
  1396.         jc = data.n - 1;
  1397.         for (i = 0; i < data.m; i++)
  1398.         {
  1399.                 dataj.s = data.pbi[i];
  1400.                 dataj.xj = data.px[i];
  1401.                 dataj.yj = data.py[i];
  1402.                 dataj.zj = data.pz[i];
  1403.                 dataj.t1xj = data.t1x[i];
  1404.                 dataj.t1yj = data.t1y[i];
  1405.                 dataj.t1zj = data.t1z[i];
  1406.                 dataj.t2xj = data.t2x[i];
  1407.                 dataj.t2yj = data.t2y[i];
  1408.                 dataj.t2zj = data.t2z[i];
  1409.                 jc += 3;
  1410.                 acx = dataj.t1xj * crnt.cur[jc - 2] + dataj.t1yj * crnt.cur[jc - 1] +
  1411.                       dataj.t1zj * crnt.cur[jc];
  1412.                 bcx = dataj.t2xj * crnt.cur[jc - 2] + dataj.t2yj * crnt.cur[jc - 1] +
  1413.                       dataj.t2zj * crnt.cur[jc];
  1414.  
  1415.                 for (ipa = 0; ipa < gnd.ksymp; ipa++)
  1416.                 {
  1417.                         dataj.ipgnd = ipa + 1;
  1418.                         unere (xob, yob, zob);
  1419.                         *ex = *ex + acx * dataj.exk + bcx * dataj.exs;
  1420.                         *ey = *ey + acx * dataj.eyk + bcx * dataj.eys;
  1421.                         *ez = *ez + acx * dataj.ezk + bcx * dataj.ezs;
  1422.                 }
  1423.  
  1424.         } /* for( i = 0; i < m; i++ ) */
  1425.  
  1426.         return;
  1427. }
  1428.  
  1429. /*-----------------------------------------------------------------------*/
  1430.  
  1431. /* compute near e or h fields over a range of points */
  1432. void nfpat (void)
  1433. {
  1434.         int i, j, kk;
  1435.         double znrt, cth = 0., sth = 0., ynrt, cph = 0., sph = 0., xnrt, xob, yob;
  1436.         double zob, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, xxx;
  1437.         complex double ex, ey, ez;
  1438.  
  1439.         if (fpat.nfeh != 1)
  1440.         {
  1441.                 fprintf (
  1442.                     output_fp,
  1443.                     "\n\n\n"
  1444.                     "                             "
  1445.                     "-------- NEAR ELECTRIC FIELDS --------\n"
  1446.                     "     ------- LOCATION -------     ------- EX ------    ------- EY ------ "
  1447.                     "   ------- EZ ------\n"
  1448.                     "      X         Y         Z       MAGNITUDE   PHASE    MAGNITUDE   PHASE "
  1449.                     "   MAGNITUDE   PHASE\n"
  1450.                     "    METERS    METERS    METERS     VOLTS/M  DEGREES    VOLTS/M   DEGREES "
  1451.                     "    VOLTS/M  DEGREES");
  1452.         }
  1453.         else
  1454.         {
  1455.                 fprintf (
  1456.                     output_fp,
  1457.                     "\n\n\n"
  1458.                     "                                   "
  1459.                     "-------- NEAR MAGNETIC FIELDS ---------\n\n"
  1460.                     "     ------- LOCATION -------     ------- HX ------    ------- HY ------ "
  1461.                     "   ------- HZ ------\n"
  1462.                     "      X         Y         Z       MAGNITUDE   PHASE    MAGNITUDE   PHASE "
  1463.                     "   MAGNITUDE   PHASE\n"
  1464.                     "    METERS    METERS    METERS      AMPS/M  DEGREES      AMPS/M  DEGREES "
  1465.                     "     AMPS/M  DEGREES");
  1466.         }
  1467.  
  1468.         znrt = fpat.znr - fpat.dznr;
  1469.         for (i = 0; i < fpat.nrz; i++)
  1470.         {
  1471.                 znrt += fpat.dznr;
  1472.                 if (fpat.near != 0)
  1473.                 {
  1474.                         cth = cosl (TA * znrt);
  1475.                         sth = sinl (TA * znrt);
  1476.                 }
  1477.  
  1478.                 ynrt = fpat.ynr - fpat.dynr;
  1479.                 for (j = 0; j < fpat.nry; j++)
  1480.                 {
  1481.                         ynrt += fpat.dynr;
  1482.                         if (fpat.near != 0)
  1483.                         {
  1484.                                 cph = cosl (TA * ynrt);
  1485.                                 sph = sinl (TA * ynrt);
  1486.                         }
  1487.  
  1488.                         xnrt = fpat.xnr - fpat.dxnr;
  1489.                         for (kk = 0; kk < fpat.nrx; kk++)
  1490.                         {
  1491.                                 xnrt += fpat.dxnr;
  1492.                                 if (fpat.near != 0)
  1493.                                 {
  1494.                                         xob = xnrt * sth * cph;
  1495.                                         yob = xnrt * sth * sph;
  1496.                                         zob = xnrt * cth;
  1497.                                 }
  1498.                                 else
  1499.                                 {
  1500.                                         xob = xnrt;
  1501.                                         yob = ynrt;
  1502.                                         zob = znrt;
  1503.                                 }
  1504.  
  1505.                                 tmp1 = xob / data.wlam;
  1506.                                 tmp2 = yob / data.wlam;
  1507.                                 tmp3 = zob / data.wlam;
  1508.  
  1509.                                 if (fpat.nfeh != 1)
  1510.                                         nefld (tmp1, tmp2, tmp3, &ex, &ey, &ez);
  1511.                                 else
  1512.                                         nhfld (tmp1, tmp2, tmp3, &ex, &ey, &ez);
  1513.  
  1514.                                 tmp1 = cabsl (ex);
  1515.                                 tmp2 = cang (ex);
  1516.                                 tmp3 = cabsl (ey);
  1517.                                 tmp4 = cang (ey);
  1518.                                 tmp5 = cabsl (ez);
  1519.                                 tmp6 = cang (ez);
  1520.  
  1521.                                 fprintf (
  1522.                                     output_fp,
  1523.                                     "\n"
  1524.                                     " %9.4f %9.4f %9.4f  %11.4E %7.2f  %11.4E %7.2f  %11.4E "
  1525.                                     "%7.2f",
  1526.                                     xob,
  1527.                                     yob,
  1528.                                     zob,
  1529.                                     tmp1,
  1530.                                     tmp2,
  1531.                                     tmp3,
  1532.                                     tmp4,
  1533.                                     tmp5,
  1534.                                     tmp6);
  1535.  
  1536.                                 if (plot.iplp1 != 2)
  1537.                                         continue;
  1538.  
  1539.                                 if (plot.iplp4 < 0)
  1540.                                         xxx = xob;
  1541.                                 else if (plot.iplp4 == 0)
  1542.                                         xxx = yob;
  1543.                                 else
  1544.                                         xxx = zob;
  1545.  
  1546.                                 if (plot.iplp2 == 2)
  1547.                                 {
  1548.                                         switch (plot.iplp3)
  1549.                                         {
  1550.                                         case 1:
  1551.                                                 fprintf (
  1552.                                                     plot_fp,
  1553.                                                     "%12.4E %12.4E %12.4E\n",
  1554.                                                     xxx,
  1555.                                                     tmp1,
  1556.                                                     tmp2);
  1557.                                                 break;
  1558.                                         case 2:
  1559.                                                 fprintf (
  1560.                                                     plot_fp,
  1561.                                                     "%12.4E %12.4E %12.4E\n",
  1562.                                                     xxx,
  1563.                                                     tmp3,
  1564.                                                     tmp4);
  1565.                                                 break;
  1566.                                         case 3:
  1567.                                                 fprintf (
  1568.                                                     plot_fp,
  1569.                                                     "%12.4E %12.4E %12.4E\n",
  1570.                                                     xxx,
  1571.                                                     tmp5,
  1572.                                                     tmp6);
  1573.                                                 break;
  1574.                                         case 4:
  1575.                                                 fprintf (
  1576.                                                     plot_fp,
  1577.                                                     "%12.4E %12.4E %12.4E %12.4E %12.4E "
  1578.                                                     "%12.4E %12.4E\n",
  1579.                                                     xxx,
  1580.                                                     tmp1,
  1581.                                                     tmp2,
  1582.                                                     tmp3,
  1583.                                                     tmp4,
  1584.                                                     tmp5,
  1585.                                                     tmp6);
  1586.                                         }
  1587.                                         continue;
  1588.                                 }
  1589.  
  1590.                                 if (plot.iplp2 != 1)
  1591.                                         continue;
  1592.  
  1593.                                 switch (plot.iplp3)
  1594.                                 {
  1595.                                 case 1:
  1596.                                         fprintf (
  1597.                                             plot_fp,
  1598.                                             "%12.4E %12.4E %12.4E\n",
  1599.                                             xxx,
  1600.                                             print_creall (ex),
  1601.                                             print_cimagl (ex));
  1602.                                         break;
  1603.                                 case 2:
  1604.                                         fprintf (
  1605.                                             plot_fp,
  1606.                                             "%12.4E %12.4E %12.4E\n",
  1607.                                             xxx,
  1608.                                             print_creall (ey),
  1609.                                             print_cimagl (ey));
  1610.                                         break;
  1611.                                 case 3:
  1612.                                         fprintf (
  1613.                                             plot_fp,
  1614.                                             "%12.4E %12.4E %12.4E\n",
  1615.                                             xxx,
  1616.                                             print_creall (ez),
  1617.                                             print_cimagl (ez));
  1618.                                         break;
  1619.                                 case 4:
  1620.                                         fprintf (
  1621.                                             plot_fp,
  1622.                                             "%12.4E %12.4E %12.4E %12.4E %12.4E %12.4E "
  1623.                                             "%12.4E\n",
  1624.                                             xxx,
  1625.                                             print_creall (ex),
  1626.                                             print_cimagl (ex),
  1627.                                             print_creall (ey),
  1628.                                             print_cimagl (ey),
  1629.                                             print_creall (ez),
  1630.                                             print_cimagl (ez));
  1631.                                 }
  1632.                         } /* for( kk = 0; kk < fpat.nrx; kk++ ) */
  1633.  
  1634.                 } /* for( j = 0; j < fpat.nry; j++ ) */
  1635.  
  1636.         } /* for( i = 0; i < fpat.nrz; i++ ) */
  1637.  
  1638.         return;
  1639. }
  1640.  
  1641. /*-----------------------------------------------------------------------*/
  1642.  
  1643. /* nhfld computes the near field at specified points in space after */
  1644. /* the structure currents have been computed. */
  1645.  
  1646. void nhfld (
  1647.     double xob,
  1648.     double yob,
  1649.     double zob,
  1650.     complex double *hx,
  1651.     complex double *hy,
  1652.     complex double *hz)
  1653. {
  1654.         int i, jc;
  1655.         double ax, zp;
  1656.         complex double acx, bcx, ccx;
  1657.  
  1658.         *hx = CPLX_00;
  1659.         *hy = CPLX_00;
  1660.         *hz = CPLX_00;
  1661.         ax = 0.;
  1662.  
  1663.         if (data.n != 0)
  1664.         {
  1665.                 for (i = 0; i < data.n; i++)
  1666.                 {
  1667.                         dataj.xj = xob - data.x[i];
  1668.                         dataj.yj = yob - data.y[i];
  1669.                         dataj.zj = zob - data.z[i];
  1670.                         zp = data.cab[i] * dataj.xj + data.sab[i] * dataj.yj +
  1671.                              data.salp[i] * dataj.zj;
  1672.  
  1673.                         if (fabsl (zp) > 0.5001 * data.si[i])
  1674.                                 continue;
  1675.  
  1676.                         zp = dataj.xj * dataj.xj + dataj.yj * dataj.yj + dataj.zj * dataj.zj -
  1677.                              zp * zp;
  1678.                         dataj.xj = data.bi[i];
  1679.  
  1680.                         if (zp > 0.9 * dataj.xj * dataj.xj)
  1681.                                 continue;
  1682.  
  1683.                         ax = dataj.xj;
  1684.                         break;
  1685.                 }
  1686.  
  1687.                 for (i = 0; i < data.n; i++)
  1688.                 {
  1689.                         dataj.s = data.si[i];
  1690.                         dataj.b = data.bi[i];
  1691.                         dataj.xj = data.x[i];
  1692.                         dataj.yj = data.y[i];
  1693.                         dataj.zj = data.z[i];
  1694.                         dataj.cabj = data.cab[i];
  1695.                         dataj.sabj = data.sab[i];
  1696.                         dataj.salpj = data.salp[i];
  1697.                         hsfld (xob, yob, zob, ax);
  1698.                         acx = cmplx (crnt.air[i], crnt.aii[i]);
  1699.                         bcx = cmplx (crnt.bir[i], crnt.bii[i]);
  1700.                         ccx = cmplx (crnt.cir[i], crnt.cii[i]);
  1701.                         *hx += dataj.exk * acx + dataj.exs * bcx + dataj.exc * ccx;
  1702.                         *hy += dataj.eyk * acx + dataj.eys * bcx + dataj.eyc * ccx;
  1703.                         *hz += dataj.ezk * acx + dataj.ezs * bcx + dataj.ezc * ccx;
  1704.                 }
  1705.  
  1706.                 if (data.m == 0)
  1707.                         return;
  1708.  
  1709.         } /* if( data.n != 0) */
  1710.  
  1711.         jc = data.n - 1;
  1712.         for (i = 0; i < data.m; i++)
  1713.         {
  1714.                 dataj.s = data.pbi[i];
  1715.                 dataj.xj = data.px[i];
  1716.                 dataj.yj = data.py[i];
  1717.                 dataj.zj = data.pz[i];
  1718.                 dataj.t1xj = data.t1x[i];
  1719.                 dataj.t1yj = data.t1y[i];
  1720.                 dataj.t1zj = data.t1z[i];
  1721.                 dataj.t2xj = data.t2x[i];
  1722.                 dataj.t2yj = data.t2y[i];
  1723.                 dataj.t2zj = data.t2z[i];
  1724.                 hintg (xob, yob, zob);
  1725.                 jc += 3;
  1726.                 acx = dataj.t1xj * crnt.cur[jc - 2] + dataj.t1yj * crnt.cur[jc - 1] +
  1727.                       dataj.t1zj * crnt.cur[jc];
  1728.                 bcx = dataj.t2xj * crnt.cur[jc - 2] + dataj.t2yj * crnt.cur[jc - 1] +
  1729.                       dataj.t2zj * crnt.cur[jc];
  1730.                 *hx = *hx + acx * dataj.exk + bcx * dataj.exs;
  1731.                 *hy = *hy + acx * dataj.eyk + bcx * dataj.eys;
  1732.                 *hz = *hz + acx * dataj.ezk + bcx * dataj.ezs;
  1733.         }
  1734.  
  1735.         return;
  1736. }
  1737.  
  1738. /*-----------------------------------------------------------------------*/
  1739.  
  1740. /* integrate over patches at wire connection point */
  1741. void pcint (
  1742.     double xi, double yi, double zi, double cabi, double sabi, double salpi, complex double *e)
  1743. {
  1744.         int nint, i1, i2;
  1745.         double d, ds, da, gcon, fcon, xxj, xyj, xzj, xs, s1;
  1746.         double xss, yss, zss, s2x, s2, g1, g2, g3, g4, f2, f1;
  1747.         complex double e1, e2, e3, e4, e5, e6, e7, e8, e9;
  1748.  
  1749.         nint = 10;
  1750.         d = sqrtl (dataj.s) * .5;
  1751.         ds = 4. * d / (double) nint;
  1752.         da = ds * ds;
  1753.         gcon = 1. / dataj.s;
  1754.         fcon = 1. / (2. * TP * d);
  1755.         xxj = dataj.xj;
  1756.         xyj = dataj.yj;
  1757.         xzj = dataj.zj;
  1758.         xs = dataj.s;
  1759.         dataj.s = da;
  1760.         s1 = d + ds * .5;
  1761.         xss = dataj.xj + s1 * (dataj.t1xj + dataj.t2xj);
  1762.         yss = dataj.yj + s1 * (dataj.t1yj + dataj.t2yj);
  1763.         zss = dataj.zj + s1 * (dataj.t1zj + dataj.t2zj);
  1764.         s1 = s1 + d;
  1765.         s2x = s1;
  1766.         e1 = CPLX_00;
  1767.         e2 = CPLX_00;
  1768.         e3 = CPLX_00;
  1769.         e4 = CPLX_00;
  1770.         e5 = CPLX_00;
  1771.         e6 = CPLX_00;
  1772.         e7 = CPLX_00;
  1773.         e8 = CPLX_00;
  1774.         e9 = CPLX_00;
  1775.  
  1776.         for (i1 = 0; i1 < nint; i1++)
  1777.         {
  1778.                 s1 = s1 - ds;
  1779.                 s2 = s2x;
  1780.                 xss = xss - ds * dataj.t1xj;
  1781.                 yss = yss - ds * dataj.t1yj;
  1782.                 zss = zss - ds * dataj.t1zj;
  1783.                 dataj.xj = xss;
  1784.                 dataj.yj = yss;
  1785.                 dataj.zj = zss;
  1786.  
  1787.                 for (i2 = 0; i2 < nint; i2++)
  1788.                 {
  1789.                         s2 = s2 - ds;
  1790.                         dataj.xj = dataj.xj - ds * dataj.t2xj;
  1791.                         dataj.yj = dataj.yj - ds * dataj.t2yj;
  1792.                         dataj.zj = dataj.zj - ds * dataj.t2zj;
  1793.                         unere (xi, yi, zi);
  1794.                         dataj.exk = dataj.exk * cabi + dataj.eyk * sabi + dataj.ezk * salpi;
  1795.                         dataj.exs = dataj.exs * cabi + dataj.eys * sabi + dataj.ezs * salpi;
  1796.                         g1 = (d + s1) * (d + s2) * gcon;
  1797.                         g2 = (d - s1) * (d + s2) * gcon;
  1798.                         g3 = (d - s1) * (d - s2) * gcon;
  1799.                         g4 = (d + s1) * (d - s2) * gcon;
  1800.                         f2 = (s1 * s1 + s2 * s2) * TP;
  1801.                         f1 = s1 / f2 - (g1 - g2 - g3 + g4) * fcon;
  1802.                         f2 = s2 / f2 - (g1 + g2 - g3 - g4) * fcon;
  1803.                         e1 = e1 + dataj.exk * g1;
  1804.                         e2 = e2 + dataj.exk * g2;
  1805.                         e3 = e3 + dataj.exk * g3;
  1806.                         e4 = e4 + dataj.exk * g4;
  1807.                         e5 = e5 + dataj.exs * g1;
  1808.                         e6 = e6 + dataj.exs * g2;
  1809.                         e7 = e7 + dataj.exs * g3;
  1810.                         e8 = e8 + dataj.exs * g4;
  1811.                         e9 = e9 + dataj.exk * f1 + dataj.exs * f2;
  1812.  
  1813.                 } /* for( i2 = 0; i2 < nint; i2++ ) */
  1814.  
  1815.         } /* for( i1 = 0; i1 < nint; i1++ ) */
  1816.  
  1817.         e[0] = e1;
  1818.         e[1] = e2;
  1819.         e[2] = e3;
  1820.         e[3] = e4;
  1821.         e[4] = e5;
  1822.         e[5] = e6;
  1823.         e[6] = e7;
  1824.         e[7] = e8;
  1825.         e[8] = e9;
  1826.         dataj.xj = xxj;
  1827.         dataj.yj = xyj;
  1828.         dataj.zj = xzj;
  1829.         dataj.s = xs;
  1830.  
  1831.         return;
  1832. }
  1833.  
  1834. /*-----------------------------------------------------------------------*/
  1835.  
  1836. /* calculates the electric field due to unit current */
  1837. /* in the t1 and t2 directions on a patch */
  1838. void unere (double xob, double yob, double zob)
  1839. {
  1840.         double zr, t1zr, t2zr, rx, ry, rz, r, tt1;
  1841.         double tt2, rt, xymag, px, py, cth, r2;
  1842.         complex double er, q1, q2, rrv, rrh, edp;
  1843.  
  1844.         zr = dataj.zj;
  1845.         t1zr = dataj.t1zj;
  1846.         t2zr = dataj.t2zj;
  1847.  
  1848.         if (dataj.ipgnd == 2)
  1849.         {
  1850.                 zr = -zr;
  1851.                 t1zr = -t1zr;
  1852.                 t2zr = -t2zr;
  1853.         }
  1854.  
  1855.         rx = xob - dataj.xj;
  1856.         ry = yob - dataj.yj;
  1857.         rz = zob - zr;
  1858.         r2 = rx * rx + ry * ry + rz * rz;
  1859.  
  1860.         if (r2 <= 1.0e-20)
  1861.         {
  1862.                 dataj.exk = CPLX_00;
  1863.                 dataj.eyk = CPLX_00;
  1864.                 dataj.ezk = CPLX_00;
  1865.                 dataj.exs = CPLX_00;
  1866.                 dataj.eys = CPLX_00;
  1867.                 dataj.ezs = CPLX_00;
  1868.                 return;
  1869.         }
  1870.  
  1871.         r = sqrtl (r2);
  1872.         tt1 = -TP * r;
  1873.         tt2 = tt1 * tt1;
  1874.         rt = r2 * r;
  1875.         er = cmplx (sinl (tt1), -cosl (tt1)) * (CONST2 * dataj.s);
  1876.         q1 = cmplx (tt2 - 1., tt1) * er / rt;
  1877.         q2 = cmplx (3. - tt2, -3. * tt1) * er / (rt * r2);
  1878.         er = q2 * (dataj.t1xj * rx + dataj.t1yj * ry + t1zr * rz);
  1879.         dataj.exk = q1 * dataj.t1xj + er * rx;
  1880.         dataj.eyk = q1 * dataj.t1yj + er * ry;
  1881.         dataj.ezk = q1 * t1zr + er * rz;
  1882.         er = q2 * (dataj.t2xj * rx + dataj.t2yj * ry + t2zr * rz);
  1883.         dataj.exs = q1 * dataj.t2xj + er * rx;
  1884.         dataj.eys = q1 * dataj.t2yj + er * ry;
  1885.         dataj.ezs = q1 * t2zr + er * rz;
  1886.  
  1887.         if (dataj.ipgnd == 1)
  1888.                 return;
  1889.  
  1890.         if (gnd.iperf == 1)
  1891.         {
  1892.                 dataj.exk = -dataj.exk;
  1893.                 dataj.eyk = -dataj.eyk;
  1894.                 dataj.ezk = -dataj.ezk;
  1895.                 dataj.exs = -dataj.exs;
  1896.                 dataj.eys = -dataj.eys;
  1897.                 dataj.ezs = -dataj.ezs;
  1898.                 return;
  1899.         }
  1900.  
  1901.         xymag = sqrtl (rx * rx + ry * ry);
  1902.         if (xymag <= 1.0e-6)
  1903.         {
  1904.                 px = 0.;
  1905.                 py = 0.;
  1906.                 cth = 1.;
  1907.                 rrv = CPLX_10;
  1908.         }
  1909.         else
  1910.         {
  1911.                 px = -ry / xymag;
  1912.                 py = rx / xymag;
  1913.                 cth = rz / sqrtl (xymag * xymag + rz * rz);
  1914.                 rrv = csqrtl (1. - gnd.zrati * gnd.zrati * (1. - cth * cth));
  1915.         }
  1916.  
  1917.         rrh = gnd.zrati * cth;
  1918.         rrh = (rrh - rrv) / (rrh + rrv);
  1919.         rrv = gnd.zrati * rrv;
  1920.         rrv = -(cth - rrv) / (cth + rrv);
  1921.         edp = (dataj.exk * px + dataj.eyk * py) * (rrh - rrv);
  1922.         dataj.exk = dataj.exk * rrv + edp * px;
  1923.         dataj.eyk = dataj.eyk * rrv + edp * py;
  1924.         dataj.ezk = dataj.ezk * rrv;
  1925.         edp = (dataj.exs * px + dataj.eys * py) * (rrh - rrv);
  1926.         dataj.exs = dataj.exs * rrv + edp * px;
  1927.         dataj.eys = dataj.eys * rrv + edp * py;
  1928.         dataj.ezs = dataj.ezs * rrv;
  1929.  
  1930.         return;
  1931. }
  1932.  
  1933. /*-----------------------------------------------------------------------*/
  1934.