Subversion Repositories Nec2c

Rev

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

  1. /* last change:  pgm   8 nov 2000    1:04 pm
  2.  program somnec(input,output,tape21)
  3.  
  4.  program to generate nec interpolation grids for fields due to
  5.  ground.  field components are computed by numerical evaluation
  6.  of modified sommerfeld integrals.
  7.  
  8.  somnec2d is a double precision version of somnec for use with
  9.  nec2d.  an alternate version (somnec2sd) is also provided in which
  10.  computation is in single precision but the output file is written
  11.  in double precision for use with nec2d.  somnec2sd runs about twic
  12.  as fast as the full double precision somnec2d.  the difference
  13.  between nec2d results using a for021 file from this code rather
  14.  than from somnec2sd was insignficant in the cases tested.
  15.  
  16.  changes made by j bergervoet, 31-5-95:
  17.  parameter 0. --> 0.d0 in calling of routine test
  18.  status of output files set to 'unknown' */
  19.  
  20. #include "nec2c.h"
  21.  
  22. /* common /evlcom/ */
  23. static int jh;
  24. static double ck2, ck2sq, tkmag, tsmag, ck1r, zph, rho;
  25. static complex double ct1, ct2, ct3, ck1, ck1sq, cksm;
  26.  
  27. /* common /cntour/ */
  28. static complex double a, b;
  29.  
  30. /*common  /ggrid/ */
  31. ggrid_t ggrid;
  32.  
  33. /*-----------------------------------------------------------------------*/
  34.  
  35. /* This is the "main" of somnec */
  36. void somnec (double epr, double sig, double fmhz)
  37. {
  38.         int k, nth, ith, irs, ir, nr;
  39.         double tim, wlam, tst, dr, dth, r, rk, thet, tfac1, tfac2;
  40.         complex double erv, ezv, erh, eph, cl1, cl2, con;
  41.  
  42.         if (sig >= 0.)
  43.         {
  44.                 wlam = CVEL / fmhz;
  45.                 ggrid.epscf = cmplx (epr, -sig * wlam * 59.96);
  46.         }
  47.         else
  48.                 ggrid.epscf = cmplx (epr, sig);
  49.  
  50.         secnds (&tst);
  51.         ck2 = TP;
  52.         ck2sq = ck2 * ck2;
  53.  
  54.         /* sommerfeld integral evaluation uses exp(-jwt), nec uses exp(+jwt), */
  55.         /* hence need conjg(ggrid.epscf).  conjugate of fields occurs in subroutine */
  56.         /* evlua. */
  57.  
  58.         ck1sq = ck2sq * conj (ggrid.epscf);
  59.         ck1 = csqrtl (ck1sq);
  60.         ck1r = creal (ck1);
  61.         tkmag = 100. * cabs (ck1);
  62.         tsmag = 100. * ck1 * conj (ck1);
  63.         cksm = ck2sq / (ck1sq + ck2sq);
  64.         ct1 = .5 * (ck1sq - ck2sq);
  65.         erv = ck1sq * ck1sq;
  66.         ezv = ck2sq * ck2sq;
  67.         ct2 = .125 * (erv - ezv);
  68.         erv *= ck1sq;
  69.         ezv *= ck2sq;
  70.         ct3 = .0625 * (erv - ezv);
  71.  
  72.         /* loop over 3 grid regions */
  73.         for (k = 0; k < 3; k++)
  74.         {
  75.                 nr = ggrid.nxa[k];
  76.                 nth = ggrid.nya[k];
  77.                 dr = ggrid.dxa[k];
  78.                 dth = ggrid.dya[k];
  79.                 r = ggrid.xsa[k] - dr;
  80.                 irs = 1;
  81.                 if (k == 0)
  82.                 {
  83.                         r = ggrid.xsa[k];
  84.                         irs = 2;
  85.                 }
  86.  
  87.                 /*  loop over r.  (r=sqrtl(rho**2 + (z+h)**2)) */
  88.                 for (ir = irs - 1; ir < nr; ir++)
  89.                 {
  90.                         r += dr;
  91.                         thet = ggrid.ysa[k] - dth;
  92.  
  93.                         /* loop over theta.  (theta=atan((z+h)/rho)) */
  94.                         for (ith = 0; ith < nth; ith++)
  95.                         {
  96.                                 thet += dth;
  97.                                 rho = r * cosl (thet);
  98.                                 zph = r * sinl (thet);
  99.                                 if (rho < 1.e-7)
  100.                                         rho = 1.e-8;
  101.                                 if (zph < 1.e-7)
  102.                                         zph = 0.;
  103.  
  104.                                 evlua (&erv, &ezv, &erh, &eph);
  105.  
  106.                                 rk = ck2 * r;
  107.                                 con = -CONST1 * r / cmplx (cosl (rk), -sinl (rk));
  108.  
  109.                                 switch (k)
  110.                                 {
  111.                                 case 0:
  112.                                         ggrid.ar1[ir + ith * 11 + 0] = erv * con;
  113.                                         ggrid.ar1[ir + ith * 11 + 110] = ezv * con;
  114.                                         ggrid.ar1[ir + ith * 11 + 220] = erh * con;
  115.                                         ggrid.ar1[ir + ith * 11 + 330] = eph * con;
  116.                                         break;
  117.  
  118.                                 case 1:
  119.                                         ggrid.ar2[ir + ith * 17 + 0] = erv * con;
  120.                                         ggrid.ar2[ir + ith * 17 + 85] = ezv * con;
  121.                                         ggrid.ar2[ir + ith * 17 + 170] = erh * con;
  122.                                         ggrid.ar2[ir + ith * 17 + 255] = eph * con;
  123.                                         break;
  124.  
  125.                                 case 2:
  126.                                         ggrid.ar3[ir + ith * 9 + 0] = erv * con;
  127.                                         ggrid.ar3[ir + ith * 9 + 72] = ezv * con;
  128.                                         ggrid.ar3[ir + ith * 9 + 144] = erh * con;
  129.                                         ggrid.ar3[ir + ith * 9 + 216] = eph * con;
  130.  
  131.                                 } /* switch( k ) */
  132.  
  133.                         } /* for( ith = 0; ith < nth; ith++ ) */
  134.  
  135.                 } /* for( ir = irs-1; ir < nr; ir++; ) */
  136.  
  137.         } /* for( k = 0; k < 3; k++; ) */
  138.  
  139.         /* fill grid 1 for r equal to zero. */
  140.         cl2 = -CONST4 * (ggrid.epscf - 1.) / (ggrid.epscf + 1.);
  141.         cl1 = cl2 / (ggrid.epscf + 1.);
  142.         ezv = ggrid.epscf * cl1;
  143.         thet = -dth;
  144.         nth = ggrid.nya[0];
  145.  
  146.         for (ith = 0; ith < nth; ith++)
  147.         {
  148.                 thet += dth;
  149.                 if ((ith + 1) != nth)
  150.                 {
  151.                         tfac2 = cosl (thet);
  152.                         tfac1 = (1. - sinl (thet)) / tfac2;
  153.                         tfac2 = tfac1 / tfac2;
  154.                         erv = ggrid.epscf * cl1 * tfac1;
  155.                         erh = cl1 * (tfac2 - 1.) + cl2;
  156.                         eph = cl1 * tfac2 - cl2;
  157.                 }
  158.                 else
  159.                 {
  160.                         erv = 0.;
  161.                         erh = cl2 - .5 * cl1;
  162.                         eph = -erh;
  163.                 }
  164.  
  165.                 ggrid.ar1[0 + ith * 11 + 0] = erv;
  166.                 ggrid.ar1[0 + ith * 11 + 110] = ezv;
  167.                 ggrid.ar1[0 + ith * 11 + 220] = erh;
  168.                 ggrid.ar1[0 + ith * 11 + 330] = eph;
  169.         }
  170.  
  171.         secnds (&tim);
  172.         tim -= tst;
  173.  
  174.         return;
  175. }
  176.  
  177. /*-----------------------------------------------------------------------*/
  178.  
  179. /* bessel evaluates the zero-order bessel function */
  180. /* and its derivative for complex argument z. */
  181. void bessel (complex double z, complex double *j0, complex double *j0p)
  182. {
  183.         int k, i, ib, iz, miz;
  184.         static int m[101], init = FALSE;
  185.         static double a1[25], a2[25];
  186.         double tst, zms;
  187.         complex double p0z, p1z, q0z, q1z, zi, zi2, zk, cz, sz, j0x = CPLX_00, j0px = CPLX_00;
  188.  
  189.         /* initialization of constants */
  190.         if (!init)
  191.         {
  192.                 for (k = 1; k <= 25; k++)
  193.                 {
  194.                         i = k - 1;
  195.                         a1[i] = -.25 / (k * k);
  196.                         a2[i] = 1.0 / (k + 1.0);
  197.                 }
  198.  
  199.                 for (i = 1; i <= 101; i++)
  200.                 {
  201.                         tst = 1.0;
  202.                         for (k = 0; k < 24; k++)
  203.                         {
  204.                                 init = k;
  205.                                 tst *= -i * a1[k];
  206.                                 if (tst < 1.0e-6)
  207.                                         break;
  208.                         }
  209.  
  210.                         m[i - 1] = init + 1;
  211.                 } /* for( i = 1; i<= 101; i++ ) */
  212.  
  213.                 init = TRUE;
  214.         } /* if(init == 0) */
  215.  
  216.         zms = z * conj (z);
  217.         if (zms <= 1.e-12)
  218.         {
  219.                 *j0 = CPLX_10;
  220.                 *j0p = -.5 * z;
  221.                 return;
  222.         }
  223.  
  224.         ib = 0;
  225.         if (zms <= 37.21)
  226.         {
  227.                 if (zms > 36.)
  228.                         ib = 1;
  229.  
  230.                 /* series expansion */
  231.                 iz = zms;
  232.                 miz = m[iz];
  233.                 *j0 = CPLX_10;
  234.                 *j0p = *j0;
  235.                 zk = *j0;
  236.                 zi = z * z;
  237.  
  238.                 for (k = 0; k < miz; k++)
  239.                 {
  240.                         zk *= a1[k] * zi;
  241.                         *j0 += zk;
  242.                         *j0p += a2[k] * zk;
  243.                 }
  244.                 *j0p *= -.5 * z;
  245.  
  246.                 if (ib == 0)
  247.                         return;
  248.  
  249.                 j0x = *j0;
  250.                 j0px = *j0p;
  251.         }
  252.  
  253.         /* asymptotic expansion */
  254.         zi = 1. / z;
  255.         zi2 = zi * zi;
  256.         p0z = 1. + (P20 * zi2 - P10) * zi2;
  257.         p1z = 1. + (P11 - P21 * zi2) * zi2;
  258.         q0z = (Q20 * zi2 - Q10) * zi;
  259.         q1z = (Q11 - Q21 * zi2) * zi;
  260.         zk = cexp (CPLX_01 * (z - POF));
  261.         zi2 = 1. / zk;
  262.         cz = .5 * (zk + zi2);
  263.         sz = CPLX_01 * .5 * (zi2 - zk);
  264.         zk = C3 * csqrtl (zi);
  265.         *j0 = zk * (p0z * cz - q0z * sz);
  266.         *j0p = -zk * (p1z * sz + q1z * cz);
  267.  
  268.         if (ib == 0)
  269.                 return;
  270.  
  271.         zms = cosl ((sqrtl (zms) - 6.) * PI10);
  272.         *j0 = .5 * (j0x * (1. + zms) + *j0 * (1. - zms));
  273.         *j0p = .5 * (j0px * (1. + zms) + *j0p * (1. - zms));
  274.  
  275.         return;
  276. }
  277.  
  278. /*-----------------------------------------------------------------------*/
  279.  
  280. /* evlua controls the integration contour in the complex */
  281. /* lambda plane for evaluation of the sommerfeld integrals */
  282. void evlua (complex double *erv, complex double *ezv, complex double *erh, complex double *eph)
  283. {
  284.         int i, jump;
  285.         static double del, slope, rmis;
  286.         static complex double cp1, cp2, cp3, bk, delta, delta2, sum[6], ans[6];
  287.  
  288.         del = zph;
  289.         if (rho > del)
  290.                 del = rho;
  291.  
  292.         if (zph >= 2. * rho)
  293.         {
  294.                 /* bessel function form of sommerfeld integrals */
  295.                 jh = 0;
  296.                 a = CPLX_00;
  297.                 del = 1. / del;
  298.  
  299.                 if (del > tkmag)
  300.                 {
  301.                         b = cmplx (.1 * tkmag, -.1 * tkmag);
  302.                         rom1 (6, sum, 2);
  303.                         a = b;
  304.                         b = cmplx (del, -del);
  305.                         rom1 (6, ans, 2);
  306.                         for (i = 0; i < 6; i++)
  307.                                 sum[i] += ans[i];
  308.                 }
  309.                 else
  310.                 {
  311.                         b = cmplx (del, -del);
  312.                         rom1 (6, sum, 2);
  313.                 }
  314.  
  315.                 delta = PTP * del;
  316.                 gshank (b, delta, ans, 6, sum, 0, b, b);
  317.                 ans[5] *= ck1;
  318.  
  319.                 /* conjugate since nec uses exp(+jwt) */
  320.                 *erv = conj (ck1sq * ans[2]);
  321.                 *ezv = conj (ck1sq * (ans[1] + ck2sq * ans[4]));
  322.                 *erh = conj (ck2sq * (ans[0] + ans[5]));
  323.                 *eph = -conj (ck2sq * (ans[3] + ans[5]));
  324.  
  325.                 return;
  326.  
  327.         } /* if(zph >= 2.*rho) */
  328.  
  329.         /* hankel function form of sommerfeld integrals */
  330.         jh = 1;
  331.         cp1 = cmplx (0.0, .4 * ck2);
  332.         cp2 = cmplx (.6 * ck2, -.2 * ck2);
  333.         cp3 = cmplx (1.02 * ck2, -.2 * ck2);
  334.         a = cp1;
  335.         b = cp2;
  336.         rom1 (6, sum, 2);
  337.         a = cp2;
  338.         b = cp3;
  339.         rom1 (6, ans, 2);
  340.  
  341.         for (i = 0; i < 6; i++)
  342.                 sum[i] = -(sum[i] + ans[i]);
  343.  
  344.         /* path from imaginary axis to -infinity */
  345.         if (zph > .001 * rho)
  346.                 slope = rho / zph;
  347.         else
  348.                 slope = 1000.;
  349.  
  350.         del = PTP / del;
  351.         delta = cmplx (-1.0, slope) * del / sqrtl (1. + slope * slope);
  352.         delta2 = -conj (delta);
  353.         gshank (cp1, delta, ans, 6, sum, 0, bk, bk);
  354.         rmis = rho * (creal (ck1) - ck2);
  355.  
  356.         jump = FALSE;
  357.         if ((rmis >= 2. * ck2) && (rho >= 1.e-10))
  358.         {
  359.                 if (zph >= 1.e-10)
  360.                 {
  361.                         bk = cmplx (-zph, rho) * (ck1 - cp3);
  362.                         rmis = -creal (bk) / fabsl (cimag (bk));
  363.                         if (rmis > 4. * rho / zph)
  364.                                 jump = TRUE;
  365.                 }
  366.  
  367.                 if (!jump)
  368.                 {
  369.                         /* integrate up between branch cuts, then to + infinity */
  370.                         cp1 = ck1 - (.1 + .2fj);
  371.                         cp2 = cp1 + .2;
  372.                         bk = cmplx (0., del);
  373.                         gshank (cp1, bk, sum, 6, ans, 0, bk, bk);
  374.                         a = cp1;
  375.                         b = cp2;
  376.                         rom1 (6, ans, 1);
  377.                         for (i = 0; i < 6; i++)
  378.                                 ans[i] -= sum[i];
  379.  
  380.                         gshank (cp3, bk, sum, 6, ans, 0, bk, bk);
  381.                         gshank (cp2, delta2, ans, 6, sum, 0, bk, bk);
  382.                 }
  383.  
  384.                 jump = TRUE;
  385.  
  386.         } /* if( (rmis >= 2.*ck2) || (rho >= 1.e-10) ) */
  387.         else
  388.                 jump = FALSE;
  389.  
  390.         if (!jump)
  391.         {
  392.                 /* integrate below branch points, then to + infinity */
  393.                 for (i = 0; i < 6; i++)
  394.                         sum[i] = -ans[i];
  395.  
  396.                 rmis = creal (ck1) * 1.01;
  397.                 if ((ck2 + 1.) > rmis)
  398.                         rmis = ck2 + 1.;
  399.  
  400.                 bk = cmplx (rmis, .99 * cimag (ck1));
  401.                 delta = bk - cp3;
  402.                 delta *= del / cabs (delta);
  403.                 gshank (cp3, delta, ans, 6, sum, 1, bk, delta2);
  404.  
  405.         } /* if( ! jump ) */
  406.  
  407.         ans[5] *= ck1;
  408.  
  409.         /* conjugate since nec uses exp(+jwt) */
  410.         *erv = conj (ck1sq * ans[2]);
  411.         *ezv = conj (ck1sq * (ans[1] + ck2sq * ans[4]));
  412.         *erh = conj (ck2sq * (ans[0] + ans[5]));
  413.         *eph = -conj (ck2sq * (ans[3] + ans[5]));
  414.  
  415.         return;
  416. }
  417.  
  418. /*-----------------------------------------------------------------------*/
  419.  
  420. /* fbar is sommerfeld attenuation function for numerical distance p */
  421. void fbar (complex double p, complex double *fbar)
  422. {
  423.         int i, minus;
  424.         double tms, sms;
  425.         complex double z, zs, sum, pow, term;
  426.  
  427.         z = CPLX_01 * csqrtl (p);
  428.         if (cabs (z) <= 3.)
  429.         {
  430.                 /* series expansion */
  431.                 zs = z * z;
  432.                 sum = z;
  433.                 pow = z;
  434.  
  435.                 for (i = 1; i <= 100; i++)
  436.                 {
  437.                         pow = -pow * zs / (double) i;
  438.                         term = pow / (2. * i + 1.);
  439.                         sum = sum + term;
  440.                         tms = creal (term * conj (term));
  441.                         sms = creal (sum * conj (sum));
  442.                         if (tms / sms < ACCS)
  443.                                 break;
  444.                 }
  445.  
  446.                 *fbar = 1. - (1. - sum * TOSP) * z * cexp (zs) * SP;
  447.  
  448.         } /* if( cabs( z) <= 3.) */
  449.  
  450.         /* asymptotic expansion */
  451.         if (creal (z) < 0.)
  452.         {
  453.                 minus = 1;
  454.                 z = -z;
  455.         }
  456.         else
  457.                 minus = 0;
  458.  
  459.         zs = .5 / (z * z);
  460.         sum = CPLX_00;
  461.         term = CPLX_10;
  462.  
  463.         for (i = 1; i <= 6; i++)
  464.         {
  465.                 term = -term * (2. * i - 1.) * zs;
  466.                 sum += term;
  467.         }
  468.  
  469.         if (minus == 1)
  470.                 sum -= 2. * SP * z * cexp (z * z);
  471.         *fbar = -sum;
  472. }
  473.  
  474. /*-----------------------------------------------------------------------*/
  475.  
  476. /* gshank integrates the 6 sommerfeld integrals from start to */
  477. /* infinity (until convergence) in lambda.  at the break point, bk, */
  478. /* the step increment may be changed from dela to delb.  shank's */
  479. /* algorithm to accelerate convergence of a slowly converging series */
  480. /* is used */
  481. void gshank (
  482.     complex double start,
  483.     complex double dela,
  484.     complex double *sum,
  485.     int nans,
  486.     complex double *seed,
  487.     int ibk,
  488.     complex double bk,
  489.     complex double delb)
  490. {
  491.         int ibx, j, i, jm, intx, inx, brk = 0;
  492.         static double rbk, amg, den, denm;
  493.         complex double a1, a2, as1, as2, del, aa;
  494.         complex double q1[6][20], q2[6][20], ans1[6], ans2[6];
  495.  
  496.         rbk = creal (bk);
  497.         del = dela;
  498.         if (ibk == 0)
  499.                 ibx = 1;
  500.         else
  501.                 ibx = 0;
  502.  
  503.         for (i = 0; i < nans; i++)
  504.                 ans2[i] = seed[i];
  505.  
  506.         b = start;
  507.         for (intx = 1; intx <= MAXH; intx++)
  508.         {
  509.                 inx = intx - 1;
  510.                 a = b;
  511.                 b += del;
  512.  
  513.                 if ((ibx == 0) && (creal (b) >= rbk))
  514.                 {
  515.                         /* hit break point.  reset seed and start over. */
  516.                         ibx = 1;
  517.                         b = bk;
  518.                         del = delb;
  519.                         rom1 (nans, sum, 2);
  520.                         if (ibx != 2)
  521.                         {
  522.                                 for (i = 0; i < nans; i++)
  523.                                         ans2[i] += sum[i];
  524.                                 intx = 0;
  525.                                 continue;
  526.                         }
  527.  
  528.                         for (i = 0; i < nans; i++)
  529.                                 ans2[i] = ans1[i] + sum[i];
  530.                         intx = 0;
  531.                         continue;
  532.  
  533.                 } /* if( (ibx == 0) && (creal(b) >= rbk) ) */
  534.  
  535.                 rom1 (nans, sum, 2);
  536.                 for (i = 0; i < nans; i++)
  537.                         ans1[i] = ans2[i] + sum[i];
  538.                 a = b;
  539.                 b += del;
  540.  
  541.                 if ((ibx == 0) && (creal (b) >= rbk))
  542.                 {
  543.                         /* hit break point.  reset seed and start over. */
  544.                         ibx = 2;
  545.                         b = bk;
  546.                         del = delb;
  547.                         rom1 (nans, sum, 2);
  548.                         if (ibx != 2)
  549.                         {
  550.                                 for (i = 0; i < nans; i++)
  551.                                         ans2[i] += sum[i];
  552.                                 intx = 0;
  553.                                 continue;
  554.                         }
  555.  
  556.                         for (i = 0; i < nans; i++)
  557.                                 ans2[i] = ans1[i] + sum[i];
  558.                         intx = 0;
  559.                         continue;
  560.  
  561.                 } /* if( (ibx == 0) && (creal(b) >= rbk) ) */
  562.  
  563.                 rom1 (nans, sum, 2);
  564.                 for (i = 0; i < nans; i++)
  565.                         ans2[i] = ans1[i] + sum[i];
  566.  
  567.                 den = 0.;
  568.                 for (i = 0; i < nans; i++)
  569.                 {
  570.                         as1 = ans1[i];
  571.                         as2 = ans2[i];
  572.  
  573.                         if (intx >= 2)
  574.                         {
  575.                                 for (j = 1; j < intx; j++)
  576.                                 {
  577.                                         jm = j - 1;
  578.                                         aa = q2[i][jm];
  579.                                         a1 = q1[i][jm] + as1 - 2. * aa;
  580.  
  581.                                         if ((creal (a1) != 0.) || (cimag (a1) != 0.))
  582.                                         {
  583.                                                 a2 = aa - q1[i][jm];
  584.                                                 a1 = q1[i][jm] - a2 * a2 / a1;
  585.                                         }
  586.                                         else
  587.                                                 a1 = q1[i][jm];
  588.  
  589.                                         a2 = aa + as2 - 2. * as1;
  590.                                         if ((creal (a2) != 0.) || (cimag (a2) != 0.))
  591.                                                 a2 = aa - (as1 - aa) * (as1 - aa) / a2;
  592.                                         else
  593.                                                 a2 = aa;
  594.  
  595.                                         q1[i][jm] = as1;
  596.                                         q2[i][jm] = as2;
  597.                                         as1 = a1;
  598.                                         as2 = a2;
  599.  
  600.                                 } /* for( j = 1; i < intx; i++ ) */
  601.  
  602.                         } /* if(intx >= 2) */
  603.  
  604.                         q1[i][intx - 1] = as1;
  605.                         q2[i][intx - 1] = as2;
  606.                         amg = fabsl (creal (as2)) + fabsl (cimag (as2));
  607.                         if (amg > den)
  608.                                 den = amg;
  609.  
  610.                 } /* for( i = 0; i < nans; i++ ) */
  611.  
  612.                 denm = 1.e-3 * den * CRIT;
  613.                 jm = intx - 3;
  614.                 if (jm < 1)
  615.                         jm = 1;
  616.  
  617.                 for (j = jm - 1; j < intx; j++)
  618.                 {
  619.                         brk = FALSE;
  620.                         for (i = 0; i < nans; i++)
  621.                         {
  622.                                 a1 = q2[i][j];
  623.                                 den = (fabsl (creal (a1)) + fabsl (cimag (a1))) * CRIT;
  624.                                 if (den < denm)
  625.                                         den = denm;
  626.                                 a1 = q1[i][j] - a1;
  627.                                 amg = fabsl (creal (a1) + fabsl (cimag (a1)));
  628.                                 if (amg > den)
  629.                                 {
  630.                                         brk = TRUE;
  631.                                         break;
  632.                                 }
  633.  
  634.                         } /* for( i = 0; i < nans; i++ ) */
  635.  
  636.                         if (brk)
  637.                                 break;
  638.  
  639.                 } /* for( j = jm-1; j < intx; j++ ) */
  640.  
  641.                 if (!brk)
  642.                 {
  643.                         for (i = 0; i < nans; i++)
  644.                                 sum[i] = .5 * (q1[i][inx] + q2[i][inx]);
  645.                         return;
  646.                 }
  647.  
  648.         } /* for( intx = 1; intx <= maxh; intx++ ) */
  649.  
  650.         /* No convergence */
  651.         abort_on_error (-6);
  652. }
  653.  
  654. /*-----------------------------------------------------------------------*/
  655.  
  656. /* hankel evaluates hankel function of the first kind,   */
  657. /* order zero, and its derivative for complex argument z */
  658. void hankel (complex double z, complex double *h0, complex double *h0p)
  659. {
  660.         int i, k, ib, iz, miz;
  661.         static int m[101], init = FALSE;
  662.         static double a1[25], a2[25], a3[25], a4[25], psi, tst, zms;
  663.         complex double clogz, j0, j0p, p0z, p1z, q0z, q1z, y0 = CPLX_00, y0p = CPLX_00, zi,
  664.                                                            zi2, zk;
  665.  
  666.         /* initialization of constants */
  667.         if (!init)
  668.         {
  669.                 psi = -GAMMA;
  670.                 for (k = 1; k <= 25; k++)
  671.                 {
  672.                         i = k - 1;
  673.                         a1[i] = -.25 / (k * k);
  674.                         a2[i] = 1.0 / (k + 1.0);
  675.                         psi += 1.0 / k;
  676.                         a3[i] = psi + psi;
  677.                         a4[i] = (psi + psi + 1.0 / (k + 1.0)) / (k + 1.0);
  678.                 }
  679.  
  680.                 for (i = 1; i <= 101; i++)
  681.                 {
  682.                         tst = 1.0;
  683.                         for (k = 0; k < 24; k++)
  684.                         {
  685.                                 init = k;
  686.                                 tst *= -i * a1[k];
  687.                                 if (tst * a3[k] < 1.e-6)
  688.                                         break;
  689.                         }
  690.                         m[i - 1] = init + 1;
  691.                 }
  692.  
  693.                 init = TRUE;
  694.  
  695.         } /* if( ! init ) */
  696.  
  697.         zms = z * conj (z);
  698.         if (zms == 0.)
  699.                 abort_on_error (-7);
  700.  
  701.         ib = 0;
  702.         if (zms <= 16.81)
  703.         {
  704.                 if (zms > 16.)
  705.                         ib = 1;
  706.  
  707.                 /* series expansion */
  708.                 iz = zms;
  709.                 miz = m[iz];
  710.                 j0 = CPLX_10;
  711.                 j0p = j0;
  712.                 y0 = CPLX_00;
  713.                 y0p = y0;
  714.                 zk = j0;
  715.                 zi = z * z;
  716.  
  717.                 for (k = 0; k < miz; k++)
  718.                 {
  719.                         zk *= a1[k] * zi;
  720.                         j0 += zk;
  721.                         j0p += a2[k] * zk;
  722.                         y0 += a3[k] * zk;
  723.                         y0p += a4[k] * zk;
  724.                 }
  725.  
  726.                 j0p *= -.5 * z;
  727.                 clogz = clogl (.5 * z);
  728.                 y0 = (2. * j0 * clogz - y0) / PI + C2;
  729.                 y0p = (2. / z + 2. * j0p * clogz + .5 * y0p * z) / PI + C1 * z;
  730.                 *h0 = j0 + CPLX_01 * y0;
  731.                 *h0p = j0p + CPLX_01 * y0p;
  732.  
  733.                 if (ib == 0)
  734.                         return;
  735.  
  736.                 y0 = *h0;
  737.                 y0p = *h0p;
  738.  
  739.         } /* if(zms <= 16.81) */
  740.  
  741.         /* asymptotic expansion */
  742.         zi = 1. / z;
  743.         zi2 = zi * zi;
  744.         p0z = 1. + (P20 * zi2 - P10) * zi2;
  745.         p1z = 1. + (P11 - P21 * zi2) * zi2;
  746.         q0z = (Q20 * zi2 - Q10) * zi;
  747.         q1z = (Q11 - Q21 * zi2) * zi;
  748.         zk = cexp (CPLX_01 * (z - POF)) * csqrtl (zi) * C3;
  749.         *h0 = zk * (p0z + CPLX_01 * q0z);
  750.         *h0p = CPLX_01 * zk * (p1z + CPLX_01 * q1z);
  751.  
  752.         if (ib == 0)
  753.                 return;
  754.  
  755.         zms = cosl ((sqrtl (zms) - 4.) * 31.41592654);
  756.         *h0 = .5 * (y0 * (1. + zms) + *h0 * (1. - zms));
  757.         *h0p = .5 * (y0p * (1. + zms) + *h0p * (1. - zms));
  758.  
  759.         return;
  760. }
  761.  
  762. /*-----------------------------------------------------------------------*/
  763.  
  764. /* compute integration parameter xlam=lambda from parameter t. */
  765. void lambda (double t, complex double *xlam, complex double *dxlam)
  766. {
  767.         *dxlam = b - a;
  768.         *xlam = a + *dxlam * t;
  769.         return;
  770. }
  771.  
  772. /*-----------------------------------------------------------------------*/
  773.  
  774. /* rom1 integrates the 6 sommerfeld integrals from a to b in lambda. */
  775. /* the method of variable interval width romberg integration is used. */
  776. void rom1 (int n, complex double *sum, int nx)
  777. {
  778.         int jump, lstep, nogo, i, ns, nt;
  779.         static double z, ze, s, ep, zend, dz = 0., dzot = 0., tr, ti;
  780.         static complex double t00, t11, t02;
  781.         static complex double g1[6], g2[6], g3[6], g4[6], g5[6], t01[6], t10[6], t20[6];
  782.  
  783.         lstep = 0;
  784.         z = 0.;
  785.         ze = 1.;
  786.         s = 1.;
  787.         ep = s / (1.e4 * NM);
  788.         zend = ze - ep;
  789.         for (i = 0; i < n; i++)
  790.                 sum[i] = CPLX_00;
  791.         ns = nx;
  792.         nt = 0;
  793.         saoa (z, g1);
  794.  
  795.         jump = FALSE;
  796.         while (TRUE)
  797.         {
  798.                 if (!jump)
  799.                 {
  800.                         dz = s / ns;
  801.                         if ((z + dz) > ze)
  802.                         {
  803.                                 dz = ze - z;
  804.                                 if (dz <= ep)
  805.                                         return;
  806.                         }
  807.  
  808.                         dzot = dz * .5;
  809.                         saoa (z + dzot, g3);
  810.                         saoa (z + dz, g5);
  811.  
  812.                 } /* if( ! jump ) */
  813.  
  814.                 nogo = FALSE;
  815.                 for (i = 0; i < n; i++)
  816.                 {
  817.                         t00 = (g1[i] + g5[i]) * dzot;
  818.                         t01[i] = (t00 + dz * g3[i]) * .5;
  819.                         t10[i] = (4. * t01[i] - t00) / 3.;
  820.  
  821.                         /* test convergence of 3 point romberg result */
  822.                         test (
  823.                             creal (t01[i]),
  824.                             creal (t10[i]),
  825.                             &tr,
  826.                             cimag (t01[i]),
  827.                             cimag (t10[i]),
  828.                             &ti,
  829.                             0.);
  830.                         if ((tr > CRIT) || (ti > CRIT))
  831.                                 nogo = TRUE;
  832.                 }
  833.  
  834.                 if (!nogo)
  835.                 {
  836.                         for (i = 0; i < n; i++)
  837.                                 sum[i] += t10[i];
  838.  
  839.                         nt += 2;
  840.                         z += dz;
  841.                         if (z > zend)
  842.                                 return;
  843.  
  844.                         for (i = 0; i < n; i++)
  845.                                 g1[i] = g5[i];
  846.  
  847.                         if ((nt >= NTS) && (ns > nx))
  848.                         {
  849.                                 ns = ns / 2;
  850.                                 nt = 1;
  851.                         }
  852.  
  853.                         jump = FALSE;
  854.                         continue;
  855.  
  856.                 } /* if( ! nogo ) */
  857.  
  858.                 saoa (z + dz * .25, g2);
  859.                 saoa (z + dz * .75, g4);
  860.                 nogo = FALSE;
  861.                 for (i = 0; i < n; i++)
  862.                 {
  863.                         t02 = (t01[i] + dzot * (g2[i] + g4[i])) * .5;
  864.                         t11 = (4. * t02 - t01[i]) / 3.;
  865.                         t20[i] = (16. * t11 - t10[i]) / 15.;
  866.  
  867.                         /* test convergence of 5 point romberg result */
  868.                         test (
  869.                             creal (t11),
  870.                             creal (t20[i]),
  871.                             &tr,
  872.                             cimag (t11),
  873.                             cimag (t20[i]),
  874.                             &ti,
  875.                             0.);
  876.                         if ((tr > CRIT) || (ti > CRIT))
  877.                                 nogo = TRUE;
  878.                 }
  879.  
  880.                 if (!nogo)
  881.                 {
  882.                         for (i = 0; i < n; i++)
  883.                                 sum[i] += t20[i];
  884.  
  885.                         nt++;
  886.                         z += dz;
  887.                         if (z > zend)
  888.                                 return;
  889.  
  890.                         for (i = 0; i < n; i++)
  891.                                 g1[i] = g5[i];
  892.  
  893.                         if ((nt >= NTS) && (ns > nx))
  894.                         {
  895.                                 ns = ns / 2;
  896.                                 nt = 1;
  897.                         }
  898.  
  899.                         jump = FALSE;
  900.                         continue;
  901.  
  902.                 } /* if( ! nogo ) */
  903.  
  904.                 nt = 0;
  905.                 if (ns < NM)
  906.                 {
  907.                         ns *= 2;
  908.                         dz = s / ns;
  909.                         dzot = dz * .5;
  910.  
  911.                         for (i = 0; i < n; i++)
  912.                         {
  913.                                 g5[i] = g3[i];
  914.                                 g3[i] = g2[i];
  915.                         }
  916.  
  917.                         jump = TRUE;
  918.                         continue;
  919.  
  920.                 } /* if(ns < nm) */
  921.  
  922.                 if (!lstep)
  923.                 {
  924.                         lstep = TRUE;
  925.                         lambda (z, &t00, &t11);
  926.                 }
  927.  
  928.                 for (i = 0; i < n; i++)
  929.                         sum[i] += t20[i];
  930.  
  931.                 nt++;
  932.                 z += dz;
  933.                 if (z > zend)
  934.                         return;
  935.  
  936.                 for (i = 0; i < n; i++)
  937.                         g1[i] = g5[i];
  938.  
  939.                 if ((nt >= NTS) && (ns > nx))
  940.                 {
  941.                         ns /= 2;
  942.                         nt = 1;
  943.                 }
  944.  
  945.                 jump = FALSE;
  946.  
  947.         } /* while( TRUE ) */
  948. }
  949.  
  950. /*-----------------------------------------------------------------------*/
  951.  
  952. /* saoa computes the integrand for each of the 6 sommerfeld */
  953. /* integrals for source and observer above ground */
  954. void saoa (double t, complex double *ans)
  955. {
  956.         double xlr, sign;
  957.         static complex double xl, dxl, cgam1, cgam2, b0, b0p, com, dgam, den1, den2;
  958.  
  959.         lambda (t, &xl, &dxl);
  960.         if (jh == 0)
  961.         {
  962.                 /* bessel function form */
  963.                 bessel (xl * rho, &b0, &b0p);
  964.                 b0 *= 2.;
  965.                 b0p *= 2.;
  966.                 cgam1 = csqrtl (xl * xl - ck1sq);
  967.                 cgam2 = csqrtl (xl * xl - ck2sq);
  968.                 if (creal (cgam1) == 0.)
  969.                         cgam1 = cmplx (0., -fabsl (cimag (cgam1)));
  970.                 if (creal (cgam2) == 0.)
  971.                         cgam2 = cmplx (0., -fabsl (cimag (cgam2)));
  972.         }
  973.         else
  974.         {
  975.                 /* hankel function form */
  976.                 hankel (xl * rho, &b0, &b0p);
  977.                 com = xl - ck1;
  978.                 cgam1 = csqrtl (xl + ck1) * csqrtl (com);
  979.                 if (creal (com) < 0. && cimag (com) >= 0.)
  980.                         cgam1 = -cgam1;
  981.                 com = xl - ck2;
  982.                 cgam2 = csqrtl (xl + ck2) * csqrtl (com);
  983.                 if (creal (com) < 0. && cimag (com) >= 0.)
  984.                         cgam2 = -cgam2;
  985.         }
  986.  
  987.         xlr = xl * conj (xl);
  988.         if (xlr >= tsmag)
  989.         {
  990.                 if (cimag (xl) >= 0.)
  991.                 {
  992.                         xlr = creal (xl);
  993.                         if (xlr >= ck2)
  994.                         {
  995.                                 if (xlr <= ck1r)
  996.                                         dgam = cgam2 - cgam1;
  997.                                 else
  998.                                 {
  999.                                         sign = 1.;
  1000.                                         dgam = 1. / (xl * xl);
  1001.                                         dgam = sign * ((ct3 * dgam + ct2) * dgam + ct1) / xl;
  1002.                                 }
  1003.                         }
  1004.                         else
  1005.                         {
  1006.                                 sign = -1.;
  1007.                                 dgam = 1. / (xl * xl);
  1008.                                 dgam = sign * ((ct3 * dgam + ct2) * dgam + ct1) / xl;
  1009.                         } /* if(xlr >= ck2) */
  1010.  
  1011.                 } /* if(cimag(xl) >= 0.) */
  1012.                 else
  1013.                 {
  1014.                         sign = 1.;
  1015.                         dgam = 1. / (xl * xl);
  1016.                         dgam = sign * ((ct3 * dgam + ct2) * dgam + ct1) / xl;
  1017.                 }
  1018.  
  1019.         } /* if(xlr < tsmag) */
  1020.         else
  1021.                 dgam = cgam2 - cgam1;
  1022.  
  1023.         den2 = cksm * dgam / (cgam2 * (ck1sq * cgam2 + ck2sq * cgam1));
  1024.         den1 = 1. / (cgam1 + cgam2) - cksm / cgam2;
  1025.         com = dxl * xl * cexp (-cgam2 * zph);
  1026.         ans[5] = com * b0 * den1 / ck1;
  1027.         com *= den2;
  1028.  
  1029.         if (rho != 0.)
  1030.         {
  1031.                 b0p = b0p / rho;
  1032.                 ans[0] = -com * xl * (b0p + b0 * xl);
  1033.                 ans[3] = com * xl * b0p;
  1034.         }
  1035.         else
  1036.         {
  1037.                 ans[0] = -com * xl * xl * .5;
  1038.                 ans[3] = ans[0];
  1039.         }
  1040.  
  1041.         ans[1] = com * cgam2 * cgam2 * b0;
  1042.         ans[2] = -ans[3] * cgam2 * rho;
  1043.         ans[4] = com * b0;
  1044.  
  1045.         return;
  1046. }
  1047.