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. #include <omp.h>
  29.  
  30. /* common  /data/ */
  31. extern data_t data;
  32.  
  33. /* common  /dataj/ */
  34. extern dataj_t dataj;
  35.  
  36. /* common  /matpar/ */
  37. extern matpar_t matpar;
  38.  
  39. /* common  /segj/ */
  40. extern segj_t segj;
  41.  
  42. /* common  /zload/ */
  43. extern zload_t zload;
  44.  
  45. /* common  /smat/ */
  46. extern smat_t smat;
  47.  
  48. /* common  /gnd/ */
  49. extern gnd_t gnd;
  50.  
  51. /* common  /vsorc/ */
  52. extern vsorc_t vsorc;
  53.  
  54. /* pointers to input/output files */
  55. extern FILE *input_fp, *output_fp, *plot_fp;
  56.  
  57. /*-------------------------------------------------------------------*/
  58.  
  59. /* cmset sets up the complex structure matrix in the array cm */
  60. void cmset (int nrow, complex double *cm, double rkhx, int iexkx)
  61. {
  62.         int mp2, neq, npeq, it, i, j, i1, i2, in2;
  63.         // int iout;
  64.         int im1, im2, ist, ij, ipr, jss, jm1, jm2, jst, k, ka, kk;
  65.         complex double zaj, deter, *scm = NULL;
  66.  
  67.         mp2 = 2 * data.mp;
  68.         npeq = data.np + mp2;
  69.         neq = data.n + 2 * data.m;
  70.         smat.nop = neq / npeq;
  71.  
  72.         dataj.rkh = rkhx;
  73.         dataj.iexk = iexkx;
  74.         // iout=2* matpar.npblk* nrow;
  75.         it = matpar.nlast;
  76.  
  77.         for (i = 0; i < nrow; i++)
  78.                 for (j = 0; j < it; j++)
  79.                         cm[i + j * nrow] = CPLX_00;
  80.  
  81.         i1 = 1;
  82.         i2 = it;
  83.         in2 = i2;
  84.  
  85.         if (in2 > data.np)
  86.                 in2 = data.np;
  87.  
  88.         im1 = i1 - data.np;
  89.         im2 = i2 - data.np;
  90.  
  91.         if (im1 < 1)
  92.                 im1 = 1;
  93.  
  94.         ist = 1;
  95.         if (i1 <= data.np)
  96.                 ist = data.np - i1 + 2;
  97.  
  98.         /* wire source loop */
  99.         if (data.n != 0)
  100.         {
  101.                 for (j = 1; j <= data.n; j++)
  102.                 {
  103.                         trio (j);
  104.                         for (i = 0; i < segj.jsno; i++)
  105.                         {
  106.                                 ij = segj.jco[i];
  107.                                 segj.jco[i] = ((ij - 1) / data.np) * mp2 + ij;
  108.                         }
  109.  
  110.                         if (i1 <= in2)
  111.                                 cmww (j, i1, in2, cm, nrow, cm, nrow, 1);
  112.  
  113.                         if (im1 <= im2)
  114.                                 cmws (j, im1, im2, &cm[(ist - 1) * nrow], nrow, cm, nrow, 1);
  115.  
  116.                         /* matrix elements modified by loading */
  117.                         if (zload.nload == 0)
  118.                                 continue;
  119.  
  120.                         if (j > data.np)
  121.                                 continue;
  122.  
  123.                         ipr = j;
  124.                         if ((ipr < 1) || (ipr > it))
  125.                                 continue;
  126.  
  127.                         zaj = zload.zarray[j - 1];
  128.  
  129.                         for (i = 0; i < segj.jsno; i++)
  130.                         {
  131.                                 jss = segj.jco[i];
  132.                                 cm[(jss - 1) + (ipr - 1) * nrow] -=
  133.                                     (segj.ax[i] + segj.cx[i]) * zaj;
  134.                         }
  135.  
  136.                 } /* for( j = 1; j <= n; j++ ) */
  137.  
  138.         } /* if( n != 0) */
  139.  
  140.         if (data.m != 0)
  141.         {
  142.                 /* matrix elements for patch current sources */
  143.                 jm1 = 1 - data.mp;
  144.                 jm2 = 0;
  145.                 jst = 1 - mp2;
  146.  
  147.                 for (i = 0; i < smat.nop; i++)
  148.                 {
  149.                         jm1 += data.mp;
  150.                         jm2 += data.mp;
  151.                         jst += npeq;
  152.  
  153.                         if (i1 <= in2)
  154.                                 cmsw (jm1, jm2, i1, in2, &cm[(jst - 1)], cm, 0, nrow, 1);
  155.  
  156.                         if (im1 <= im2)
  157.                                 cmss (
  158.                                     jm1,
  159.                                     jm2,
  160.                                     im1,
  161.                                     im2,
  162.                                     &cm[(jst - 1) + (ist - 1) * nrow],
  163.                                     nrow,
  164.                                     1);
  165.                 }
  166.  
  167.         } /* if( m != 0) */
  168.  
  169.         if (matpar.icase == 1)
  170.                 return;
  171.  
  172.         /* Allocate to scratch memory */
  173.         mem_alloc ((void *) &scm, data.np2m * sizeof (complex double));
  174.  
  175.         /* combine elements for symmetry modes */
  176.         for (i = 0; i < it; i++)
  177.         {
  178.                 for (j = 0; j < npeq; j++)
  179.                 {
  180.                         for (k = 0; k < smat.nop; k++)
  181.                         {
  182.                                 ka = j + k * npeq;
  183.                                 scm[k] = cm[ka + i * nrow];
  184.                         }
  185.  
  186.                         deter = scm[0];
  187.  
  188.                         for (kk = 1; kk < smat.nop; kk++)
  189.                                 deter += scm[kk];
  190.  
  191.                         cm[j + i * nrow] = deter;
  192.  
  193.                         for (k = 1; k < smat.nop; k++)
  194.                         {
  195.                                 ka = j + k * npeq;
  196.                                 deter = scm[0];
  197.  
  198.                                 for (kk = 1; kk < smat.nop; kk++)
  199.                                 {
  200.                                         deter += scm[kk] * smat.ssx[k + kk * smat.nop];
  201.                                         cm[ka + i * nrow] = deter;
  202.                                 }
  203.  
  204.                         } /* for( k = 1; k < smat.nop; k++ ) */
  205.  
  206.                 } /* for( j = 0; j < npeq; j++ ) */
  207.  
  208.         } /* for( i = 0; i < it; i++ ) */
  209.  
  210.         free_ptr ((void *) &scm);
  211.  
  212.         return;
  213. }
  214.  
  215. /*-----------------------------------------------------------------------*/
  216.  
  217. /* cmss computes matrix elements for surface-surface interactions. */
  218. void cmss (int j1, int j2, int im1, int im2, complex double *cm, int nrow, int itrp)
  219. {
  220.         int i1, i2, icomp, ii1, i, il, ii2, jj1, j, jl, jj2;
  221.         double t1xi, t1yi, t1zi, t2xi, t2yi, t2zi, xi, yi, zi;
  222.         complex double g11, g12, g21, g22;
  223.  
  224.         i1 = (im1 + 1) / 2;
  225.         i2 = (im2 + 1) / 2;
  226.         icomp = i1 * 2 - 3;
  227.         ii1 = -2;
  228.         if (icomp + 2 < im1)
  229.                 ii1 = -3;
  230.  
  231.         /* loop over observation patches */
  232.         il = -1;
  233.         for (i = i1; i <= i2; i++)
  234.         {
  235.                 il++;
  236.                 icomp += 2;
  237.                 ii1 += 2;
  238.                 ii2 = ii1 + 1;
  239.  
  240.                 t1xi = data.t1x[il] * data.psalp[il];
  241.                 t1yi = data.t1y[il] * data.psalp[il];
  242.                 t1zi = data.t1z[il] * data.psalp[il];
  243.                 t2xi = data.t2x[il] * data.psalp[il];
  244.                 t2yi = data.t2y[il] * data.psalp[il];
  245.                 t2zi = data.t2z[il] * data.psalp[il];
  246.                 xi = data.px[il];
  247.                 yi = data.py[il];
  248.                 zi = data.pz[il];
  249.  
  250.                 /* loop over source patches */
  251.                 jj1 = -2;
  252.                 for (j = j1; j <= j2; j++)
  253.                 {
  254.                         jl = j - 1;
  255.                         jj1 += 2;
  256.                         jj2 = jj1 + 1;
  257.  
  258.                         dataj.s = data.pbi[jl];
  259.                         dataj.xj = data.px[jl];
  260.                         dataj.yj = data.py[jl];
  261.                         dataj.zj = data.pz[jl];
  262.                         dataj.t1xj = data.t1x[jl];
  263.                         dataj.t1yj = data.t1y[jl];
  264.                         dataj.t1zj = data.t1z[jl];
  265.                         dataj.t2xj = data.t2x[jl];
  266.                         dataj.t2yj = data.t2y[jl];
  267.                         dataj.t2zj = data.t2z[jl];
  268.  
  269.                         hintg (xi, yi, zi);
  270.  
  271.                         g11 = -(t2xi * dataj.exk + t2yi * dataj.eyk + t2zi * dataj.ezk);
  272.                         g12 = -(t2xi * dataj.exs + t2yi * dataj.eys + t2zi * dataj.ezs);
  273.                         g21 = -(t1xi * dataj.exk + t1yi * dataj.eyk + t1zi * dataj.ezk);
  274.                         g22 = -(t1xi * dataj.exs + t1yi * dataj.eys + t1zi * dataj.ezs);
  275.  
  276.                         if (i == j)
  277.                         {
  278.                                 g11 -= .5;
  279.                                 g22 += .5;
  280.                         }
  281.  
  282.                         /* normal fill */
  283.                         if (itrp == 0)
  284.                         {
  285.                                 if (icomp >= im1)
  286.                                 {
  287.                                         cm[ii1 + jj1 * nrow] = g11;
  288.                                         cm[ii1 + jj2 * nrow] = g12;
  289.                                 }
  290.  
  291.                                 if (icomp >= im2)
  292.                                         continue;
  293.  
  294.                                 cm[ii2 + jj1 * nrow] = g21;
  295.                                 cm[ii2 + jj2 * nrow] = g22;
  296.                                 continue;
  297.  
  298.                         } /* if( itrp == 0) */
  299.  
  300.                         /* transposed fill */
  301.                         if (icomp >= im1)
  302.                         {
  303.                                 cm[jj1 + ii1 * nrow] = g11;
  304.                                 cm[jj2 + ii1 * nrow] = g12;
  305.                         }
  306.  
  307.                         if (icomp >= im2)
  308.                                 continue;
  309.  
  310.                         cm[jj1 + ii2 * nrow] = g21;
  311.                         cm[jj2 + ii2 * nrow] = g22;
  312.  
  313.                 } /* for( j = j1; j <= j2; j++ ) */
  314.  
  315.         } /* for( i = i1; i <= i2; i++ ) */
  316.  
  317.         return;
  318. }
  319.  
  320. /*-----------------------------------------------------------------------*/
  321.  
  322. /* computes matrix elements for e along wires due to patch current */
  323. void cmsw (
  324.     int j1,
  325.     int j2,
  326.     int i1,
  327.     int i2,
  328.     complex double *cm,
  329.     complex double *cw,
  330.     int ncw,
  331.     int nrow,
  332.     int itrp)
  333. {
  334.         // int neqs;
  335.         int k, icgo, i, ipch, jl, j, js, il, ip;
  336.         int jsnox; /* -1 offset to "jsno" for array indexing */
  337.         double xi, yi, zi, cabi, sabi, salpi, fsign = 1., pyl, pxl;
  338.         complex double emel[9];
  339.  
  340.         // neqs= data.np2m;
  341.         jsnox = segj.jsno - 1;
  342.  
  343.         if (itrp >= 0)
  344.         {
  345.                 k = -1;
  346.                 icgo = 0;
  347.  
  348.                 /* observation loop */
  349.                 for (i = i1 - 1; i < i2; i++)
  350.                 {
  351.                         k++;
  352.                         xi = data.x[i];
  353.                         yi = data.y[i];
  354.                         zi = data.z[i];
  355.                         cabi = data.cab[i];
  356.                         sabi = data.sab[i];
  357.                         salpi = data.salp[i];
  358.                         ipch = 0;
  359.  
  360.                         if (data.icon1[i] >= PCHCON)
  361.                         {
  362.                                 ipch = data.icon1[i] - PCHCON;
  363.                                 fsign = -1.;
  364.                         }
  365.  
  366.                         if (data.icon2[i] >= PCHCON)
  367.                         {
  368.                                 ipch = data.icon2[i] - PCHCON;
  369.                                 fsign = 1.;
  370.                         }
  371.  
  372.                         /* source loop */
  373.                         jl = -1;
  374.                         for (j = j1; j <= j2; j++)
  375.                         {
  376.                                 jl += 2;
  377.                                 js = j - 1;
  378.                                 dataj.t1xj = data.t1x[js];
  379.                                 dataj.t1yj = data.t1y[js];
  380.                                 dataj.t1zj = data.t1z[js];
  381.                                 dataj.t2xj = data.t2x[js];
  382.                                 dataj.t2yj = data.t2y[js];
  383.                                 dataj.t2zj = data.t2z[js];
  384.                                 dataj.xj = data.px[js];
  385.                                 dataj.yj = data.py[js];
  386.                                 dataj.zj = data.pz[js];
  387.                                 dataj.s = data.pbi[js];
  388.  
  389.                                 /* ground loop */
  390.                                 for (ip = 1; ip <= gnd.ksymp; ip++)
  391.                                 {
  392.                                         dataj.ipgnd = ip;
  393.  
  394.                                         if (((ipch == j) || (icgo != 0)) && (ip != 2))
  395.                                         {
  396.                                                 if (icgo <= 0)
  397.                                                 {
  398.                                                         pcint (
  399.                                                             xi,
  400.                                                             yi,
  401.                                                             zi,
  402.                                                             cabi,
  403.                                                             sabi,
  404.                                                             salpi,
  405.                                                             emel);
  406.  
  407.                                                         pyl = PI * data.si[i] * fsign;
  408.                                                         pxl = sinl (pyl);
  409.                                                         pyl = cosl (pyl);
  410.                                                         dataj.exc = emel[8] * fsign;
  411.  
  412.                                                         trio (i + 1);
  413.  
  414.                                                         il = i - ncw;
  415.                                                         if (i < data.np)
  416.                                                                 il += (il / data.np) * 2 *
  417.                                                                       data.mp;
  418.  
  419.                                                         if (itrp == 0)
  420.                                                                 cw[k + il * nrow] +=
  421.                                                                     dataj.exc *
  422.                                                                     (segj.ax[jsnox] +
  423.                                                                      segj.bx[jsnox] * pxl +
  424.                                                                      segj.cx[jsnox] * pyl);
  425.                                                         else
  426.                                                                 cw[il + k * nrow] +=
  427.                                                                     dataj.exc *
  428.                                                                     (segj.ax[jsnox] +
  429.                                                                      segj.bx[jsnox] * pxl +
  430.                                                                      segj.cx[jsnox] * pyl);
  431.  
  432.                                                 } /* if( icgo <= 0 ) */
  433.  
  434.                                                 if (itrp == 0)
  435.                                                 {
  436.                                                         cm[k + (jl - 1) * nrow] = emel[icgo];
  437.                                                         cm[k + jl * nrow] = emel[icgo + 4];
  438.                                                 }
  439.                                                 else
  440.                                                 {
  441.                                                         cm[(jl - 1) + k * nrow] = emel[icgo];
  442.                                                         cm[jl + k * nrow] = emel[icgo + 4];
  443.                                                 }
  444.  
  445.                                                 icgo++;
  446.                                                 if (icgo == 4)
  447.                                                         icgo = 0;
  448.  
  449.                                                 continue;
  450.  
  451.                                         } /* if( ((ipch == (j+1)) || (icgo != 0)) && (ip != 2)
  452.                                              ) */
  453.  
  454.                                         unere (xi, yi, zi);
  455.  
  456.                                         /* normal fill */
  457.                                         if (itrp == 0)
  458.                                         {
  459.                                                 cm[k + (jl - 1) * nrow] += dataj.exk * cabi +
  460.                                                                            dataj.eyk * sabi +
  461.                                                                            dataj.ezk * salpi;
  462.                                                 cm[k + jl * nrow] += dataj.exs * cabi +
  463.                                                                      dataj.eys * sabi +
  464.                                                                      dataj.ezs * salpi;
  465.                                                 continue;
  466.                                         }
  467.  
  468.                                         /* transposed fill */
  469.                                         cm[(jl - 1) + k * nrow] += dataj.exk * cabi +
  470.                                                                    dataj.eyk * sabi +
  471.                                                                    dataj.ezk * salpi;
  472.                                         cm[jl + k * nrow] += dataj.exs * cabi +
  473.                                                              dataj.eys * sabi +
  474.                                                              dataj.ezs * salpi;
  475.  
  476.                                 } /* for( ip = 1; ip <= gnd.ksymp; ip++ ) */
  477.  
  478.                         } /* for( j = j1; j <= j2; j++ ) */
  479.  
  480.                 } /* for( i = i1-1; i < i2; i++ ) */
  481.  
  482.         } /* if( itrp >= 0) */
  483.  
  484.         return;
  485. }
  486.  
  487. /*-----------------------------------------------------------------------*/
  488.  
  489. /* cmws computes matrix elements for wire-surface interactions */
  490. void cmws (
  491.     int j, int i1, int i2, complex double *cm, int nr, complex double *cw, int nw, int itrp)
  492. {
  493.         int ipr, i, ipatch, ik, js = 0, ij, jx;
  494.         double xi, yi, zi, tx, ty, tz;
  495.         complex double etk, ets, etc;
  496.         i = nw;
  497.         j--;
  498.         dataj.s = data.si[j];
  499.         dataj.b = data.bi[j];
  500.         dataj.xj = data.x[j];
  501.         dataj.yj = data.y[j];
  502.         dataj.zj = data.z[j];
  503.         dataj.cabj = data.cab[j];
  504.         dataj.sabj = data.sab[j];
  505.         dataj.salpj = data.salp[j];
  506.  
  507.         /* observation loop */
  508.         ipr = -1;
  509.         for (i = i1; i <= i2; i++)
  510.         {
  511.                 ipr++;
  512.                 ipatch = (i + 1) / 2;
  513.                 ik = i - (i / 2) * 2;
  514.  
  515.                 if ((ik != 0) || (ipr == 0))
  516.                 {
  517.                         js = ipatch - 1;
  518.                         xi = data.px[js];
  519.                         yi = data.py[js];
  520.                         zi = data.pz[js];
  521.                         hsfld (xi, yi, zi, 0.);
  522.  
  523.                         if (ik != 0)
  524.                         {
  525.                                 tx = data.t2x[js];
  526.                                 ty = data.t2y[js];
  527.                                 tz = data.t2z[js];
  528.                         }
  529.                         else
  530.                         {
  531.                                 tx = data.t1x[js];
  532.                                 ty = data.t1y[js];
  533.                                 tz = data.t1z[js];
  534.                         }
  535.  
  536.                 } /* if( (ik != 0) || (ipr == 0) ) */
  537.                 else
  538.                 {
  539.                         tx = data.t1x[js];
  540.                         ty = data.t1y[js];
  541.                         tz = data.t1z[js];
  542.  
  543.                 } /* if( (ik != 0) || (ipr == 0) ) */
  544.  
  545.                 etk = -(dataj.exk * tx + dataj.eyk * ty + dataj.ezk * tz) * data.psalp[js];
  546.                 ets = -(dataj.exs * tx + dataj.eys * ty + dataj.ezs * tz) * data.psalp[js];
  547.                 etc = -(dataj.exc * tx + dataj.eyc * ty + dataj.ezc * tz) * data.psalp[js];
  548.  
  549.                 /* fill matrix elements.  element locations */
  550.                 /* determined by connection data. */
  551.  
  552.                 /* normal fill */
  553.                 if (itrp == 0)
  554.                 {
  555.                         for (ij = 0; ij < segj.jsno; ij++)
  556.                         {
  557.                                 jx = segj.jco[ij] - 1;
  558.                                 cm[ipr + jx * nr] +=
  559.                                     etk * segj.ax[ij] + ets * segj.bx[ij] + etc * segj.cx[ij];
  560.                         }
  561.  
  562.                         continue;
  563.                 } /* if( itrp == 0) */
  564.  
  565.                 /* transposed fill */
  566.                 if (itrp != 2)
  567.                 {
  568.                         for (ij = 0; ij < segj.jsno; ij++)
  569.                         {
  570.                                 jx = segj.jco[ij] - 1;
  571.                                 cm[jx + ipr * nr] +=
  572.                                     etk * segj.ax[ij] + ets * segj.bx[ij] + etc * segj.cx[ij];
  573.                         }
  574.  
  575.                         continue;
  576.                 } /* if( itrp != 2) */
  577.  
  578.                 /* transposed fill - c(ws) and d(ws)prime (=cw) */
  579.                 for (ij = 0; ij < segj.jsno; ij++)
  580.                 {
  581.                         jx = segj.jco[ij] - 1;
  582.                         if (jx < nr)
  583.                                 cm[jx + ipr * nr] +=
  584.                                     etk * segj.ax[ij] + ets * segj.bx[ij] + etc * segj.cx[ij];
  585.                         else
  586.                         {
  587.                                 jx -= nr;
  588.                                 cw[jx + ipr * nr] +=
  589.                                     etk * segj.ax[ij] + ets * segj.bx[ij] + etc * segj.cx[ij];
  590.                         }
  591.                 } /* for( ij = 0; ij < segj.jsno; ij++ ) */
  592.  
  593.         } /* for( i = i1; i <= i2; i++ ) */
  594.  
  595.         return;
  596. }
  597.  
  598. /*-----------------------------------------------------------------------*/
  599.  
  600. /* cmww computes matrix elements for wire-wire interactions */
  601. void cmww (
  602.     int j, int i1, int i2, complex double *cm, int nr, complex double *cw, int nw, int itrp)
  603. {
  604.         int ipr, iprx, i, ij, jx;
  605.         double xi, yi, zi, ai, cabi, sabi, salpi;
  606.         complex double etk, ets, etc;
  607.  
  608.         /* set source segment parameters */
  609.         jx = j;
  610.         j--;
  611.         dataj.s = data.si[j];
  612.         dataj.b = data.bi[j];
  613.         dataj.xj = data.x[j];
  614.         dataj.yj = data.y[j];
  615.         dataj.zj = data.z[j];
  616.         dataj.cabj = data.cab[j];
  617.         dataj.sabj = data.sab[j];
  618.         dataj.salpj = data.salp[j];
  619.  
  620.         /* decide whether ext. t.w. approx. can be used */
  621.         if (dataj.iexk != 0)
  622.         {
  623.                 ipr = data.icon1[j];
  624.                 if (ipr > PCHCON)
  625.                         dataj.ind1 = 0;
  626.                 else if (ipr < 0)
  627.                 {
  628.                         ipr = -ipr;
  629.                         iprx = ipr - 1;
  630.  
  631.                         if (-data.icon1[iprx] != jx)
  632.                                 dataj.ind1 = 2;
  633.                         else
  634.                         {
  635.                                 xi = fabsl (
  636.                                     dataj.cabj * data.cab[iprx] + dataj.sabj * data.sab[iprx] +
  637.                                     dataj.salpj * data.salp[iprx]);
  638.                                 if ((xi < 0.999999) ||
  639.                                     (fabsl (data.bi[iprx] / dataj.b - 1.) > 1.e-6))
  640.                                         dataj.ind1 = 2;
  641.                                 else
  642.                                         dataj.ind1 = 0;
  643.  
  644.                         } /* if( -data.icon1[iprx] != jx ) */
  645.  
  646.                 } /* if( ipr < 0 ) */
  647.                 else
  648.                 {
  649.                         iprx = ipr - 1;
  650.                         if (ipr == 0)
  651.                                 dataj.ind1 = 1;
  652.                         else
  653.                         {
  654.                                 if (ipr != jx)
  655.                                 {
  656.                                         if (data.icon2[iprx] != jx)
  657.                                                 dataj.ind1 = 2;
  658.                                         else
  659.                                         {
  660.                                                 xi = fabsl (
  661.                                                     dataj.cabj * data.cab[iprx] +
  662.                                                     dataj.sabj * data.sab[iprx] +
  663.                                                     dataj.salpj * data.salp[iprx]);
  664.                                                 if ((xi < 0.999999) ||
  665.                                                     (fabsl (data.bi[iprx] / dataj.b - 1.) >
  666.                                                      1.e-6))
  667.                                                         dataj.ind1 = 2;
  668.                                                 else
  669.                                                         dataj.ind1 = 0;
  670.  
  671.                                         } /* if( data.icon2[iprx] != jx ) */
  672.  
  673.                                 } /* if( ipr != jx ) */
  674.                                 else if (
  675.                                     dataj.cabj * dataj.cabj + dataj.sabj * dataj.sabj > 1.e-8)
  676.                                         dataj.ind1 = 2;
  677.                                 else
  678.                                         dataj.ind1 = 0;
  679.  
  680.                         } /* if( ipr == 0 ) */
  681.  
  682.                 } /* if( ipr < 0 ) */
  683.  
  684.                 ipr = data.icon2[j];
  685.                 if (ipr > PCHCON)
  686.                         dataj.ind2 = 2;
  687.                 else if (ipr < 0)
  688.                 {
  689.                         ipr = -ipr;
  690.                         iprx = ipr - 1;
  691.                         if (-data.icon2[iprx] != jx)
  692.                                 dataj.ind2 = 2;
  693.                         else
  694.                         {
  695.                                 xi = fabsl (
  696.                                     dataj.cabj * data.cab[iprx] + dataj.sabj * data.sab[iprx] +
  697.                                     dataj.salpj * data.salp[iprx]);
  698.                                 if ((xi < 0.999999) ||
  699.                                     (fabsl (data.bi[iprx] / dataj.b - 1.) > 1.e-6))
  700.                                         dataj.ind2 = 2;
  701.                                 else
  702.                                         dataj.ind2 = 0;
  703.  
  704.                         } /* if( -data.icon1[iprx] != jx ) */
  705.  
  706.                 } /* if( ipr < 0 ) */
  707.                 else
  708.                 {
  709.                         iprx = ipr - 1;
  710.                         if (ipr == 0)
  711.                                 dataj.ind2 = 1;
  712.                         else
  713.                         {
  714.                                 if (ipr != jx)
  715.                                 {
  716.                                         if (data.icon1[iprx] != jx)
  717.                                                 dataj.ind2 = 2;
  718.                                         else
  719.                                         {
  720.                                                 xi = fabsl (
  721.                                                     dataj.cabj * data.cab[iprx] +
  722.                                                     dataj.sabj * data.sab[iprx] +
  723.                                                     dataj.salpj * data.salp[iprx]);
  724.                                                 if ((xi < 0.999999) ||
  725.                                                     (fabsl (data.bi[iprx] / dataj.b - 1.) >
  726.                                                      1.e-6))
  727.                                                         dataj.ind2 = 2;
  728.                                                 else
  729.                                                         dataj.ind2 = 0;
  730.  
  731.                                         } /* if( data.icon2[iprx] != jx ) */
  732.  
  733.                                 } /* if( ipr != jx ) */
  734.                                 else if (
  735.                                     dataj.cabj * dataj.cabj + dataj.sabj * dataj.sabj > 1.e-8)
  736.                                         dataj.ind2 = 2;
  737.                                 else
  738.                                         dataj.ind2 = 0;
  739.  
  740.                         } /* if( ipr == 0 ) */
  741.  
  742.                 } /* if( ipr < 0 ) */
  743.  
  744.         } /* if( dataj.iexk != 0) */
  745.  
  746.         /* observation loop */
  747.         ipr = -1;
  748.         for (i = i1 - 1; i < i2; i++)
  749.         {
  750.                 ipr++;
  751.                 ij = i - j;
  752.                 xi = data.x[i];
  753.                 yi = data.y[i];
  754.                 zi = data.z[i];
  755.                 ai = data.bi[i];
  756.                 cabi = data.cab[i];
  757.                 sabi = data.sab[i];
  758.                 salpi = data.salp[i];
  759.  
  760.                 efld (xi, yi, zi, ai, ij);
  761.  
  762.                 etk = dataj.exk * cabi + dataj.eyk * sabi + dataj.ezk * salpi;
  763.                 ets = dataj.exs * cabi + dataj.eys * sabi + dataj.ezs * salpi;
  764.                 etc = dataj.exc * cabi + dataj.eyc * sabi + dataj.ezc * salpi;
  765.  
  766.                 /* fill matrix elements. element locations */
  767.                 /* determined by connection data. */
  768.  
  769.                 /* normal fill */
  770.                 if (itrp == 0)
  771.                 {
  772.                         for (ij = 0; ij < segj.jsno; ij++)
  773.                         {
  774.                                 jx = segj.jco[ij] - 1;
  775.                                 cm[ipr + jx * nr] +=
  776.                                     etk * segj.ax[ij] + ets * segj.bx[ij] + etc * segj.cx[ij];
  777.                         }
  778.                         continue;
  779.                 }
  780.  
  781.                 /* transposed fill */
  782.                 if (itrp != 2)
  783.                 {
  784.                         for (ij = 0; ij < segj.jsno; ij++)
  785.                         {
  786.                                 jx = segj.jco[ij] - 1;
  787.                                 cm[jx + ipr * nr] +=
  788.                                     etk * segj.ax[ij] + ets * segj.bx[ij] + etc * segj.cx[ij];
  789.                         }
  790.                         continue;
  791.                 }
  792.  
  793.                 /* trans. fill for c(ww) - test for elements for d(ww)prime.  (=cw) */
  794.                 for (ij = 0; ij < segj.jsno; ij++)
  795.                 {
  796.                         jx = segj.jco[ij] - 1;
  797.                         if (jx < nr)
  798.                                 cm[jx + ipr * nr] +=
  799.                                     etk * segj.ax[ij] + ets * segj.bx[ij] + etc * segj.cx[ij];
  800.                         else
  801.                         {
  802.                                 jx -= nr;
  803.                                 cw[jx * ipr * nw] +=
  804.                                     etk * segj.ax[ij] + ets * segj.bx[ij] + etc * segj.cx[ij];
  805.                         }
  806.  
  807.                 } /* for( ij = 0; ij < segj.jsno; ij++ ) */
  808.  
  809.         } /* for( i = i1-1; i < i2; i++ ) */
  810.  
  811.         return;
  812. }
  813.  
  814. /*-----------------------------------------------------------------------*/
  815.  
  816. /* etmns fills the array e with the negative of the */
  817. /* electric field incident on the structure. e is the */
  818. /* right hand side of the matrix equation. */
  819. void etmns (
  820.     double p1,
  821.     double p2,
  822.     double p3,
  823.     double p4,
  824.     double p5,
  825.     double p6,
  826.     int ipr,
  827.     complex double *e)
  828. {
  829.         int i, is, i1, i2 = 0, neq;
  830.         double cth, sth, cph, sph, cet, set, pxl, pyl, pzl, wx;
  831.         double wy, wz, qx, qy, qz, arg, ds, dsh, rs, r;
  832.         complex double cx, cy, cz, er, et, ezh, erh, rrv = CPLX_00, rrh = CPLX_00, tt1, tt2;
  833.  
  834.         neq = data.n + 2 * data.m;
  835.         vsorc.nqds = 0;
  836.  
  837.         /* applied field of voltage sources for transmitting case */
  838.         if ((ipr == 0) || (ipr == 5))
  839.         {
  840.                 for (i = 0; i < neq; i++)
  841.                         e[i] = CPLX_00;
  842.  
  843.                 if (vsorc.nsant != 0)
  844.                 {
  845.                         for (i = 0; i < vsorc.nsant; i++)
  846.                         {
  847.                                 is = vsorc.isant[i] - 1;
  848.                                 e[is] = -vsorc.vsant[i] / (data.si[is] * data.wlam);
  849.                         }
  850.                 }
  851.  
  852.                 if (vsorc.nvqd == 0)
  853.                         return;
  854.  
  855.                 for (i = 0; i < vsorc.nvqd; i++)
  856.                 {
  857.                         is = vsorc.ivqd[i];
  858.                         qdsrc (is, vsorc.vqd[i], e);
  859.                 }
  860.                 return;
  861.  
  862.         } /* if( (ipr == 0) || (ipr == 5) ) */
  863.  
  864.         /* incident plane wave, linearly polarized. */
  865.         if (ipr <= 3)
  866.         {
  867.                 cth = cosl (p1);
  868.                 sth = sinl (p1);
  869.                 cph = cosl (p2);
  870.                 sph = sinl (p2);
  871.                 cet = cosl (p3);
  872.                 set = sinl (p3);
  873.                 pxl = cth * cph * cet - sph * set;
  874.                 pyl = cth * sph * cet + cph * set;
  875.                 pzl = -sth * cet;
  876.                 wx = -sth * cph;
  877.                 wy = -sth * sph;
  878.                 wz = -cth;
  879.                 qx = wy * pzl - wz * pyl;
  880.                 qy = wz * pxl - wx * pzl;
  881.                 qz = wx * pyl - wy * pxl;
  882.  
  883.                 if (gnd.ksymp != 1)
  884.                 {
  885.                         if (gnd.iperf != 1)
  886.                         {
  887.                                 rrv = csqrtl (1. - gnd.zrati * gnd.zrati * sth * sth);
  888.                                 rrh = gnd.zrati * cth;
  889.                                 rrh = (rrh - rrv) / (rrh + rrv);
  890.                                 rrv = gnd.zrati * rrv;
  891.                                 rrv = -(cth - rrv) / (cth + rrv);
  892.                         }
  893.                         else
  894.                         {
  895.                                 rrv = -CPLX_10;
  896.                                 rrh = -CPLX_10;
  897.                         } /* if( gnd.iperf != 1) */
  898.  
  899.                 } /* if( gnd.ksymp != 1) */
  900.  
  901.                 if (ipr == 1)
  902.                 {
  903.                         if (data.n != 0)
  904.                         {
  905.                                 for (i = 0; i < data.n; i++)
  906.                                 {
  907.                                         arg = -TP * (wx * data.x[i] + wy * data.y[i] +
  908.                                                      wz * data.z[i]);
  909.                                         e[i] = -(pxl * data.cab[i] + pyl * data.sab[i] +
  910.                                                  pzl * data.salp[i]) *
  911.                                                cmplx (cosl (arg), sinl (arg));
  912.                                 }
  913.  
  914.                                 if (gnd.ksymp != 1)
  915.                                 {
  916.                                         tt1 = (pyl * cph - pxl * sph) * (rrh - rrv);
  917.                                         cx = rrv * pxl - tt1 * sph;
  918.                                         cy = rrv * pyl + tt1 * cph;
  919.                                         cz = -rrv * pzl;
  920.  
  921.                                         for (i = 0; i < data.n; i++)
  922.                                         {
  923.                                                 arg = -TP * (wx * data.x[i] + wy * data.y[i] -
  924.                                                              wz * data.z[i]);
  925.                                                 e[i] = e[i] -
  926.                                                        (cx * data.cab[i] + cy * data.sab[i] +
  927.                                                         cz * data.salp[i]) *
  928.                                                            cmplx (cosl (arg), sinl (arg));
  929.                                         }
  930.  
  931.                                 } /* if( gnd.ksymp != 1) */
  932.  
  933.                         } /* if( data.n != 0) */
  934.  
  935.                         if (data.m == 0)
  936.                                 return;
  937.  
  938.                         i = -1;
  939.                         i1 = data.n - 2;
  940.                         for (is = 0; is < data.m; is++)
  941.                         {
  942.                                 i++;
  943.                                 i1 += 2;
  944.                                 i2 = i1 + 1;
  945.                                 arg = -TP *
  946.                                       (wx * data.px[i] + wy * data.py[i] + wz * data.pz[i]);
  947.                                 tt1 = cmplx (cosl (arg), sinl (arg)) * data.psalp[i] * RETA;
  948.                                 e[i2] =
  949.                                     (qx * data.t1x[i] + qy * data.t1y[i] + qz * data.t1z[i]) *
  950.                                     tt1;
  951.                                 e[i1] =
  952.                                     (qx * data.t2x[i] + qy * data.t2y[i] + qz * data.t2z[i]) *
  953.                                     tt1;
  954.                         }
  955.  
  956.                         if (gnd.ksymp == 1)
  957.                                 return;
  958.  
  959.                         tt1 = (qy * cph - qx * sph) * (rrv - rrh);
  960.                         cx = -(rrh * qx - tt1 * sph);
  961.                         cy = -(rrh * qy + tt1 * cph);
  962.                         cz = rrh * qz;
  963.  
  964.                         i = -1;
  965.                         i1 = data.n - 2;
  966.                         for (is = 0; is < data.m; is++)
  967.                         {
  968.                                 i++;
  969.                                 i1 += 2;
  970.                                 i2 = i1 + 1;
  971.                                 arg = -TP *
  972.                                       (wx * data.px[i] + wy * data.py[i] - wz * data.pz[i]);
  973.                                 tt1 = cmplx (cosl (arg), sinl (arg)) * data.psalp[i] * RETA;
  974.                                 e[i2] = e[i2] + (cx * data.t1x[i] + cy * data.t1y[i] +
  975.                                                  cz * data.t1z[i]) *
  976.                                                     tt1;
  977.                                 e[i1] = e[i1] + (cx * data.t2x[i] + cy * data.t2y[i] +
  978.                                                  cz * data.t2z[i]) *
  979.                                                     tt1;
  980.                         }
  981.                         return;
  982.  
  983.                 } /* if( ipr == 1) */
  984.  
  985.                 /* incident plane wave, elliptic polarization. */
  986.                 tt1 = -(CPLX_01) *p6;
  987.                 if (ipr == 3)
  988.                         tt1 = -tt1;
  989.  
  990.                 if (data.n != 0)
  991.                 {
  992.                         cx = pxl + tt1 * qx;
  993.                         cy = pyl + tt1 * qy;
  994.                         cz = pzl + tt1 * qz;
  995.  
  996.                         for (i = 0; i < data.n; i++)
  997.                         {
  998.                                 arg = -TP * (wx * data.x[i] + wy * data.y[i] + wz * data.z[i]);
  999.                                 e[i] = -(cx * data.cab[i] + cy * data.sab[i] +
  1000.                                          cz * data.salp[i]) *
  1001.                                        cmplx (cosl (arg), sinl (arg));
  1002.                         }
  1003.  
  1004.                         if (gnd.ksymp != 1)
  1005.                         {
  1006.                                 tt2 = (cy * cph - cx * sph) * (rrh - rrv);
  1007.                                 cx = rrv * cx - tt2 * sph;
  1008.                                 cy = rrv * cy + tt2 * cph;
  1009.                                 cz = -rrv * cz;
  1010.  
  1011.                                 for (i = 0; i < data.n; i++)
  1012.                                 {
  1013.                                         arg = -TP * (wx * data.x[i] + wy * data.y[i] -
  1014.                                                      wz * data.z[i]);
  1015.                                         e[i] = e[i] - (cx * data.cab[i] + cy * data.sab[i] +
  1016.                                                        cz * data.salp[i]) *
  1017.                                                           cmplx (cosl (arg), sinl (arg));
  1018.                                 }
  1019.  
  1020.                         } /* if( gnd.ksymp != 1) */
  1021.  
  1022.                 } /* if( n != 0) */
  1023.  
  1024.                 if (data.m == 0)
  1025.                         return;
  1026.  
  1027.                 cx = qx - tt1 * pxl;
  1028.                 cy = qy - tt1 * pyl;
  1029.                 cz = qz - tt1 * pzl;
  1030.  
  1031.                 i = -1;
  1032.                 i1 = data.n - 2;
  1033.                 for (is = 0; is < data.m; is++)
  1034.                 {
  1035.                         i++;
  1036.                         i1 += 2;
  1037.                         i2 = i1 + 1;
  1038.                         arg = -TP * (wx * data.px[i] + wy * data.py[i] + wz * data.pz[i]);
  1039.                         tt2 = cmplx (cosl (arg), sinl (arg)) * data.psalp[i] * RETA;
  1040.                         e[i2] = (cx * data.t1x[i] + cy * data.t1y[i] + cz * data.t1z[i]) * tt2;
  1041.                         e[i1] = (cx * data.t2x[i] + cy * data.t2y[i] + cz * data.t2z[i]) * tt2;
  1042.                 }
  1043.  
  1044.                 if (gnd.ksymp == 1)
  1045.                         return;
  1046.  
  1047.                 tt1 = (cy * cph - cx * sph) * (rrv - rrh);
  1048.                 cx = -(rrh * cx - tt1 * sph);
  1049.                 cy = -(rrh * cy + tt1 * cph);
  1050.                 cz = rrh * cz;
  1051.  
  1052.                 i = -1;
  1053.                 i1 = data.n - 2;
  1054.                 for (is = 0; is < data.m; is++)
  1055.                 {
  1056.                         i++;
  1057.                         i1 += 2;
  1058.                         i2 = i1 + 1;
  1059.                         arg = -TP * (wx * data.px[i] + wy * data.py[i] - wz * data.pz[i]);
  1060.                         tt1 = cmplx (cosl (arg), sinl (arg)) * data.psalp[i] * RETA;
  1061.                         e[i2] = e[i2] +
  1062.                                 (cx * data.t1x[i] + cy * data.t1y[i] + cz * data.t1z[i]) * tt1;
  1063.                         e[i1] = e[i1] +
  1064.                                 (cx * data.t2x[i] + cy * data.t2y[i] + cz * data.t2z[i]) * tt1;
  1065.                 }
  1066.  
  1067.                 return;
  1068.  
  1069.         } /* if( ipr <= 3) */
  1070.  
  1071.         /* incident field of an elementary current source. */
  1072.         wz = cosl (p4);
  1073.         wx = wz * cosl (p5);
  1074.         wy = wz * sinl (p5);
  1075.         wz = sinl (p4);
  1076.         ds = p6 * 59.958;
  1077.         dsh = p6 / (2. * TP);
  1078.  
  1079.         is = 0;
  1080.         i1 = data.n - 2;
  1081.         for (i = 0; i < data.npm; i++)
  1082.         {
  1083.                 if (i >= data.n)
  1084.                 {
  1085.                         i1 += 2;
  1086.                         i2 = i1 + 1;
  1087.                         pxl = data.px[is] - p1;
  1088.                         pyl = data.py[is] - p2;
  1089.                         pzl = data.pz[is] - p3;
  1090.                 }
  1091.                 else
  1092.                 {
  1093.                         pxl = data.x[i] - p1;
  1094.                         pyl = data.y[i] - p2;
  1095.                         pzl = data.z[i] - p3;
  1096.                 }
  1097.  
  1098.                 rs = pxl * pxl + pyl * pyl + pzl * pzl;
  1099.                 if (rs < 1.0e-30)
  1100.                         continue;
  1101.  
  1102.                 r = sqrtl (rs);
  1103.                 pxl = pxl / r;
  1104.                 pyl = pyl / r;
  1105.                 pzl = pzl / r;
  1106.                 cth = pxl * wx + pyl * wy + pzl * wz;
  1107.                 sth = sqrtl (1. - cth * cth);
  1108.                 qx = pxl - wx * cth;
  1109.                 qy = pyl - wy * cth;
  1110.                 qz = pzl - wz * cth;
  1111.  
  1112.                 arg = sqrtl (qx * qx + qy * qy + qz * qz);
  1113.                 if (arg >= 1.e-30)
  1114.                 {
  1115.                         qx = qx / arg;
  1116.                         qy = qy / arg;
  1117.                         qz = qz / arg;
  1118.                 }
  1119.                 else
  1120.                 {
  1121.                         qx = 1.;
  1122.                         qy = 0.;
  1123.                         qz = 0.;
  1124.  
  1125.                 } /* if( arg >= 1.e-30) */
  1126.  
  1127.                 arg = -TP * r;
  1128.                 tt1 = cmplx (cosl (arg), sinl (arg));
  1129.  
  1130.                 if (i < data.n)
  1131.                 {
  1132.                         tt2 = cmplx (1.0, -1.0 / (r * TP)) / rs;
  1133.                         er = ds * tt1 * tt2 * cth;
  1134.                         et = .5 * ds * tt1 * ((CPLX_01) *TP / r + tt2) * sth;
  1135.                         ezh = er * cth - et * sth;
  1136.                         erh = er * sth + et * cth;
  1137.                         cx = ezh * wx + erh * qx;
  1138.                         cy = ezh * wy + erh * qy;
  1139.                         cz = ezh * wz + erh * qz;
  1140.                         e[i] = -(cx * data.cab[i] + cy * data.sab[i] + cz * data.salp[i]);
  1141.                 }
  1142.                 else
  1143.                 {
  1144.                         pxl = wy * qz - wz * qy;
  1145.                         pyl = wz * qx - wx * qz;
  1146.                         pzl = wx * qy - wy * qx;
  1147.                         tt2 = dsh * tt1 * cmplx (1. / r, TP) / r * sth * data.psalp[is];
  1148.                         cx = tt2 * pxl;
  1149.                         cy = tt2 * pyl;
  1150.                         cz = tt2 * pzl;
  1151.                         e[i2] = cx * data.t1x[is] + cy * data.t1y[is] + cz * data.t1z[is];
  1152.                         e[i1] = cx * data.t2x[is] + cy * data.t2y[is] + cz * data.t2z[is];
  1153.                         is++;
  1154.                 } /* if( i < data.n) */
  1155.  
  1156.         } /* for( i = 0; i < npm; i++ ) */
  1157.  
  1158.         return;
  1159. }
  1160.  
  1161. /*-----------------------------------------------------------------------*/
  1162.  
  1163. /* subroutine to factor a matrix into a unit lower triangular matrix */
  1164. /* and an upper triangular matrix using the gauss-doolittle algorithm */
  1165. /* presented on pages 411-416 of a. ralston--a first course in */
  1166. /* numerical analysis.  comments below refer to comments in ralstons */
  1167. /* text.    (matrix transposed.) */
  1168.  
  1169. void factr (int n, complex double *a, int *ip, int ndim)
  1170. {
  1171.         int r, rm1, rp1, pj, pr, iflg, k, j, jp1, i;
  1172.         double dmax, elmag;
  1173.         complex double arj, *scm = NULL;
  1174.  
  1175.         /* Allocate to scratch memory */
  1176.         mem_alloc ((void *) &scm, data.np2m * sizeof (complex double));
  1177.  
  1178.         /* Un-transpose the matrix for Gauss elimination */
  1179.         for (i = 1; i < n; i++)
  1180.                 for (j = 0; j < i; j++)
  1181.                 {
  1182.                         arj = a[i + j * ndim];
  1183.                         a[i + j * ndim] = a[j + i * ndim];
  1184.                         a[j + i * ndim] = arj;
  1185.                 }
  1186.  
  1187.         iflg = FALSE;
  1188.         /* step 1 */
  1189.         for (r = 0; r < n; r++)
  1190.         {
  1191.                 for (k = 0; k < n; k++)
  1192.                         scm[k] = a[k + r * ndim];
  1193.  
  1194.                 /* steps 2 and 3 */
  1195.                 rm1 = r;
  1196.                 if (rm1 > 0)
  1197.                 {
  1198.                         for (j = 0; j < rm1; j++)
  1199.                         {
  1200.                                 pj = ip[j] - 1;
  1201.                                 arj = scm[pj];
  1202.                                 a[j + r * ndim] = arj;
  1203.                                 scm[pj] = scm[j];
  1204.                                 jp1 = j + 1;
  1205.  
  1206.                                 for (i = jp1; i < n; i++)
  1207.                                         scm[i] -= a[i + j * ndim] * arj;
  1208.  
  1209.                         } /* for( j = 0; j < rm1; j++ ) */
  1210.  
  1211.                 } /* if( rm1 >= 0.) */
  1212.  
  1213.                 /* step 4 */
  1214.                 dmax = creal (scm[r] * conjl (scm[r]));
  1215.  
  1216.                 rp1 = r + 1;
  1217.                 ip[r] = rp1;
  1218.                 if (rp1 < n)
  1219.                 {
  1220.                         for (i = rp1; i < n; i++)
  1221.                         {
  1222.                                 elmag = creal (scm[i] * conjl (scm[i]));
  1223.                                 if (elmag >= dmax)
  1224.                                 {
  1225.                                         dmax = elmag;
  1226.                                         ip[r] = i + 1;
  1227.                                 }
  1228.                         }
  1229.                 } /* if( rp1 < n) */
  1230.  
  1231.                 if (dmax < 1.e-10)
  1232.                         iflg = TRUE;
  1233.  
  1234.                 pr = ip[r] - 1;
  1235.                 a[r + r * ndim] = scm[pr];
  1236.                 scm[pr] = scm[r];
  1237.  
  1238.                 /* step 5 */
  1239.                 if (rp1 < n)
  1240.                 {
  1241.                         arj = 1. / a[r + r * ndim];
  1242.  
  1243.                         for (i = rp1; i < n; i++)
  1244.                                 a[i + r * ndim] = scm[i] * arj;
  1245.                 }
  1246.  
  1247.                 if (iflg == TRUE)
  1248.                 {
  1249.                         fprintf (output_fp, "\n  PIVOT(%d)= %16.8E", r, dmax);
  1250.                         iflg = FALSE;
  1251.                 }
  1252.  
  1253.         } /* for( r=0; r < n; r++ ) */
  1254.  
  1255.         free_ptr ((void *) &scm);
  1256.  
  1257.         return;
  1258. }
  1259.  
  1260. /*-----------------------------------------------------------------------*/
  1261.  
  1262. /* factrs, for symmetric structure, transforms submatricies to form */
  1263. /* matricies of the symmetric modes and calls routine to factor */
  1264. /* matricies.  if no symmetry, the routine is called to factor the */
  1265. /* complete matrix. */
  1266. void factrs (int np, int nrow, complex double *a, int *ip)
  1267. {
  1268.         int kk, ka;
  1269.  
  1270.         smat.nop = nrow / np;
  1271.         for (kk = 0; kk < smat.nop; kk++)
  1272.         {
  1273.                 ka = kk * np;
  1274.                 factr (np, &a[ka], &ip[ka], nrow);
  1275.         }
  1276.         return;
  1277. }
  1278.  
  1279. /*-----------------------------------------------------------------------*/
  1280.  
  1281. /* fblock sets parameters for out-of-core */
  1282. /* solution for the primary matrix (a) */
  1283. void fblock (int nrow, int ncol, int imax, int ipsym)
  1284. {
  1285.         int i, j, k, ka, kk;
  1286.         double phaz, arg;
  1287.         complex double deter;
  1288.  
  1289.         if (nrow * ncol <= imax)
  1290.         {
  1291.                 matpar.npblk = nrow;
  1292.                 matpar.nlast = nrow;
  1293.                 matpar.imat = nrow * ncol;
  1294.  
  1295.                 if (nrow == ncol)
  1296.                 {
  1297.                         matpar.icase = 1;
  1298.                         return;
  1299.                 }
  1300.                 else
  1301.                         matpar.icase = 2;
  1302.  
  1303.         } /* if( nrow*ncol <= imax) */
  1304.  
  1305.         smat.nop = ncol / nrow;
  1306.         if (smat.nop * nrow != ncol)
  1307.         {
  1308.                 fprintf (output_fp, "\n  SYMMETRY ERROR - NROW: %d NCOL: %d", nrow, ncol);
  1309.                 stop (-1);
  1310.         }
  1311.  
  1312.         /* set up smat.ssx matrix for rotational symmetry. */
  1313.         if (ipsym <= 0)
  1314.         {
  1315.                 phaz = TP / smat.nop;
  1316.  
  1317.                 for (i = 1; i < smat.nop; i++)
  1318.                 {
  1319.                         for (j = i; j < smat.nop; j++)
  1320.                         {
  1321.                                 arg = phaz * (double) i * (double) j;
  1322.                                 smat.ssx[i + j * smat.nop] = cmplx (cosl (arg), sinl (arg));
  1323.                                 smat.ssx[j + i * smat.nop] = smat.ssx[i + j * smat.nop];
  1324.                         }
  1325.                 }
  1326.                 return;
  1327.  
  1328.         } /* if( ipsym <= 0) */
  1329.  
  1330.         /* set up smat.ssx matrix for plane symmetry */
  1331.         kk = 1;
  1332.         smat.ssx[0] = CPLX_10;
  1333.  
  1334.         k = 2;
  1335.         for (ka = 1; k != smat.nop; ka++)
  1336.                 k *= 2;
  1337.  
  1338.         for (k = 0; k < ka; k++)
  1339.         {
  1340.                 for (i = 0; i < kk; i++)
  1341.                 {
  1342.                         for (j = 0; j < kk; j++)
  1343.                         {
  1344.                                 deter = smat.ssx[i + j * smat.nop];
  1345.                                 smat.ssx[i + (j + kk) * smat.nop] = deter;
  1346.                                 smat.ssx[i + kk + (j + kk) * smat.nop] = -deter;
  1347.                                 smat.ssx[i + kk + j * smat.nop] = deter;
  1348.                         }
  1349.                 }
  1350.                 kk *= 2;
  1351.  
  1352.         } /* for( k = 0; k < ka; k++ ) */
  1353.  
  1354.         return;
  1355. }
  1356.  
  1357. /*-----------------------------------------------------------------------*/
  1358.  
  1359. /* subroutine to solve the matrix equation lu*x=b where l is a unit */
  1360. /* lower triangular matrix and u is an upper triangular matrix both */
  1361. /* of which are stored in a.  the rhs vector b is input and the */
  1362. /* solution is returned through vector b.   (matrix transposed. */
  1363. void solve (int n, complex double *a, int *ip, complex double *b, int ndim)
  1364. {
  1365.         int i, ip1, j, k, pia;
  1366.         complex double sum, *scm = NULL;
  1367.  
  1368.         /* Allocate to scratch memory */
  1369.         mem_alloc ((void *) &scm, data.np2m * sizeof (complex double));
  1370.  
  1371.         /* forward substitution */
  1372.         for (i = 0; i < n; i++)
  1373.         {
  1374.                 pia = ip[i] - 1;
  1375.                 scm[i] = b[pia];
  1376.                 b[pia] = b[i];
  1377.                 ip1 = i + 1;
  1378.  
  1379.                 if (ip1 < n)
  1380.                         for (j = ip1; j < n; j++)
  1381.                                 b[j] -= a[j + i * ndim] * scm[i];
  1382.         }
  1383.  
  1384.         /* backward substitution */
  1385.         for (k = 0; k < n; k++)
  1386.         {
  1387.                 i = n - k - 1;
  1388.                 sum = CPLX_00;
  1389.                 ip1 = i + 1;
  1390.  
  1391.                 if (ip1 < n)
  1392.                         for (j = ip1; j < n; j++)
  1393.                                 sum += a[i + j * ndim] * b[j];
  1394.  
  1395.                 b[i] = (scm[i] - sum) / a[i + i * ndim];
  1396.         }
  1397.  
  1398.         free_ptr ((void *) &scm);
  1399.  
  1400.         return;
  1401. }
  1402.  
  1403. /*-----------------------------------------------------------------------*/
  1404.  
  1405. /* subroutine solves, for symmetric structures, handles the */
  1406. /* transformation of the right hand side vector and solution */
  1407. /* of the matrix eq. */
  1408. void solves (
  1409.     complex double *a,
  1410.     int *ip,
  1411.     complex double *b,
  1412.     int neq,
  1413.     int nrh,
  1414.     int np,
  1415.     int n,
  1416.     int mp,
  1417.     int m)
  1418. {
  1419.         int npeq, nrow, ic, i, kk, ia, ib, j, k;
  1420.         double fnop, fnorm;
  1421.         complex double sum, *scm = NULL;
  1422.  
  1423.         npeq = np + 2 * mp;
  1424.         smat.nop = neq / npeq;
  1425.         fnop = smat.nop;
  1426.         fnorm = 1. / fnop;
  1427.         nrow = neq;
  1428.  
  1429.         /* Allocate to scratch memory */
  1430.         mem_alloc ((void *) &scm, data.np2m * sizeof (complex double));
  1431.  
  1432.         if (smat.nop != 1)
  1433.         {
  1434.                 for (ic = 0; ic < nrh; ic++)
  1435.                 {
  1436.                         if ((n != 0) && (m != 0))
  1437.                         {
  1438.                                 for (i = 0; i < neq; i++)
  1439.                                         scm[i] = b[i + ic * neq];
  1440.  
  1441.                                 kk = 2 * mp;
  1442.                                 ia = np - 1;
  1443.                                 ib = n - 1;
  1444.                                 j = np - 1;
  1445.  
  1446.                                 for (k = 0; k < smat.nop; k++)
  1447.                                 {
  1448.                                         if (k != 0)
  1449.                                         {
  1450.                                                 for (i = 0; i < np; i++)
  1451.                                                 {
  1452.                                                         ia++;
  1453.                                                         j++;
  1454.                                                         b[j + ic * neq] = scm[ia];
  1455.                                                 }
  1456.  
  1457.                                                 if (k == (smat.nop - 1))
  1458.                                                         continue;
  1459.  
  1460.                                         } /* if( k != 0 ) */
  1461.  
  1462.                                         for (i = 0; i < kk; i++)
  1463.                                         {
  1464.                                                 ib++;
  1465.                                                 j++;
  1466.                                                 b[j + ic * neq] = scm[ib];
  1467.                                         }
  1468.  
  1469.                                 } /* for( k = 0; k < smat.nop; k++ ) */
  1470.  
  1471.                         } /* if( (n != 0) && (m != 0) ) */
  1472.  
  1473.                         /* transform matrix eq. rhs vector according to symmetry modes */
  1474.                         for (i = 0; i < npeq; i++)
  1475.                         {
  1476.                                 for (k = 0; k < smat.nop; k++)
  1477.                                 {
  1478.                                         ia = i + k * npeq;
  1479.                                         scm[k] = b[ia + ic * neq];
  1480.                                 }
  1481.  
  1482.                                 sum = scm[0];
  1483.                                 for (k = 1; k < smat.nop; k++)
  1484.                                         sum += scm[k];
  1485.  
  1486.                                 b[i + ic * neq] = sum * fnorm;
  1487.  
  1488.                                 for (k = 1; k < smat.nop; k++)
  1489.                                 {
  1490.                                         ia = i + k * npeq;
  1491.                                         sum = scm[0];
  1492.  
  1493.                                         for (j = 1; j < smat.nop; j++)
  1494.                                                 sum += scm[j] *
  1495.                                                        conjl (smat.ssx[k + j * smat.nop]);
  1496.  
  1497.                                         b[ia + ic * neq] = sum * fnorm;
  1498.                                 }
  1499.  
  1500.                         } /* for( i = 0; i < npeq; i++ ) */
  1501.  
  1502.                 } /* for( ic = 0; ic < nrh; ic++ ) */
  1503.  
  1504.         } /* if( smat.nop != 1) */
  1505.  
  1506.         /* solve each mode equation */
  1507.         for (kk = 0; kk < smat.nop; kk++)
  1508.         {
  1509.                 ia = kk * npeq;
  1510.                 ib = ia;
  1511.  
  1512.                 for (ic = 0; ic < nrh; ic++)
  1513.                         solve (npeq, &a[ib], &ip[ia], &b[ia + ic * neq], nrow);
  1514.  
  1515.         } /* for( kk = 0; kk < smat.nop; kk++ ) */
  1516.  
  1517.         if (smat.nop == 1)
  1518.         {
  1519.                 free_ptr ((void *) &scm);
  1520.                 return;
  1521.         }
  1522.  
  1523.         /* inverse transform the mode solutions */
  1524.         for (ic = 0; ic < nrh; ic++)
  1525.         {
  1526.                 for (i = 0; i < npeq; i++)
  1527.                 {
  1528.                         for (k = 0; k < smat.nop; k++)
  1529.                         {
  1530.                                 ia = i + k * npeq;
  1531.                                 scm[k] = b[ia + ic * neq];
  1532.                         }
  1533.  
  1534.                         sum = scm[0];
  1535.                         for (k = 1; k < smat.nop; k++)
  1536.                                 sum += scm[k];
  1537.  
  1538.                         b[i + ic * neq] = sum;
  1539.                         for (k = 1; k < smat.nop; k++)
  1540.                         {
  1541.                                 ia = i + k * npeq;
  1542.                                 sum = scm[0];
  1543.  
  1544.                                 for (j = 1; j < smat.nop; j++)
  1545.                                         sum += scm[j] * smat.ssx[k + j * smat.nop];
  1546.  
  1547.                                 b[ia + ic * neq] = sum;
  1548.                         }
  1549.  
  1550.                 } /* for( i = 0; i < npeq; i++ ) */
  1551.  
  1552.                 if ((n == 0) || (m == 0))
  1553.                         continue;
  1554.  
  1555.                 for (i = 0; i < neq; i++)
  1556.                         scm[i] = b[i + ic * neq];
  1557.  
  1558.                 kk = 2 * mp;
  1559.                 ia = np - 1;
  1560.                 ib = n - 1;
  1561.                 j = np - 1;
  1562.  
  1563.                 for (k = 0; k < smat.nop; k++)
  1564.                 {
  1565.                         if (k != 0)
  1566.                         {
  1567.                                 for (i = 0; i < np; i++)
  1568.                                 {
  1569.                                         ia++;
  1570.                                         j++;
  1571.                                         b[ia + ic * neq] = scm[j];
  1572.                                 }
  1573.  
  1574.                                 if (k == smat.nop)
  1575.                                         continue;
  1576.  
  1577.                         } /* if( k != 0 ) */
  1578.  
  1579.                         for (i = 0; i < kk; i++)
  1580.                         {
  1581.                                 ib++;
  1582.                                 j++;
  1583.                                 b[ib + ic * neq] = scm[j];
  1584.                         }
  1585.  
  1586.                 } /* for( k = 0; k < smat.nop; k++ ) */
  1587.  
  1588.         } /* for( ic = 0; ic < nrh; ic++ ) */
  1589.  
  1590.         free_ptr ((void *) &scm);
  1591.  
  1592.         return;
  1593. }
  1594.  
  1595. /*-----------------------------------------------------------------------*/
  1596.