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  /tmi/ */
  29. extern tmi_t tmi;
  30.  
  31. /*common  /ggrid/ */
  32. extern ggrid_t ggrid;
  33.  
  34. /* common  /data/ */
  35. extern data_t data;
  36.  
  37. /* common  /crnt/ */
  38. extern crnt_t crnt;
  39.  
  40. /* common  /vsorc/ */
  41. extern vsorc_t vsorc;
  42.  
  43. /* common  /segj/ */
  44. extern segj_t segj;
  45.  
  46. /* common  /yparm/ */
  47. extern yparm_t yparm;
  48.  
  49. /* common  /zload/ */
  50. extern zload_t zload;
  51.  
  52. /* common  /smat/ */
  53. extern smat_t smat;
  54.  
  55. /* pointers to input/output files */
  56. extern FILE *input_fp, *output_fp, *plot_fp;
  57.  
  58. /*-----------------------------------------------------------------------*/
  59.  
  60. /* cabc computes coefficients of the constant (a), sine (b), and */
  61. /* cosine (c) terms in the current interpolation functions for the */
  62. /* current vector cur. */
  63. void cabc (complex double *curx)
  64. {
  65.         int i, is, j, jx, jco1, jco2;
  66.         double ar, ai, sh;
  67.         complex double curd, cs1, cs2;
  68.  
  69.         if (data.n != 0)
  70.         {
  71.                 for (i = 0; i < data.n; i++)
  72.                 {
  73.                         crnt.air[i] = 0.;
  74.                         crnt.aii[i] = 0.;
  75.                         crnt.bir[i] = 0.;
  76.                         crnt.bii[i] = 0.;
  77.                         crnt.cir[i] = 0.;
  78.                         crnt.cii[i] = 0.;
  79.                 }
  80.  
  81.                 for (i = 0; i < data.n; i++)
  82.                 {
  83.                         ar = creall (curx[i]);
  84.                         ai = cimagl (curx[i]);
  85.                         tbf (i + 1, 1);
  86.  
  87.                         for (jx = 0; jx < segj.jsno; jx++)
  88.                         {
  89.                                 j = segj.jco[jx] - 1;
  90.                                 crnt.air[j] += segj.ax[jx] * ar;
  91.                                 crnt.aii[j] += segj.ax[jx] * ai;
  92.                                 crnt.bir[j] += segj.bx[jx] * ar;
  93.                                 crnt.bii[j] += segj.bx[jx] * ai;
  94.                                 crnt.cir[j] += segj.cx[jx] * ar;
  95.                                 crnt.cii[j] += segj.cx[jx] * ai;
  96.                         }
  97.  
  98.                 } /* for( i = 0; i < n; i++ ) */
  99.  
  100.                 if (vsorc.nqds != 0)
  101.                 {
  102.                         for (is = 0; is < vsorc.nqds; is++)
  103.                         {
  104.                                 i = vsorc.iqds[is] - 1;
  105.                                 jx = data.icon1[i];
  106.                                 data.icon1[i] = 0;
  107.                                 tbf (i + 1, 0);
  108.                                 data.icon1[i] = jx;
  109.                                 sh = data.si[i] * .5;
  110.                                 curd = CCJ * vsorc.vqds[is] /
  111.                                        ((logl (2. * sh / data.bi[i]) - 1.) *
  112.                                         (segj.bx[segj.jsno - 1] * cosl (TP * sh) +
  113.                                          segj.cx[segj.jsno - 1] * sinl (TP * sh)) *
  114.                                         data.wlam);
  115.                                 ar = print_creall (curd);
  116.                                 ai = print_cimagl (curd);
  117.  
  118.                                 for (jx = 0; jx < segj.jsno; jx++)
  119.                                 {
  120.                                         j = segj.jco[jx] - 1;
  121.                                         crnt.air[j] = crnt.air[j] + segj.ax[jx] * ar;
  122.                                         crnt.aii[j] = crnt.aii[j] + segj.ax[jx] * ai;
  123.                                         crnt.bir[j] = crnt.bir[j] + segj.bx[jx] * ar;
  124.                                         crnt.bii[j] = crnt.bii[j] + segj.bx[jx] * ai;
  125.                                         crnt.cir[j] = crnt.cir[j] + segj.cx[jx] * ar;
  126.                                         crnt.cii[j] = crnt.cii[j] + segj.cx[jx] * ai;
  127.                                 }
  128.  
  129.                         } /* for( is = 0; is < vsorc.nqds; is++ ) */
  130.  
  131.                 } /* if( vsorc.nqds != 0) */
  132.  
  133.                 for (i = 0; i < data.n; i++)
  134.                         curx[i] = cmplx (crnt.air[i] + crnt.cir[i], crnt.aii[i] + crnt.cii[i]);
  135.  
  136.         } /* if( n != 0) */
  137.  
  138.         if (data.m == 0)
  139.                 return;
  140.  
  141.         /* convert surface currents from */
  142.         /* t1,t2 components to x,y,z components */
  143.         jco1 = data.np2m;
  144.         jco2 = jco1 + data.m;
  145.         for (i = 1; i <= data.m; i++)
  146.         {
  147.                 jco1 -= 2;
  148.                 jco2 -= 3;
  149.                 cs1 = curx[jco1];
  150.                 cs2 = curx[jco1 + 1];
  151.                 curx[jco2] = cs1 * data.t1x[data.m - i] + cs2 * data.t2x[data.m - i];
  152.                 curx[jco2 + 1] = cs1 * data.t1y[data.m - i] + cs2 * data.t2y[data.m - i];
  153.                 curx[jco2 + 2] = cs1 * data.t1z[data.m - i] + cs2 * data.t2z[data.m - i];
  154.         }
  155.  
  156.         return;
  157. }
  158.  
  159. /*-----------------------------------------------------------------------*/
  160.  
  161. /* couple computes the maximum coupling between pairs of segments. */
  162. void couple (complex double *cur, double wlam)
  163. {
  164.         int j, j1, j2, l1, i, k, itt1, itt2, its1, its2, isg1, isg2, npm1;
  165.         double dbc, c, gmax;
  166.         complex double y11, y12, y22, yl, yin, zl, zin, rho;
  167.  
  168.         if ((vsorc.nsant != 1) || (vsorc.nvqd != 0))
  169.                 return;
  170.  
  171.         j = isegno (yparm.nctag[yparm.icoup], yparm.ncseg[yparm.icoup]);
  172.         if (j != vsorc.isant[0])
  173.                 return;
  174.  
  175.         zin = vsorc.vsant[0];
  176.         yparm.icoup++;
  177.         mem_realloc ((void *) &yparm.y11a, yparm.icoup * sizeof (complex double));
  178.         yparm.y11a[yparm.icoup - 1] = cur[j - 1] * wlam / zin;
  179.  
  180.         l1 = (yparm.icoup - 1) * (yparm.ncoup - 1);
  181.         for (i = 0; i < yparm.ncoup; i++)
  182.         {
  183.                 if ((i + 1) == yparm.icoup)
  184.                         continue;
  185.  
  186.                 l1++;
  187.                 mem_realloc ((void *) &yparm.y12a, l1 * sizeof (complex double));
  188.                 k = isegno (yparm.nctag[i], yparm.ncseg[i]);
  189.                 yparm.y12a[l1 - 1] = cur[k - 1] * wlam / zin;
  190.         }
  191.  
  192.         if (yparm.icoup < yparm.ncoup)
  193.                 return;
  194.  
  195.         fprintf (
  196.             output_fp,
  197.             "\n\n\n"
  198.             "                        -----------"
  199.             " ISOLATION DATA -----------\n\n"
  200.             " ------- COUPLING BETWEEN ------     MAXIMUM    "
  201.             " ---------- FOR MAXIMUM COUPLING ----------\n"
  202.             "            SEG              SEG    COUPLING  LOAD"
  203.             " IMPEDANCE (2ND SEG)         INPUT IMPEDANCE \n"
  204.             " TAG  SEG   NO.   TAG  SEG   NO.      (DB)       "
  205.             " REAL     IMAGINARY         REAL       IMAGINARY");
  206.  
  207.         npm1 = yparm.ncoup - 1;
  208.  
  209.         for (i = 0; i < npm1; i++)
  210.         {
  211.                 itt1 = yparm.nctag[i];
  212.                 its1 = yparm.ncseg[i];
  213.                 isg1 = isegno (itt1, its1);
  214.                 l1 = i + 1;
  215.  
  216.                 for (j = l1; j < yparm.ncoup; j++)
  217.                 {
  218.                         itt2 = yparm.nctag[j];
  219.                         its2 = yparm.ncseg[j];
  220.                         isg2 = isegno (itt2, its2);
  221.                         j1 = j + i * npm1 - 1;
  222.                         j2 = i + j * npm1;
  223.                         y11 = yparm.y11a[i];
  224.                         y22 = yparm.y11a[j];
  225.                         y12 = .5 * (yparm.y12a[j1] + yparm.y12a[j2]);
  226.                         yin = y12 * y12;
  227.                         dbc = cabsl (yin);
  228.                         c = dbc / (2. * creall (y11) * creall (y22) - creall (yin));
  229.  
  230.                         if ((c >= 0.0) && (c <= 1.0))
  231.                         {
  232.                                 if (c >= .01)
  233.                                         gmax = (1. - sqrtl (1. - c * c)) / c;
  234.                                 else
  235.                                         gmax = .5 * (c + .25 * c * c * c);
  236.  
  237.                                 rho = gmax * conjl (yin) / dbc;
  238.                                 yl = ((1. - rho) / (1. + rho) + 1.) * creall (y22) - y22;
  239.                                 zl = 1. / yl;
  240.                                 yin = y11 - yin / (y22 + yl);
  241.                                 zin = 1. / yin;
  242.                                 dbc = db10 (gmax);
  243.  
  244.                                 fprintf (
  245.                                     output_fp,
  246.                                     "\n"
  247.                                     " %4d %4d %5d  %4d %4d %5d  %9.3f"
  248.                                     "  %12.5f %12.5f  %12.5f %12.5f",
  249.                                     itt1,
  250.                                     its1,
  251.                                     isg1,
  252.                                     itt2,
  253.                                     its2,
  254.                                     isg2,
  255.                                     dbc,
  256.                                     print_creall (zl),
  257.                                     print_cimagl (zl),
  258.                                     print_creall (zin),
  259.                                     print_cimagl (zin));
  260.  
  261.                                 continue;
  262.  
  263.                         } /* if( (c >= 0.0) && (c <= 1.0) ) */
  264.  
  265.                         fprintf (
  266.                             output_fp,
  267.                             "\n"
  268.                             " %4d %4d %5d   %4d %4d %5d  **ERROR** "
  269.                             "COUPLING IS NOT BETWEEN 0 AND 1. (= %12.5f)",
  270.                             itt1,
  271.                             its1,
  272.                             isg1,
  273.                             itt2,
  274.                             its2,
  275.                             isg2,
  276.                             c);
  277.  
  278.                 } /* for( j = l1; j < yparm.ncoup; j++ ) */
  279.  
  280.         } /* for( i = 0; i < npm1; i++ ) */
  281.  
  282.         return;
  283. }
  284.  
  285. /*-----------------------------------------------------------------------*/
  286.  
  287. /* load calculates the impedance of specified */
  288. /* segments for various types of loading */
  289. void load (
  290.     int *ldtyp, int *ldtag, int *ldtagf, int *ldtagt, double *zlr, double *zli, double *zlc)
  291. {
  292.         int i, iwarn, istep, istepx, l1, l2, ldtags, jump, ichk;
  293.         complex double zt = CPLX_00, tpcj;
  294.  
  295.         tpcj = (0.0 + 1.8836989e-9j);
  296.         fprintf (
  297.             output_fp,
  298.             "\n"
  299.             "  LOCATION        RESISTANCE  INDUCTANCE  CAPACITANCE   "
  300.             "  IMPEDANCE (OHMS)   CONDUCTIVITY  CIRCUIT\n"
  301.             "       OHMS       HENRYS      FARADS     "
  302.             "  REAL     IMAGINARY   MHOS/METER      TYPE");
  303.  
  304.         /* initialize d array, used for temporary */
  305.         /* storage of loading information. */
  306.         mem_realloc ((void *) &zload.zarray, data.npm * sizeof (complex double));
  307.         for (i = 0; i < data.n; i++)
  308.                 zload.zarray[i] = CPLX_00;
  309.  
  310.         iwarn = FALSE;
  311.         istep = 0;
  312.  
  313.         /* cycle over loading cards */
  314.         while (TRUE)
  315.         {
  316.                 istepx = istep;
  317.                 istep++;
  318.  
  319.                 if (istep > zload.nload)
  320.                 {
  321.                         if (iwarn == TRUE)
  322.                                 fprintf (
  323.                                     output_fp,
  324.                                     "\n  NOTE, SOME OF THE ABOVE SEGMENTS "
  325.                                     "HAVE BEEN LOADED TWICE - IMPEDANCES ADDED");
  326.  
  327.                         smat.nop = data.n / data.np;
  328.                         if (smat.nop == 1)
  329.                                 return;
  330.  
  331.                         for (i = 0; i < data.np; i++)
  332.                         {
  333.                                 zt = zload.zarray[i];
  334.                                 l1 = i;
  335.  
  336.                                 for (l2 = 1; l2 < smat.nop; l2++)
  337.                                 {
  338.                                         l1 += data.np;
  339.                                         zload.zarray[l1] = zt;
  340.                                 }
  341.                         }
  342.                         return;
  343.  
  344.                 } /* if( istep > zload.nload) */
  345.  
  346.                 if (ldtyp[istepx] > 5)
  347.                 {
  348.                         fprintf (
  349.                             output_fp,
  350.                             "\n  IMPROPER LOAD TYPE CHOSEN,"
  351.                             " REQUESTED TYPE IS %d",
  352.                             ldtyp[istepx]);
  353.                         stop (-1);
  354.                 }
  355.  
  356.                 /* search segments for proper itags */
  357.                 ldtags = ldtag[istepx];
  358.                 jump = ldtyp[istepx] + 1;
  359.                 ichk = 0;
  360.                 l1 = 1;
  361.                 l2 = data.n;
  362.  
  363.                 if (ldtags == 0)
  364.                 {
  365.                         if ((ldtagf[istepx] != 0) || (ldtagt[istepx] != 0))
  366.                         {
  367.                                 l1 = ldtagf[istepx];
  368.                                 l2 = ldtagt[istepx];
  369.  
  370.                         } /* if( (ldtagf[istepx] != 0) || (ldtagt[istepx] != 0) ) */
  371.  
  372.                 } /* if( ldtags == 0) */
  373.  
  374.                 for (i = l1 - 1; i < l2; i++)
  375.                 {
  376.                         if (ldtags != 0)
  377.                         {
  378.                                 if (ldtags != data.itag[i])
  379.                                         continue;
  380.  
  381.                                 if (ldtagf[istepx] != 0)
  382.                                 {
  383.                                         ichk++;
  384.                                         if ((ichk < ldtagf[istepx]) || (ichk > ldtagt[istepx]))
  385.                                                 continue;
  386.                                 }
  387.                                 else
  388.                                         ichk = 1;
  389.  
  390.                         } /* if( ldtags != 0) */
  391.                         else
  392.                                 ichk = 1;
  393.  
  394.                         /* calculation of lamda*imped. per unit length, */
  395.                         /* jump to appropriate section for loading type */
  396.                         switch (jump)
  397.                         {
  398.                         case 1:
  399.                                 zt = zlr[istepx] / data.si[i] +
  400.                                      tpcj * zli[istepx] / (data.si[i] * data.wlam);
  401.                                 if (fabsl (zlc[istepx]) > 1.0e-20)
  402.                                         zt += data.wlam / (tpcj * data.si[i] * zlc[istepx]);
  403.                                 break;
  404.  
  405.                         case 2:
  406.                                 zt = tpcj * data.si[i] * zlc[istepx] / data.wlam;
  407.                                 if (fabsl (zli[istepx]) > 1.0e-20)
  408.                                         zt += data.si[i] * data.wlam / (tpcj * zli[istepx]);
  409.                                 if (fabsl (zlr[istepx]) > 1.0e-20)
  410.                                         zt += data.si[i] / zlr[istepx];
  411.                                 zt = 1. / zt;
  412.                                 break;
  413.  
  414.                         case 3:
  415.                                 zt = zlr[istepx] * data.wlam + tpcj * zli[istepx];
  416.                                 if (fabsl (zlc[istepx]) > 1.0e-20)
  417.                                         zt += 1. /
  418.                                               (tpcj * data.si[i] * data.si[i] * zlc[istepx]);
  419.                                 break;
  420.  
  421.                         case 4:
  422.                                 zt = tpcj * data.si[i] * data.si[i] * zlc[istepx];
  423.                                 if (fabsl (zli[istepx]) > 1.0e-20)
  424.                                         zt += 1. / (tpcj * zli[istepx]);
  425.                                 if (fabsl (zlr[istepx]) > 1.0e-20)
  426.                                         zt += 1. / (zlr[istepx] * data.wlam);
  427.                                 zt = 1. / zt;
  428.                                 break;
  429.  
  430.                         case 5:
  431.                                 zt = cmplx (zlr[istepx], zli[istepx]) / data.si[i];
  432.                                 break;
  433.  
  434.                         case 6:
  435.                                 zint (zlr[istepx] * data.wlam, data.bi[i], &zt);
  436.  
  437.                         } /* switch( jump ) */
  438.  
  439.                         if ((fabsl (creall (zload.zarray[i])) +
  440.                              fabsl (cimagl (zload.zarray[i]))) > 1.0e-20)
  441.                                 iwarn = TRUE;
  442.                         zload.zarray[i] += zt;
  443.  
  444.                 } /* for( i = l1-1; i < l2; i++ ) */
  445.  
  446.                 if (ichk == 0)
  447.                 {
  448.                         fprintf (
  449.                             output_fp,
  450.                             "\n  LOADING DATA CARD ERROR,"
  451.                             " NO SEGMENT HAS AN ITAG = %d",
  452.                             ldtags);
  453.                         stop (-1);
  454.                 }
  455.  
  456.                 /* printing the segment loading data, jump to proper print */
  457.                 switch (jump)
  458.                 {
  459.                 case 1:
  460.                         prnt (
  461.                             ldtags,
  462.                             ldtagf[istepx],
  463.                             ldtagt[istepx],
  464.                             zlr[istepx],
  465.                             zli[istepx],
  466.                             zlc[istepx],
  467.                             0.,
  468.                             0.,
  469.                             0.,
  470.                             " SERIES ",
  471.                             2);
  472.                         break;
  473.  
  474.                 case 2:
  475.                         prnt (
  476.                             ldtags,
  477.                             ldtagf[istepx],
  478.                             ldtagt[istepx],
  479.                             zlr[istepx],
  480.                             zli[istepx],
  481.                             zlc[istepx],
  482.                             0.,
  483.                             0.,
  484.                             0.,
  485.                             "PARALLEL",
  486.                             2);
  487.                         break;
  488.  
  489.                 case 3:
  490.                         prnt (
  491.                             ldtags,
  492.                             ldtagf[istepx],
  493.                             ldtagt[istepx],
  494.                             zlr[istepx],
  495.                             zli[istepx],
  496.                             zlc[istepx],
  497.                             0.,
  498.                             0.,
  499.                             0.,
  500.                             "SERIES (PER METER)",
  501.                             5);
  502.                         break;
  503.  
  504.                 case 4:
  505.                         prnt (
  506.                             ldtags,
  507.                             ldtagf[istepx],
  508.                             ldtagt[istepx],
  509.                             zlr[istepx],
  510.                             zli[istepx],
  511.                             zlc[istepx],
  512.                             0.,
  513.                             0.,
  514.                             0.,
  515.                             "PARALLEL (PER METER)",
  516.                             5);
  517.                         break;
  518.  
  519.                 case 5:
  520.                         prnt (
  521.                             ldtags,
  522.                             ldtagf[istepx],
  523.                             ldtagt[istepx],
  524.                             0.,
  525.                             0.,
  526.                             0.,
  527.                             zlr[istepx],
  528.                             zli[istepx],
  529.                             0.,
  530.                             "FIXED IMPEDANCE ",
  531.                             4);
  532.                         break;
  533.  
  534.                 case 6:
  535.                         prnt (
  536.                             ldtags,
  537.                             ldtagf[istepx],
  538.                             ldtagt[istepx],
  539.                             0.,
  540.                             0.,
  541.                             0.,
  542.                             0.,
  543.                             0.,
  544.                             zlr[istepx],
  545.                             "  WIRE  ",
  546.                             2);
  547.  
  548.                 } /* switch( jump ) */
  549.  
  550.         } /* while( TRUE ) */
  551.  
  552.         return;
  553. }
  554.  
  555. /*-----------------------------------------------------------------------*/
  556.  
  557. /* gf computes the integrand exp(jkr)/(kr) for numerical integration. */
  558. void gf (double zk, double *co, double *si)
  559. {
  560.         double zdk, rk, rks;
  561.  
  562.         zdk = zk - tmi.zpk;
  563.         rk = sqrtl (tmi.rkb2 + zdk * zdk);
  564.         *si = sinl (rk) / rk;
  565.  
  566.         if (tmi.ij != 0)
  567.         {
  568.                 *co = cosl (rk) / rk;
  569.                 return;
  570.         }
  571.  
  572.         if (rk >= .2)
  573.         {
  574.                 *co = (cosl (rk) - 1.) / rk;
  575.                 return;
  576.         }
  577.  
  578.         rks = rk * rk;
  579.         *co = ((-1.38888889e-3 * rks + 4.16666667e-2) * rks - .5) * rk;
  580.  
  581.         return;
  582. }
  583.  
  584. /*-----------------------------------------------------------------------*/
  585.  
  586. /* function db10 returns db for magnitude (field) */
  587. double db10 (double x)
  588. {
  589.         if (x < 1.e-20)
  590.                 return (-999.99);
  591.  
  592.         return (10. * log10l (x));
  593. }
  594.  
  595. /*-----------------------------------------------------------------------*/
  596.  
  597. /* function db20 returns db for mag**2 (power) i */
  598. double db20 (double x)
  599. {
  600.         if (x < 1.e-20)
  601.                 return (-999.99);
  602.  
  603.         return (20. * log10l (x));
  604. }
  605.  
  606. /*-----------------------------------------------------------------------*/
  607.  
  608. /* intrp uses bivariate cubic interpolation to obtain */
  609. /* the values of 4 functions at the point (x,y). */
  610. void intrp (
  611.     double x,
  612.     double y,
  613.     complex double *f1,
  614.     complex double *f2,
  615.     complex double *f3,
  616.     complex double *f4)
  617. {
  618.         static int ix, iy, ixs = -10, iys = -10, igrs = -10, ixeg = 0, iyeg = 0;
  619.         static int nxm2, nym2, nxms, nyms, nd, ndp;
  620.         int nda[3] = {11, 17, 9}, ndpa[3] = {110, 85, 72};
  621.         int igr, iadd, iadz, i, k, jump;
  622.         static double dx = 1., dy = 1., xs = 0., ys = 0., xz, yz;
  623.         double xx, yy;
  624.         static complex double a[4][4], b[4][4], c[4][4], d[4][4];
  625.         complex double p1 = CPLX_00, p2 = CPLX_00, p3 = CPLX_00, p4 = CPLX_00;
  626.         complex double fx1, fx2, fx3, fx4;
  627.  
  628.         jump = FALSE;
  629.         if ((x < xs) || (y < ys))
  630.                 jump = TRUE;
  631.         else
  632.         {
  633.                 ix = (int) ((x - xs) / dx) + 1;
  634.                 iy = (int) ((y - ys) / dy) + 1;
  635.         }
  636.  
  637.         /* if point lies in same 4 by 4 point region */
  638.         /* as previous point, old values are reused. */
  639.         if ((ix < ixeg) || (iy < iyeg) || (abs (ix - ixs) >= 2) || (abs (iy - iys) >= 2) ||
  640.             jump)
  641.         {
  642.                 /* determine correct grid and grid region */
  643.                 if (x <= ggrid.xsa[1])
  644.                         igr = 0;
  645.                 else
  646.                 {
  647.                         if (y > ggrid.ysa[2])
  648.                                 igr = 2;
  649.                         else
  650.                                 igr = 1;
  651.                 }
  652.  
  653.                 if (igr != igrs)
  654.                 {
  655.                         igrs = igr;
  656.                         dx = ggrid.dxa[igrs];
  657.                         dy = ggrid.dya[igrs];
  658.                         xs = ggrid.xsa[igrs];
  659.                         ys = ggrid.ysa[igrs];
  660.                         nxm2 = ggrid.nxa[igrs] - 2;
  661.                         nym2 = ggrid.nya[igrs] - 2;
  662.                         nxms = ((nxm2 + 1) / 3) * 3 + 1;
  663.                         nyms = ((nym2 + 1) / 3) * 3 + 1;
  664.                         nd = nda[igrs];
  665.                         ndp = ndpa[igrs];
  666.                         ix = (int) ((x - xs) / dx) + 1;
  667.                         iy = (int) ((y - ys) / dy) + 1;
  668.  
  669.                 } /* if( igr != igrs) */
  670.  
  671.                 ixs = ((ix - 1) / 3) * 3 + 2;
  672.                 if (ixs < 2)
  673.                         ixs = 2;
  674.                 ixeg = -10000;
  675.  
  676.                 if (ixs > nxm2)
  677.                 {
  678.                         ixs = nxm2;
  679.                         ixeg = nxms;
  680.                 }
  681.  
  682.                 iys = ((iy - 1) / 3) * 3 + 2;
  683.                 if (iys < 2)
  684.                         iys = 2;
  685.                 iyeg = -10000;
  686.  
  687.                 if (iys > nym2)
  688.                 {
  689.                         iys = nym2;
  690.                         iyeg = nyms;
  691.                 }
  692.  
  693.                 /* compute coefficients of 4 cubic polynomials in x for */
  694.                 /* the 4 grid values of y for each of the 4 functions */
  695.                 iadz = ixs + (iys - 3) * nd - ndp;
  696.                 for (k = 0; k < 4; k++)
  697.                 {
  698.                         iadz += ndp;
  699.                         iadd = iadz;
  700.  
  701.                         for (i = 0; i < 4; i++)
  702.                         {
  703.                                 iadd += nd;
  704.  
  705.                                 switch (igrs)
  706.                                 {
  707.                                 case 0:
  708.                                         p1 = ggrid.ar1[iadd - 2];
  709.                                         p2 = ggrid.ar1[iadd - 1];
  710.                                         p3 = ggrid.ar1[iadd];
  711.                                         p4 = ggrid.ar1[iadd + 1];
  712.                                         break;
  713.  
  714.                                 case 1:
  715.                                         p1 = ggrid.ar2[iadd - 2];
  716.                                         p2 = ggrid.ar2[iadd - 1];
  717.                                         p3 = ggrid.ar2[iadd];
  718.                                         p4 = ggrid.ar2[iadd + 1];
  719.                                         break;
  720.  
  721.                                 case 2:
  722.                                         p1 = ggrid.ar3[iadd - 2];
  723.                                         p2 = ggrid.ar3[iadd - 1];
  724.                                         p3 = ggrid.ar3[iadd];
  725.                                         p4 = ggrid.ar3[iadd + 1];
  726.  
  727.                                 } /* switch( igrs ) */
  728.  
  729.                                 a[i][k] = (p4 - p1 + 3. * (p2 - p3)) * .1666666667;
  730.                                 b[i][k] = (p1 - 2. * p2 + p3) * .5;
  731.                                 c[i][k] = p3 - (2. * p1 + 3. * p2 + p4) * .1666666667;
  732.                                 d[i][k] = p2;
  733.  
  734.                         } /* for( i = 0; i < 4; i++ ) */
  735.  
  736.                 } /* for( k = 0; k < 4; k++ ) */
  737.  
  738.                 xz = (ixs - 1) * dx + xs;
  739.                 yz = (iys - 1) * dy + ys;
  740.  
  741.         } /* if( (abs(ix- ixs) >= 2) || */
  742.  
  743.         /* evaluate polymomials in x and use cubic */
  744.         /* interpolation in y for each of the 4 functions. */
  745.         xx = (x - xz) / dx;
  746.         yy = (y - yz) / dy;
  747.         fx1 = ((a[0][0] * xx + b[0][0]) * xx + c[0][0]) * xx + d[0][0];
  748.         fx2 = ((a[1][0] * xx + b[1][0]) * xx + c[1][0]) * xx + d[1][0];
  749.         fx3 = ((a[2][0] * xx + b[2][0]) * xx + c[2][0]) * xx + d[2][0];
  750.         fx4 = ((a[3][0] * xx + b[3][0]) * xx + c[3][0]) * xx + d[3][0];
  751.         p1 = fx4 - fx1 + 3. * (fx2 - fx3);
  752.         p2 = 3. * (fx1 - 2. * fx2 + fx3);
  753.         p3 = 6. * fx3 - 2. * fx1 - 3. * fx2 - fx4;
  754.         *f1 = ((p1 * yy + p2) * yy + p3) * yy * .1666666667 + fx2;
  755.         fx1 = ((a[0][1] * xx + b[0][1]) * xx + c[0][1]) * xx + d[0][1];
  756.         fx2 = ((a[1][1] * xx + b[1][1]) * xx + c[1][1]) * xx + d[1][1];
  757.         fx3 = ((a[2][1] * xx + b[2][1]) * xx + c[2][1]) * xx + d[2][1];
  758.         fx4 = ((a[3][1] * xx + b[3][1]) * xx + c[3][1]) * xx + d[3][1];
  759.         p1 = fx4 - fx1 + 3. * (fx2 - fx3);
  760.         p2 = 3. * (fx1 - 2. * fx2 + fx3);
  761.         p3 = 6. * fx3 - 2. * fx1 - 3. * fx2 - fx4;
  762.         *f2 = ((p1 * yy + p2) * yy + p3) * yy * .1666666667 + fx2;
  763.         fx1 = ((a[0][2] * xx + b[0][2]) * xx + c[0][2]) * xx + d[0][2];
  764.         fx2 = ((a[1][2] * xx + b[1][2]) * xx + c[1][2]) * xx + d[1][2];
  765.         fx3 = ((a[2][2] * xx + b[2][2]) * xx + c[2][2]) * xx + d[2][2];
  766.         fx4 = ((a[3][2] * xx + b[3][2]) * xx + c[3][2]) * xx + d[3][2];
  767.         p1 = fx4 - fx1 + 3. * (fx2 - fx3);
  768.         p2 = 3. * (fx1 - 2. * fx2 + fx3);
  769.         p3 = 6. * fx3 - 2. * fx1 - 3. * fx2 - fx4;
  770.         *f3 = ((p1 * yy + p2) * yy + p3) * yy * .1666666667 + fx2;
  771.         fx1 = ((a[0][3] * xx + b[0][3]) * xx + c[0][3]) * xx + d[0][3];
  772.         fx2 = ((a[1][3] * xx + b[1][3]) * xx + c[1][3]) * xx + d[1][3];
  773.         fx3 = ((a[2][3] * xx + b[2][3]) * xx + c[2][3]) * xx + d[2][3];
  774.         fx4 = ((a[3][3] * xx + b[3][3]) * xx + c[3][3]) * xx + d[3][3];
  775.         p1 = fx4 - fx1 + 3. * (fx2 - fx3);
  776.         p2 = 3. * (fx1 - 2. * fx2 + fx3);
  777.         p3 = 6. * fx3 - 2. * fx1 - 3. * fx2 - fx4;
  778.         *f4 = ((p1 * yy + p2) * yy + p3) * yy * .16666666670 + fx2;
  779.  
  780.         return;
  781. }
  782.  
  783. /*-----------------------------------------------------------------------*/
  784.  
  785. /* intx performs numerical integration of exp(jkr)/r by the method of */
  786. /* variable interval width romberg integration.  the integrand value */
  787. /* is supplied by subroutine gf. */
  788. void intx (double el1, double el2, double b, int ij, double *sgr, double *sgi)
  789. {
  790.         int ns, nt;
  791.         int nx = 1, nma = 65536, nts = 4;
  792.         int flag = TRUE;
  793.         double z, s, ze, fnm, ep, zend, fns, dz = 0., zp, dzot = 0., t00r, g1r, g5r, t00i;
  794.         double g1i, g5i, t01r, g3r, t01i, g3i, t10r, t10i, te1i, te1r, t02r;
  795.         double g2r, g4r, t02i, g2i, g4i, t11r, t11i, t20r, t20i, te2i, te2r;
  796.         double rx = 1.0e-4;
  797.  
  798.         z = el1;
  799.         ze = el2;
  800.         if (ij == 0)
  801.                 ze = 0.;
  802.         s = ze - z;
  803.         fnm = nma;
  804.         ep = s / (10. * fnm);
  805.         zend = ze - ep;
  806.         *sgr = 0.;
  807.         *sgi = 0.;
  808.         ns = nx;
  809.         nt = 0;
  810.         gf (z, &g1r, &g1i);
  811.  
  812.         while (TRUE)
  813.         {
  814.                 if (flag)
  815.                 {
  816.                         fns = ns;
  817.                         dz = s / fns;
  818.                         zp = z + dz;
  819.  
  820.                         if (zp > ze)
  821.                         {
  822.                                 dz = ze - z;
  823.                                 if (fabsl (dz) <= ep)
  824.                                 {
  825.                                         /* add contribution of near singularity for diagonal
  826.                                          * term */
  827.                                         if (ij == 0)
  828.                                         {
  829.                                                 *sgr =
  830.                                                     2. *
  831.                                                     (*sgr +
  832.                                                      logl ((sqrtl (b * b + s * s) + s) / b));
  833.                                                 *sgi = 2. * *sgi;
  834.                                         }
  835.                                         return;
  836.                                 }
  837.  
  838.                         } /* if( zp > ze) */
  839.  
  840.                         dzot = dz * .5;
  841.                         zp = z + dzot;
  842.                         gf (zp, &g3r, &g3i);
  843.                         zp = z + dz;
  844.                         gf (zp, &g5r, &g5i);
  845.  
  846.                 } /* if( flag ) */
  847.  
  848.                 t00r = (g1r + g5r) * dzot;
  849.                 t00i = (g1i + g5i) * dzot;
  850.                 t01r = (t00r + dz * g3r) * 0.5;
  851.                 t01i = (t00i + dz * g3i) * 0.5;
  852.                 t10r = (4.0 * t01r - t00r) / 3.0;
  853.                 t10i = (4.0 * t01i - t00i) / 3.0;
  854.  
  855.                 /* test convergence of 3 point romberg result. */
  856.                 test (t01r, t10r, &te1r, t01i, t10i, &te1i, 0.);
  857.                 if ((te1i <= rx) && (te1r <= rx))
  858.                 {
  859.                         *sgr = *sgr + t10r;
  860.                         *sgi = *sgi + t10i;
  861.                         nt += 2;
  862.  
  863.                         z += dz;
  864.                         if (z >= zend)
  865.                         {
  866.                                 /* add contribution of near singularity for diagonal term */
  867.                                 if (ij == 0)
  868.                                 {
  869.                                         *sgr = 2. *
  870.                                                (*sgr + logl ((sqrtl (b * b + s * s) + s) / b));
  871.                                         *sgi = 2. * *sgi;
  872.                                 }
  873.                                 return;
  874.                         }
  875.  
  876.                         g1r = g5r;
  877.                         g1i = g5i;
  878.                         if (nt >= nts)
  879.                                 if (ns > nx)
  880.                                 {
  881.                                         /* Double step size */
  882.                                         ns = ns / 2;
  883.                                         nt = 1;
  884.                                 }
  885.                         flag = TRUE;
  886.                         continue;
  887.  
  888.                 } /* if( (te1i <= rx) && (te1r <= rx) ) */
  889.  
  890.                 zp = z + dz * 0.25;
  891.                 gf (zp, &g2r, &g2i);
  892.                 zp = z + dz * 0.75;
  893.                 gf (zp, &g4r, &g4i);
  894.                 t02r = (t01r + dzot * (g2r + g4r)) * 0.5;
  895.                 t02i = (t01i + dzot * (g2i + g4i)) * 0.5;
  896.                 t11r = (4.0 * t02r - t01r) / 3.0;
  897.                 t11i = (4.0 * t02i - t01i) / 3.0;
  898.                 t20r = (16.0 * t11r - t10r) / 15.0;
  899.                 t20i = (16.0 * t11i - t10i) / 15.0;
  900.  
  901.                 /* test convergence of 5 point romberg result. */
  902.                 test (t11r, t20r, &te2r, t11i, t20i, &te2i, 0.);
  903.                 if ((te2i > rx) || (te2r > rx))
  904.                 {
  905.                         nt = 0;
  906.                         if (ns >= nma)
  907.                                 fprintf (output_fp, "\n  STEP SIZE LIMITED AT Z= %10.5G", z);
  908.                         else
  909.                         {
  910.                                 /* halve step size */
  911.                                 ns = ns * 2;
  912.                                 fns = ns;
  913.                                 dz = s / fns;
  914.                                 dzot = dz * 0.5;
  915.                                 g5r = g3r;
  916.                                 g5i = g3i;
  917.                                 g3r = g2r;
  918.                                 g3i = g2i;
  919.  
  920.                                 flag = FALSE;
  921.                                 continue;
  922.                         }
  923.  
  924.                 } /* if( (te2i > rx) || (te2r > rx) ) */
  925.  
  926.                 *sgr = *sgr + t20r;
  927.                 *sgi = *sgi + t20i;
  928.                 nt++;
  929.  
  930.                 z += dz;
  931.                 if (z >= zend)
  932.                 {
  933.                         /* add contribution of near singularity for diagonal term */
  934.                         if (ij == 0)
  935.                         {
  936.                                 *sgr = 2. * (*sgr + logl ((sqrtl (b * b + s * s) + s) / b));
  937.                                 *sgi = 2. * *sgi;
  938.                         }
  939.                         return;
  940.                 }
  941.  
  942.                 g1r = g5r;
  943.                 g1i = g5i;
  944.                 if (nt >= nts)
  945.                         if (ns > nx)
  946.                         {
  947.                                 /* Double step size */
  948.                                 ns = ns / 2;
  949.                                 nt = 1;
  950.                         }
  951.                 flag = TRUE;
  952.  
  953.         } /* while( TRUE ) */
  954. }
  955.  
  956. /*-----------------------------------------------------------------------*/
  957.  
  958. /* returns smallest of two arguments */
  959. int min (int a, int b)
  960. {
  961.         if (a < b)
  962.                 return (a);
  963.         else
  964.                 return (b);
  965. }
  966.  
  967. /*-----------------------------------------------------------------------*/
  968.  
  969. /* test for convergence in numerical integration */
  970. void test (double f1r, double f2r, double *tr, double f1i, double f2i, double *ti, double dmin)
  971. {
  972.         double den;
  973.  
  974.         den = fabsl (f2r);
  975.         *tr = fabsl (f2i);
  976.  
  977.         if (den < *tr)
  978.                 den = *tr;
  979.         if (den < dmin)
  980.                 den = dmin;
  981.  
  982.         if (den < 1.0e-37)
  983.         {
  984.                 *tr = 0.;
  985.                 *ti = 0.;
  986.                 return;
  987.         }
  988.  
  989.         *tr = fabsl ((f1r - f2r) / den);
  990.         *ti = fabsl ((f1i - f2i) / den);
  991.  
  992.         return;
  993. }
  994.  
  995. /*-----------------------------------------------------------------------*/
  996.  
  997. /* compute component of basis function i on segment is. */
  998. void sbf (int i, int is, double *aa, double *bb, double *cc)
  999. {
  1000.         int ix, jsno, june, jcox, jcoxx, jend, iend, njun1 = 0, njun2;
  1001.         double d, sig, pp, sdh, cdh, sd, omc, aj, pm = 0, cd, ap, qp, qm, xxi;
  1002.  
  1003.         *aa = 0.;
  1004.         *bb = 0.;
  1005.         *cc = 0.;
  1006.         june = 0;
  1007.         jsno = 0;
  1008.         pp = 0.;
  1009.         ix = i - 1;
  1010.  
  1011.         jcox = data.icon1[ix];
  1012.         if (jcox > PCHCON)
  1013.                 jcox = i;
  1014.         jcoxx = jcox - 1;
  1015.  
  1016.         jend = -1;
  1017.         iend = -1;
  1018.         sig = -1.;
  1019.  
  1020.         do
  1021.         {
  1022.                 if (jcox != 0)
  1023.                 {
  1024.                         if (jcox < 0)
  1025.                                 jcox = -jcox;
  1026.                         else
  1027.                         {
  1028.                                 sig = -sig;
  1029.                                 jend = -jend;
  1030.                         }
  1031.  
  1032.                         jcoxx = jcox - 1;
  1033.                         jsno++;
  1034.                         d = PI * data.si[jcoxx];
  1035.                         sdh = sinl (d);
  1036.                         cdh = cosl (d);
  1037.                         sd = 2. * sdh * cdh;
  1038.  
  1039.                         if (d <= 0.015)
  1040.                         {
  1041.                                 omc = 4. * d * d;
  1042.                                 omc =
  1043.                                     ((1.3888889e-3 * omc - 4.1666666667e-2) * omc + .5) * omc;
  1044.                         }
  1045.                         else
  1046.                                 omc = 1. - cdh * cdh + sdh * sdh;
  1047.  
  1048.                         aj = 1. / (logl (1. / (PI * data.bi[jcoxx])) - .577215664);
  1049.                         pp -= omc / sd * aj;
  1050.  
  1051.                         if (jcox == is)
  1052.                         {
  1053.                                 *aa = aj / sd * sig;
  1054.                                 *bb = aj / (2. * cdh);
  1055.                                 *cc = -aj / (2. * sdh) * sig;
  1056.                                 june = iend;
  1057.                         }
  1058.  
  1059.                         if (jcox != i)
  1060.                         {
  1061.                                 if (jend != 1)
  1062.                                         jcox = data.icon1[jcoxx];
  1063.                                 else
  1064.                                         jcox = data.icon2[jcoxx];
  1065.  
  1066.                                 if (abs (jcox) != i)
  1067.                                 {
  1068.                                         if (jcox == 0)
  1069.                                         {
  1070.                                                 fprintf (
  1071.                                                     output_fp,
  1072.                                                     "\n  SBF - SEGMENT CONNECTION ERROR FOR "
  1073.                                                     "SEGMENT %d",
  1074.                                                     i);
  1075.                                                 stop (-1);
  1076.                                         }
  1077.                                         else
  1078.                                                 continue;
  1079.                                 }
  1080.  
  1081.                         } /* if( jcox != i ) */
  1082.                         else if (jcox == is)
  1083.                                 *bb = -*bb;
  1084.  
  1085.                         if (iend == 1)
  1086.                                 break;
  1087.  
  1088.                 } /* if( jcox != 0 ) */
  1089.  
  1090.                 pm = -pp;
  1091.                 pp = 0.;
  1092.                 njun1 = jsno;
  1093.  
  1094.                 jcox = data.icon2[ix];
  1095.                 if (jcox > PCHCON)
  1096.                         jcox = i;
  1097.  
  1098.                 jend = 1;
  1099.                 iend = 1;
  1100.                 sig = -1.;
  1101.  
  1102.         } /* do */
  1103.         while (jcox != 0);
  1104.  
  1105.         njun2 = jsno - njun1;
  1106.         d = PI * data.si[ix];
  1107.         sdh = sinl (d);
  1108.         cdh = cosl (d);
  1109.         sd = 2. * sdh * cdh;
  1110.         cd = cdh * cdh - sdh * sdh;
  1111.  
  1112.         if (d <= 0.015)
  1113.         {
  1114.                 omc = 4. * d * d;
  1115.                 omc = ((1.3888889e-3 * omc - 4.1666666667e-2) * omc + .5) * omc;
  1116.         }
  1117.         else
  1118.                 omc = 1. - cd;
  1119.  
  1120.         ap = 1. / (logl (1. / (PI * data.bi[ix])) - .577215664);
  1121.         aj = ap;
  1122.  
  1123.         if (njun1 == 0)
  1124.         {
  1125.                 if (njun2 == 0)
  1126.                 {
  1127.                         *aa = -1.;
  1128.                         qp = PI * data.bi[ix];
  1129.                         xxi = qp * qp;
  1130.                         xxi = qp * (1. - .5 * xxi) / (1. - xxi);
  1131.                         *cc = 1. / (cdh - xxi * sdh);
  1132.                         return;
  1133.                 }
  1134.  
  1135.                 qp = PI * data.bi[ix];
  1136.                 xxi = qp * qp;
  1137.                 xxi = qp * (1. - .5 * xxi) / (1. - xxi);
  1138.                 qp = -(omc + xxi * sd) / (sd * (ap + xxi * pp) + cd * (xxi * ap - pp));
  1139.  
  1140.                 if (june == 1)
  1141.                 {
  1142.                         *aa = -*aa * qp;
  1143.                         *bb = *bb * qp;
  1144.                         *cc = -*cc * qp;
  1145.                         if (i != is)
  1146.                                 return;
  1147.                 }
  1148.  
  1149.                 *aa -= 1.;
  1150.                 d = cd - xxi * sd;
  1151.                 *bb += (sdh + ap * qp * (cdh - xxi * sdh)) / d;
  1152.                 *cc += (cdh + ap * qp * (sdh + xxi * cdh)) / d;
  1153.                 return;
  1154.  
  1155.         } /* if( njun1 == 0) */
  1156.  
  1157.         if (njun2 == 0)
  1158.         {
  1159.                 qm = PI * data.bi[ix];
  1160.                 xxi = qm * qm;
  1161.                 xxi = qm * (1. - .5 * xxi) / (1. - xxi);
  1162.                 qm = (omc + xxi * sd) / (sd * (aj - xxi * pm) + cd * (pm + xxi * aj));
  1163.  
  1164.                 if (june == -1)
  1165.                 {
  1166.                         *aa = *aa * qm;
  1167.                         *bb = *bb * qm;
  1168.                         *cc = *cc * qm;
  1169.                         if (i != is)
  1170.                                 return;
  1171.                 }
  1172.  
  1173.                 *aa -= 1.;
  1174.                 d = cd - xxi * sd;
  1175.                 *bb += (aj * qm * (cdh - xxi * sdh) - sdh) / d;
  1176.                 *cc += (cdh - aj * qm * (sdh + xxi * cdh)) / d;
  1177.                 return;
  1178.  
  1179.         } /* if( njun2 == 0) */
  1180.  
  1181.         qp = sd * (pm * pp + aj * ap) + cd * (pm * ap - pp * aj);
  1182.         qm = (ap * omc - pp * sd) / qp;
  1183.         qp = -(aj * omc + pm * sd) / qp;
  1184.  
  1185.         if (june != 0)
  1186.         {
  1187.                 if (june < 0)
  1188.                 {
  1189.                         *aa = *aa * qm;
  1190.                         *bb = *bb * qm;
  1191.                         *cc = *cc * qm;
  1192.                 }
  1193.                 else
  1194.                 {
  1195.                         *aa = -*aa * qp;
  1196.                         *bb = *bb * qp;
  1197.                         *cc = -*cc * qp;
  1198.                 }
  1199.  
  1200.                 if (i != is)
  1201.                         return;
  1202.  
  1203.         } /* if( june != 0 ) */
  1204.  
  1205.         *aa -= 1.;
  1206.         *bb += (aj * qm + ap * qp) * sdh / sd;
  1207.         *cc += (aj * qm - ap * qp) * cdh / sd;
  1208.  
  1209.         return;
  1210. }
  1211.  
  1212. /*-----------------------------------------------------------------------*/
  1213.  
  1214. /* compute basis function i */
  1215. void tbf (int i, int icap)
  1216. {
  1217.         int ix, jcox, jcoxx, jend, iend, njun1 = 0, njun2, jsnop, jsnox;
  1218.         double pp, sdh, cdh, sd, omc, aj, pm = 0, cd, ap, qp, qm, xxi;
  1219.         double d, sig; /*** also global ***/
  1220.  
  1221.         segj.jsno = 0;
  1222.         pp = 0.;
  1223.         ix = i - 1;
  1224.         jcox = data.icon1[ix];
  1225.  
  1226.         if (jcox > PCHCON)
  1227.                 jcox = i;
  1228.  
  1229.         jend = -1;
  1230.         iend = -1;
  1231.         sig = -1.;
  1232.  
  1233.         do
  1234.         {
  1235.                 if (jcox != 0)
  1236.                 {
  1237.                         if (jcox < 0)
  1238.                                 jcox = -jcox;
  1239.                         else
  1240.                         {
  1241.                                 sig = -sig;
  1242.                                 jend = -jend;
  1243.                         }
  1244.  
  1245.                         jcoxx = jcox - 1;
  1246.                         segj.jsno++;
  1247.                         jsnox = segj.jsno - 1;
  1248.                         segj.jco[jsnox] = jcox;
  1249.                         d = PI * data.si[jcoxx];
  1250.                         sdh = sinl (d);
  1251.                         cdh = cosl (d);
  1252.                         sd = 2. * sdh * cdh;
  1253.  
  1254.                         if (d <= 0.015)
  1255.                         {
  1256.                                 omc = 4. * d * d;
  1257.                                 omc =
  1258.                                     ((1.3888889e-3 * omc - 4.1666666667e-2) * omc + .5) * omc;
  1259.                         }
  1260.                         else
  1261.                                 omc = 1. - cdh * cdh + sdh * sdh;
  1262.  
  1263.                         aj = 1. / (logl (1. / (PI * data.bi[jcoxx])) - .577215664);
  1264.                         pp = pp - omc / sd * aj;
  1265.                         segj.ax[jsnox] = aj / sd * sig;
  1266.                         segj.bx[jsnox] = aj / (2. * cdh);
  1267.                         segj.cx[jsnox] = -aj / (2. * sdh) * sig;
  1268.  
  1269.                         if (jcox != i)
  1270.                         {
  1271.                                 if (jend == 1)
  1272.                                         jcox = data.icon2[jcoxx];
  1273.                                 else
  1274.                                         jcox = data.icon1[jcoxx];
  1275.  
  1276.                                 if (abs (jcox) != i)
  1277.                                 {
  1278.                                         if (jcox != 0)
  1279.                                                 continue;
  1280.                                         else
  1281.                                         {
  1282.                                                 fprintf (
  1283.                                                     output_fp,
  1284.                                                     "\n  TBF - SEGMENT CONNECTION ERROR FOR "
  1285.                                                     "SEGMENT %5d",
  1286.                                                     i);
  1287.                                                 stop (-1);
  1288.                                         }
  1289.                                 }
  1290.  
  1291.                         } /* if( jcox != i) */
  1292.                         else
  1293.                                 segj.bx[jsnox] = -segj.bx[jsnox];
  1294.  
  1295.                         if (iend == 1)
  1296.                                 break;
  1297.  
  1298.                 } /* if( jcox != 0 ) */
  1299.  
  1300.                 pm = -pp;
  1301.                 pp = 0.;
  1302.                 njun1 = segj.jsno;
  1303.  
  1304.                 jcox = data.icon2[ix];
  1305.                 if (jcox > PCHCON)
  1306.                         jcox = i;
  1307.  
  1308.                 jend = 1;
  1309.                 iend = 1;
  1310.                 sig = -1.;
  1311.  
  1312.         } /* do */
  1313.         while (jcox != 0);
  1314.  
  1315.         njun2 = segj.jsno - njun1;
  1316.         jsnop = segj.jsno;
  1317.         segj.jco[jsnop] = i;
  1318.         d = PI * data.si[ix];
  1319.         sdh = sinl (d);
  1320.         cdh = cosl (d);
  1321.         sd = 2. * sdh * cdh;
  1322.         cd = cdh * cdh - sdh * sdh;
  1323.  
  1324.         if (d <= 0.015)
  1325.         {
  1326.                 omc = 4. * d * d;
  1327.                 omc = ((1.3888889e-3 * omc - 4.1666666667e-2) * omc + .5) * omc;
  1328.         }
  1329.         else
  1330.                 omc = 1. - cd;
  1331.  
  1332.         ap = 1. / (logl (1. / (PI * data.bi[ix])) - .577215664);
  1333.         aj = ap;
  1334.  
  1335.         if (njun1 == 0)
  1336.         {
  1337.                 if (njun2 == 0)
  1338.                 {
  1339.                         segj.bx[jsnop] = 0.;
  1340.  
  1341.                         if (icap == 0)
  1342.                                 xxi = 0.;
  1343.                         else
  1344.                         {
  1345.                                 qp = PI * data.bi[ix];
  1346.                                 xxi = qp * qp;
  1347.                                 xxi = qp * (1. - .5 * xxi) / (1. - xxi);
  1348.                         }
  1349.  
  1350.                         segj.cx[jsnop] = 1. / (cdh - xxi * sdh);
  1351.                         segj.jsno = jsnop + 1;
  1352.                         segj.ax[jsnop] = -1.;
  1353.                         return;
  1354.  
  1355.                 } /* if( njun2 == 0) */
  1356.  
  1357.                 if (icap == 0)
  1358.                         xxi = 0.;
  1359.                 else
  1360.                 {
  1361.                         qp = PI * data.bi[ix];
  1362.                         xxi = qp * qp;
  1363.                         xxi = qp * (1. - .5 * xxi) / (1. - xxi);
  1364.                 }
  1365.  
  1366.                 qp = -(omc + xxi * sd) / (sd * (ap + xxi * pp) + cd * (xxi * ap - pp));
  1367.                 d = cd - xxi * sd;
  1368.                 segj.bx[jsnop] = (sdh + ap * qp * (cdh - xxi * sdh)) / d;
  1369.                 segj.cx[jsnop] = (cdh + ap * qp * (sdh + xxi * cdh)) / d;
  1370.  
  1371.                 for (iend = 0; iend < njun2; iend++)
  1372.                 {
  1373.                         segj.ax[iend] = -segj.ax[iend] * qp;
  1374.                         segj.bx[iend] = segj.bx[iend] * qp;
  1375.                         segj.cx[iend] = -segj.cx[iend] * qp;
  1376.                 }
  1377.  
  1378.                 segj.jsno = jsnop + 1;
  1379.                 segj.ax[jsnop] = -1.;
  1380.                 return;
  1381.  
  1382.         } /* if( njun1 == 0) */
  1383.  
  1384.         if (njun2 == 0)
  1385.         {
  1386.                 if (icap == 0)
  1387.                         xxi = 0.;
  1388.                 else
  1389.                 {
  1390.                         qm = PI * data.bi[ix];
  1391.                         xxi = qm * qm;
  1392.                         xxi = qm * (1. - .5 * xxi) / (1. - xxi);
  1393.                 }
  1394.  
  1395.                 qm = (omc + xxi * sd) / (sd * (aj - xxi * pm) + cd * (pm + xxi * aj));
  1396.                 d = cd - xxi * sd;
  1397.                 segj.bx[jsnop] = (aj * qm * (cdh - xxi * sdh) - sdh) / d;
  1398.                 segj.cx[jsnop] = (cdh - aj * qm * (sdh + xxi * cdh)) / d;
  1399.  
  1400.                 for (iend = 0; iend < njun1; iend++)
  1401.                 {
  1402.                         segj.ax[iend] = segj.ax[iend] * qm;
  1403.                         segj.bx[iend] = segj.bx[iend] * qm;
  1404.                         segj.cx[iend] = segj.cx[iend] * qm;
  1405.                 }
  1406.  
  1407.                 segj.jsno = jsnop + 1;
  1408.                 segj.ax[jsnop] = -1.;
  1409.                 return;
  1410.  
  1411.         } /* if( njun2 == 0) */
  1412.  
  1413.         qp = sd * (pm * pp + aj * ap) + cd * (pm * ap - pp * aj);
  1414.         qm = (ap * omc - pp * sd) / qp;
  1415.         qp = -(aj * omc + pm * sd) / qp;
  1416.         segj.bx[jsnop] = (aj * qm + ap * qp) * sdh / sd;
  1417.         segj.cx[jsnop] = (aj * qm - ap * qp) * cdh / sd;
  1418.  
  1419.         for (iend = 0; iend < njun1; iend++)
  1420.         {
  1421.                 segj.ax[iend] = segj.ax[iend] * qm;
  1422.                 segj.bx[iend] = segj.bx[iend] * qm;
  1423.                 segj.cx[iend] = segj.cx[iend] * qm;
  1424.         }
  1425.  
  1426.         jend = njun1;
  1427.         for (iend = jend; iend < segj.jsno; iend++)
  1428.         {
  1429.                 segj.ax[iend] = -segj.ax[iend] * qp;
  1430.                 segj.bx[iend] = segj.bx[iend] * qp;
  1431.                 segj.cx[iend] = -segj.cx[iend] * qp;
  1432.         }
  1433.  
  1434.         segj.jsno = jsnop + 1;
  1435.         segj.ax[jsnop] = -1.;
  1436. }
  1437.  
  1438. /*-----------------------------------------------------------------------*/
  1439.  
  1440. /* compute the components of all basis functions on segment j */
  1441. void trio (int j)
  1442. {
  1443.         int jcox, jcoxx, jsnox, jx, jend = 0, iend = 0;
  1444.  
  1445.         segj.jsno = 0;
  1446.         jx = j - 1;
  1447.         jcox = data.icon1[jx];
  1448.         jcoxx = jcox - 1;
  1449.  
  1450.         if (jcox <= PCHCON)
  1451.         {
  1452.                 jend = -1;
  1453.                 iend = -1;
  1454.         }
  1455.  
  1456.         if ((jcox == 0) || (jcox > PCHCON))
  1457.         {
  1458.                 jcox = data.icon2[jx];
  1459.                 jcoxx = jcox - 1;
  1460.  
  1461.                 if (jcox <= PCHCON)
  1462.                 {
  1463.                         jend = 1;
  1464.                         iend = 1;
  1465.                 }
  1466.  
  1467.                 if (jcox == 0 || (jcox > PCHCON))
  1468.                 {
  1469.                         jsnox = segj.jsno;
  1470.                         segj.jsno++;
  1471.  
  1472.                         /* Allocate to connections buffers */
  1473.                         if (segj.jsno >= segj.maxcon)
  1474.                         {
  1475.                                 segj.maxcon = segj.jsno + 1;
  1476.                                 mem_realloc ((void *) &segj.jco, segj.maxcon * sizeof (int));
  1477.                                 mem_realloc ((void *) &segj.ax, segj.maxcon * sizeof (double));
  1478.                                 mem_realloc ((void *) &segj.bx, segj.maxcon * sizeof (double));
  1479.                                 mem_realloc ((void *) &segj.cx, segj.maxcon * sizeof (double));
  1480.                         }
  1481.  
  1482.                         sbf (j, j, &segj.ax[jsnox], &segj.bx[jsnox], &segj.cx[jsnox]);
  1483.                         segj.jco[jsnox] = j;
  1484.                         return;
  1485.                 }
  1486.  
  1487.         } /* if( (jcox == 0) || (jcox > PCHCON) ) */
  1488.  
  1489.         do
  1490.         {
  1491.                 if (jcox < 0)
  1492.                         jcox = -jcox;
  1493.                 else
  1494.                         jend = -jend;
  1495.                 jcoxx = jcox - 1;
  1496.  
  1497.                 if (jcox != j)
  1498.                 {
  1499.                         jsnox = segj.jsno;
  1500.                         segj.jsno++;
  1501.  
  1502.                         /* Allocate to connections buffers */
  1503.                         if (segj.jsno >= segj.maxcon)
  1504.                         {
  1505.                                 segj.maxcon = segj.jsno + 1;
  1506.                                 mem_realloc ((void *) &segj.jco, segj.maxcon * sizeof (int));
  1507.                                 mem_realloc ((void *) &segj.ax, segj.maxcon * sizeof (double));
  1508.                                 mem_realloc ((void *) &segj.bx, segj.maxcon * sizeof (double));
  1509.                                 mem_realloc ((void *) &segj.cx, segj.maxcon * sizeof (double));
  1510.                         }
  1511.  
  1512.                         sbf (jcox, j, &segj.ax[jsnox], &segj.bx[jsnox], &segj.cx[jsnox]);
  1513.                         segj.jco[jsnox] = jcox;
  1514.  
  1515.                         if (jend != 1)
  1516.                                 jcox = data.icon1[jcoxx];
  1517.                         else
  1518.                                 jcox = data.icon2[jcoxx];
  1519.  
  1520.                         if (jcox == 0)
  1521.                         {
  1522.                                 fprintf (
  1523.                                     output_fp,
  1524.                                     "\n  TRIO - SEGMENT CONNENTION ERROR FOR SEGMENT %5d",
  1525.                                     j);
  1526.                                 stop (-1);
  1527.                         }
  1528.                         else
  1529.                                 continue;
  1530.  
  1531.                 } /* if( jcox != j) */
  1532.  
  1533.                 if (iend == 1)
  1534.                         break;
  1535.  
  1536.                 jcox = data.icon2[jx];
  1537.  
  1538.                 if (jcox > PCHCON)
  1539.                         break;
  1540.  
  1541.                 jend = 1;
  1542.                 iend = 1;
  1543.  
  1544.         } /* do */
  1545.         while (jcox != 0);
  1546.  
  1547.         jsnox = segj.jsno;
  1548.         segj.jsno++;
  1549.  
  1550.         /* Allocate to connections buffers */
  1551.         if (segj.jsno >= segj.maxcon)
  1552.         {
  1553.                 segj.maxcon = segj.jsno + 1;
  1554.                 mem_realloc ((void *) &segj.jco, segj.maxcon * sizeof (int));
  1555.                 mem_realloc ((void *) &segj.ax, segj.maxcon * sizeof (double));
  1556.                 mem_realloc ((void *) &segj.bx, segj.maxcon * sizeof (double));
  1557.                 mem_realloc ((void *) &segj.cx, segj.maxcon * sizeof (double));
  1558.         }
  1559.  
  1560.         sbf (j, j, &segj.ax[jsnox], &segj.bx[jsnox], &segj.cx[jsnox]);
  1561.         segj.jco[jsnox] = j;
  1562.  
  1563.         return;
  1564. }
  1565.  
  1566. /*-----------------------------------------------------------------------*/
  1567.  
  1568. /* zint computes the internal impedance of a circular wire */
  1569. void zint (double sigl, double rolam, complex double *zint)
  1570. {
  1571. #define cc1 (6.0e-7 + 1.9e-6fj)
  1572. #define cc2 (-3.4e-6 + 5.1e-6fj)
  1573. #define cc3 (-2.52e-5 + 0.fj)
  1574. #define cc4 (-9.06e-5 - 9.01e-5fj)
  1575. #define cc5 (0. - 9.765e-4fj)
  1576. #define cc6 (.0110486 - .0110485fj)
  1577. #define cc7 (0. - .3926991fj)
  1578. #define cc8 (1.6e-6 - 3.2e-6fj)
  1579. #define cc9 (1.17e-5 - 2.4e-6fj)
  1580. #define cc10 (3.46e-5 + 3.38e-5fj)
  1581. #define cc11 (5.0e-7 + 2.452e-4fj)
  1582. #define cc12 (-1.3813e-3 + 1.3811e-3fj)
  1583. #define cc13 (-6.25001e-2 - 1.0e-7fj)
  1584. #define cc14 (.7071068 + .7071068fj)
  1585. #define cn cc14
  1586.  
  1587. #define th(d)                                                                                 \
  1588.         ((((((cc1 * (d) + cc2) * (d) + cc3) * (d) + cc4) * (d) + cc5) * (d) + cc6) * (d) + cc7)
  1589. #define ph(d)                                                                                 \
  1590.         ((((((cc8 * (d) + cc9) * (d) + cc10) * (d) + cc11) * (d) + cc12) * (d) + cc13) *      \
  1591.              (d) +                                                                            \
  1592.          cc14)
  1593. #define f(d) (csqrtl (POT / (d)) * cexpl (-cn * (d) + th (-8. / x)))
  1594. #define g(d) (cexpl (cn * (d) + th (8. / x)) / csqrtl (TP * (d)))
  1595.  
  1596.         double x, y, s, ber, bei;
  1597.         double tpcmu = 2.368705e+3;
  1598.         ;
  1599.         double cmotp = 60.00;
  1600.         complex double br1, br2;
  1601.  
  1602.         x = sqrtl (tpcmu * sigl) * rolam;
  1603.         if (x <= 110.)
  1604.         {
  1605.                 if (x <= 8.)
  1606.                 {
  1607.                         y = x / 8.;
  1608.                         y = y * y;
  1609.                         s = y * y;
  1610.  
  1611.                         ber = ((((((-9.01e-6 * s + 1.22552e-3) * s - .08349609) * s +
  1612.                                   2.6419140) *
  1613.                                      s -
  1614.                                  32.363456) *
  1615.                                     s +
  1616.                                 113.77778) *
  1617.                                    s -
  1618.                                64.) *
  1619.                                   s +
  1620.                               1.;
  1621.  
  1622.                         bei = ((((((1.1346e-4 * s - .01103667) * s + .52185615) * s -
  1623.                                   10.567658) *
  1624.                                      s +
  1625.                                  72.817777) *
  1626.                                     s -
  1627.                                 113.77778) *
  1628.                                    s +
  1629.                                16.) *
  1630.                               y;
  1631.  
  1632.                         br1 = cmplx (ber, bei);
  1633.  
  1634.                         ber = (((((((-3.94e-6 * s + 4.5957e-4) * s - .02609253) * s +
  1635.                                    .66047849) *
  1636.                                       s -
  1637.                                   6.0681481) *
  1638.                                      s +
  1639.                                  14.222222) *
  1640.                                     s -
  1641.                                 4.) *
  1642.                                y) *
  1643.                               x;
  1644.  
  1645.                         bei = ((((((4.609e-5 * s - 3.79386e-3) * s + .14677204) * s -
  1646.                                   2.3116751) *
  1647.                                      s +
  1648.                                  11.377778) *
  1649.                                     s -
  1650.                                 10.666667) *
  1651.                                    s +
  1652.                                .5) *
  1653.                               x;
  1654.  
  1655.                         br2 = cmplx (ber, bei);
  1656.                         br1 = br1 / br2;
  1657.                         *zint = CPLX_01 * sqrtl (cmotp / sigl) * br1 / rolam;
  1658.  
  1659.                 } /* if( x <= 8.) */
  1660.  
  1661.                 br2 = CPLX_01 * f (x) / PI;
  1662.                 br1 = g (x) + br2;
  1663.                 br2 = g (x) * ph (8. / x) - br2 * ph (-8. / x);
  1664.                 br1 = br1 / br2;
  1665.                 *zint = CPLX_01 * sqrtl (cmotp / sigl) * br1 / rolam;
  1666.  
  1667.         } /* if( x <= 110.) */
  1668.  
  1669.         br1 = cmplx (.70710678, -.70710678);
  1670.         *zint = CPLX_01 * sqrtl (cmotp / sigl) * br1 / rolam;
  1671. }
  1672.  
  1673. /*-----------------------------------------------------------------------*/
  1674.  
  1675. /* cang returns the phase angle of a complex number in degrees. */
  1676. double cang (complex double z)
  1677. {
  1678.         return (cargl (z) * TD);
  1679. }
  1680.  
  1681. /*-----------------------------------------------------------------------*/
  1682.