Subversion Repositories Nec2c

Rev

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

  1. /******* Translated to the C language by N. Kyriazis  20 Aug 2003 ******
  2.  
  3.  Program NEC(input,tape5=input,output,tape11,tape12,tape13,tape14,
  4.  tape15,tape16,tape20,tape21)
  5.  
  6.  Numerical Electromagnetics Code (NEC2)  developed at Lawrence
  7.  Livermore lab., Livermore, CA.  (contact G. Burke at 415-422-8414
  8.  for problems with the NEC code. For problems with the vax implem-
  9.  entation, contact J. Breakall at 415-422-8196 or E. Domning at 415
  10.  422-5936)
  11.  file created 4/11/80.
  12.  
  13.                 ***********Notice**********
  14.  This computer code material was prepared as an account of work
  15.  sponsored by the United States government.  Neither the United
  16.  States nor the United States Department Of Energy, nor any of
  17.  their employees, nor any of their contractors, subcontractors,
  18.  or their employees, makes any warranty, express or implied, or
  19.  assumes any legal liability or responsibility for the accuracy,
  20.  completeness or usefulness of any information, apparatus, product
  21.  or process disclosed, or represents that its use would not infringe
  22.  privately-owned rights.
  23.  
  24. *******************************************************************/
  25.  
  26. #include "nec2c.h"
  27.  
  28. /* common  /data/ */
  29. extern data_t data;
  30.  
  31. /* common  /segj/ */
  32. extern segj_t segj;
  33.  
  34. /* common  /vsorc/ */
  35. extern vsorc_t vsorc;
  36.  
  37. /* common  /dataj/ */
  38. extern dataj_t dataj;
  39.  
  40. /* common  /zload/ */
  41. extern zload_t zload;
  42.  
  43. /* pointers to input/output files */
  44. extern FILE *input_fp, *output_fp, *plot_fp;
  45.  
  46. /*-------------------------------------------------------------------*/
  47.  
  48. /* fill incident field array for charge discontinuity voltage source */
  49. void qdsrc (int is, complex double v, complex double *e)
  50. {
  51.         int i, jx, j, jp1, ipr, ij, i1;
  52.         double xi, yi, zi, ai, cabi, sabi, salpi, tx, ty, tz;
  53.         complex double curd, etk, ets, etc;
  54.  
  55.         is--;
  56.         i = data.icon1[is];
  57.         data.icon1[is] = 0;
  58.         tbf (is + 1, 0);
  59.         data.icon1[is] = i;
  60.         dataj.s = data.si[is] * .5;
  61.         curd = CCJ * v /
  62.                ((logl (2. * dataj.s / data.bi[is]) - 1.) *
  63.                 (segj.bx[segj.jsno - 1] * cosl (TP * dataj.s) +
  64.                  segj.cx[segj.jsno - 1] * sinl (TP * dataj.s)) *
  65.                 data.wlam);
  66.         vsorc.vqds[vsorc.nqds] = v;
  67.         vsorc.iqds[vsorc.nqds] = is + 1;
  68.         vsorc.nqds++;
  69.  
  70.         for (jx = 0; jx < segj.jsno; jx++)
  71.         {
  72.                 j = segj.jco[jx] - 1;
  73.                 jp1 = j + 1;
  74.                 dataj.s = data.si[j];
  75.                 dataj.b = data.bi[j];
  76.                 dataj.xj = data.x[j];
  77.                 dataj.yj = data.y[j];
  78.                 dataj.zj = data.z[j];
  79.                 dataj.cabj = data.cab[j];
  80.                 dataj.sabj = data.sab[j];
  81.                 dataj.salpj = data.salp[j];
  82.  
  83.                 if (dataj.iexk != 0)
  84.                 {
  85.                         ipr = data.icon1[j];
  86.  
  87.                         if (ipr > PCHCON)
  88.                                 dataj.ind1 = 2;
  89.                         else if (ipr < 0)
  90.                         {
  91.                                 ipr = -ipr;
  92.                                 ipr--;
  93.                                 if (-data.icon1[ipr - 1] != jp1)
  94.                                         dataj.ind1 = 2;
  95.                                 else
  96.                                 {
  97.                                         xi = fabsl (
  98.                                             dataj.cabj * data.cab[ipr] +
  99.                                             dataj.sabj * data.sab[ipr] +
  100.                                             dataj.salpj * data.salp[ipr]);
  101.                                         if ((xi < 0.999999) ||
  102.                                             (fabsl (data.bi[ipr] / dataj.b - 1.) > 1.0e-6))
  103.                                                 dataj.ind1 = 2;
  104.                                         else
  105.                                                 dataj.ind1 = 0;
  106.                                 }
  107.                         } /* if( ipr < 0 ) */
  108.                         else if (ipr == 0)
  109.                                 dataj.ind1 = 1;
  110.                         else /* ipr > 0 */
  111.                         {
  112.                                 ipr--;
  113.                                 if (ipr != j)
  114.                                 {
  115.                                         if (data.icon2[ipr] != jp1)
  116.                                                 dataj.ind1 = 2;
  117.                                         else
  118.                                         {
  119.                                                 xi = fabsl (
  120.                                                     dataj.cabj * data.cab[ipr] +
  121.                                                     dataj.sabj * data.sab[ipr] +
  122.                                                     dataj.salpj * data.salp[ipr]);
  123.                                                 if ((xi < 0.999999) ||
  124.                                                     (fabsl (data.bi[ipr] / dataj.b - 1.) >
  125.                                                      1.0e-6))
  126.                                                         dataj.ind1 = 2;
  127.                                                 else
  128.                                                         dataj.ind1 = 0;
  129.                                         }
  130.                                 } /* if( ipr != j ) */
  131.                                 else
  132.                                 {
  133.                                         if (dataj.cabj * dataj.cabj + dataj.sabj * dataj.sabj >
  134.                                             1.0e-8)
  135.                                                 dataj.ind1 = 2;
  136.                                         else
  137.                                                 dataj.ind1 = 0;
  138.                                 }
  139.                         } /* else */
  140.  
  141.                         ipr = data.icon2[j];
  142.                         if (ipr > PCHCON)
  143.                                 dataj.ind2 = 2;
  144.                         else if (ipr < 0)
  145.                         {
  146.                                 ipr = -ipr;
  147.                                 ipr--;
  148.                                 if (-data.icon2[ipr] != jp1)
  149.                                         dataj.ind1 = 2;
  150.                                 else
  151.                                 {
  152.                                         xi = fabsl (
  153.                                             dataj.cabj * data.cab[ipr] +
  154.                                             dataj.sabj * data.sab[ipr] +
  155.                                             dataj.salpj * data.salp[ipr]);
  156.                                         if ((xi < 0.999999) ||
  157.                                             (fabsl (data.bi[ipr] / dataj.b - 1.) > 1.0e-6))
  158.                                                 dataj.ind1 = 2;
  159.                                         else
  160.                                                 dataj.ind1 = 0;
  161.                                 }
  162.                         } /* if( ipr < 0 ) */
  163.                         else if (ipr == 0)
  164.                                 dataj.ind2 = 1;
  165.                         else /* ipr > 0 */
  166.                         {
  167.                                 ipr--;
  168.                                 if (ipr != j)
  169.                                 {
  170.                                         if (data.icon1[ipr] != jp1)
  171.                                                 dataj.ind2 = 2;
  172.                                         else
  173.                                         {
  174.                                                 xi = fabsl (
  175.                                                     dataj.cabj * data.cab[ipr] +
  176.                                                     dataj.sabj * data.sab[ipr] +
  177.                                                     dataj.salpj * data.salp[ipr]);
  178.                                                 if ((xi < 0.999999) ||
  179.                                                     (fabsl (data.bi[ipr] / dataj.b - 1.) >
  180.                                                      1.0e-6))
  181.                                                         dataj.ind2 = 2;
  182.                                                 else
  183.                                                         dataj.ind2 = 0;
  184.                                         }
  185.                                 } /* if( ipr != j )*/
  186.                                 else
  187.                                 {
  188.                                         if (dataj.cabj * dataj.cabj + dataj.sabj * dataj.sabj >
  189.                                             1.0e-8)
  190.                                                 dataj.ind1 = 2;
  191.                                         else
  192.                                                 dataj.ind1 = 0;
  193.                                 }
  194.                         } /* else */
  195.  
  196.                 } /* if( dataj.iexk != 0) */
  197.  
  198.                 for (i = 0; i < data.n; i++)
  199.                 {
  200.                         ij = i - j;
  201.                         xi = data.x[i];
  202.                         yi = data.y[i];
  203.                         zi = data.z[i];
  204.                         ai = data.bi[i];
  205.                         efld (xi, yi, zi, ai, ij);
  206.                         cabi = data.cab[i];
  207.                         sabi = data.sab[i];
  208.                         salpi = data.salp[i];
  209.                         etk = dataj.exk * cabi + dataj.eyk * sabi + dataj.ezk * salpi;
  210.                         ets = dataj.exs * cabi + dataj.eys * sabi + dataj.ezs * salpi;
  211.                         etc = dataj.exc * cabi + dataj.eyc * sabi + dataj.ezc * salpi;
  212.                         e[i] =
  213.                             e[i] -
  214.                             (etk * segj.ax[jx] + ets * segj.bx[jx] + etc * segj.cx[jx]) * curd;
  215.                 }
  216.  
  217.                 if (data.m != 0)
  218.                 {
  219.                         i1 = data.n - 1;
  220.                         for (i = 0; i < data.m; i++)
  221.                         {
  222.                                 xi = data.px[i];
  223.                                 yi = data.py[i];
  224.                                 zi = data.pz[i];
  225.                                 hsfld (xi, yi, zi, 0.);
  226.                                 i1++;
  227.                                 tx = data.t2x[i];
  228.                                 ty = data.t2y[i];
  229.                                 tz = data.t2z[i];
  230.                                 etk = dataj.exk * tx + dataj.eyk * ty + dataj.ezk * tz;
  231.                                 ets = dataj.exs * tx + dataj.eys * ty + dataj.ezs * tz;
  232.                                 etc = dataj.exc * tx + dataj.eyc * ty + dataj.ezc * tz;
  233.                                 e[i1] += (etk * segj.ax[jx] + ets * segj.bx[jx] +
  234.                                           etc * segj.cx[jx]) *
  235.                                          curd * data.psalp[i];
  236.                                 i1++;
  237.                                 tx = data.t1x[i];
  238.                                 ty = data.t1y[i];
  239.                                 tz = data.t1z[i];
  240.                                 etk = dataj.exk * tx + dataj.eyk * ty + dataj.ezk * tz;
  241.                                 ets = dataj.exs * tx + dataj.eys * ty + dataj.ezs * tz;
  242.                                 etc = dataj.exc * tx + dataj.eyc * ty + dataj.ezc * tz;
  243.                                 e[i1] += (etk * segj.ax[jx] + ets * segj.bx[jx] +
  244.                                           etc * segj.cx[jx]) *
  245.                                          curd * data.psalp[i];
  246.                         }
  247.  
  248.                 } /* if( m != 0) */
  249.  
  250.                 if (zload.nload > 0)
  251.                         e[j] += zload.zarray[j] * curd * (segj.ax[jx] + segj.cx[jx]);
  252.  
  253.         } /* for( jx = 0; jx < segj.jsno; jx++ ) */
  254.  
  255.         return;
  256. }
  257.  
  258. /*-----------------------------------------------------------------------*/
  259.  
  260. void readmn (
  261.     char *gm,
  262.     int *i1,
  263.     int *i2,
  264.     int *i3,
  265.     int *i4,
  266.     double *f1,
  267.     double *f2,
  268.     double *f3,
  269.     double *f4,
  270.     double *f5,
  271.     double *f6)
  272. {
  273.         char line_buf[134];
  274.         int nlin, i, line_idx;
  275.         int nint = 4, nflt = 6;
  276.         int iarr[4] = {0, 0, 0, 0};
  277.         double rarr[6] = {0., 0., 0., 0., 0., 0.};
  278.  
  279.         /* read a line from input file */
  280.         if (load_line (line_buf, input_fp) == EOF)
  281.         {
  282.                 gm[0] = 'E';
  283.                 gm[1] = 'N';
  284.                 nlin = 2;
  285.         }
  286.         else
  287.         {
  288.                 /* get line length */
  289.                 nlin = strlen (line_buf);
  290.  
  291.                 /* abort if card's mnemonic too short or missing */
  292.                 if (nlin < 2)
  293.                 {
  294.                         fprintf (
  295.                             output_fp,
  296.                             "\n  COMMAND DATA CARD ERROR:"
  297.                             "\n  CARD'S MNEMONIC CODE TOO SHORT OR MISSING.");
  298.                         stop (-1);
  299.                 }
  300.  
  301.                 /* extract card's mnemonic code */
  302.                 gm[0] = line_buf[0];
  303.                 gm[1] = line_buf[1];
  304.                 /* was strncpy (gm, line_buf, 2); */
  305.         }
  306.         gm[2] = '\0';
  307.  
  308.         /* Exit if "XT" command read (for testing) */
  309.         if (strcmp (gm, "XT") == 0)
  310.         {
  311.                 fprintf (stderr, "\nnec2c: Exiting after an \"XT\" command in readgm()\n");
  312.                 fprintf (
  313.                     output_fp, "\n\n  nec2c: Exiting after an \"XT\" command in readgm()");
  314.                 stop (0);
  315.         }
  316.  
  317.         /* Return if only mnemonic on card */
  318.         if (nlin == 2)
  319.         {
  320.                 *i1 = *i2 = *i3 = *i4 = 0;
  321.                 *f1 = *f2 = *f3 = *f4 = *f5 = *f6 = 0.0;
  322.                 return;
  323.         }
  324.  
  325.         /* read integers from line */
  326.         line_idx = 1;
  327.         for (i = 0; i < nint; i++)
  328.         {
  329.                 /* Find first numerical character */
  330.                 while (((line_buf[++line_idx] < '0') || (line_buf[line_idx] > '9')) &&
  331.                        (line_buf[line_idx] != '+') && (line_buf[line_idx] != '-'))
  332.                         if ((line_buf[line_idx] == '\0'))
  333.                         {
  334.                                 *i1 = iarr[0];
  335.                                 *i2 = iarr[1];
  336.                                 *i3 = iarr[2];
  337.                                 *i4 = iarr[3];
  338.                                 *f1 = rarr[0];
  339.                                 *f2 = rarr[1];
  340.                                 *f3 = rarr[2];
  341.                                 *f4 = rarr[3];
  342.                                 *f5 = rarr[4];
  343.                                 *f6 = rarr[5];
  344.                                 return;
  345.                         }
  346.  
  347.                 /* read an integer from line */
  348.                 iarr[i] = atoi (&line_buf[line_idx]);
  349.  
  350.                 /* traverse numerical field to next ' ' or ',' or '\0' */
  351.                 line_idx--;
  352.                 while ((line_buf[++line_idx] != ' ') && (line_buf[line_idx] != '        ') &&
  353.                        (line_buf[line_idx] != ',') && (line_buf[line_idx] != '\0'))
  354.                 {
  355.                         /* test for non-numerical characters */
  356.                         if (((line_buf[line_idx] < '0') || (line_buf[line_idx] > '9')) &&
  357.                             (line_buf[line_idx] != '+') && (line_buf[line_idx] != '-'))
  358.                         {
  359.                                 fprintf (
  360.                                     output_fp,
  361.                                     "\n  COMMAND DATA CARD \"%s\" ERROR:"
  362.                                     "\n  NON-NUMERICAL CHARACTER '%c' IN INTEGER FIELD AT "
  363.                                     "CHAR. %d\n",
  364.                                     gm,
  365.                                     line_buf[line_idx],
  366.                                     (line_idx + 1));
  367.                                 stop (-1);
  368.                         }
  369.  
  370.                 } /* while( (line_buff[++line_idx] ... */
  371.  
  372.                 /* Return on end of line */
  373.                 if (line_buf[line_idx] == '\0')
  374.                 {
  375.                         *i1 = iarr[0];
  376.                         *i2 = iarr[1];
  377.                         *i3 = iarr[2];
  378.                         *i4 = iarr[3];
  379.                         *f1 = rarr[0];
  380.                         *f2 = rarr[1];
  381.                         *f3 = rarr[2];
  382.                         *f4 = rarr[3];
  383.                         *f5 = rarr[4];
  384.                         *f6 = rarr[5];
  385.                         return;
  386.                 }
  387.  
  388.         } /* for( i = 0; i < nint; i++ ) */
  389.  
  390.         /* read doubles from line */
  391.         for (i = 0; i < nflt; i++)
  392.         {
  393.                 /* Find first numerical character */
  394.                 while (((line_buf[++line_idx] < '0') || (line_buf[line_idx] > '9')) &&
  395.                        (line_buf[line_idx] != '+') && (line_buf[line_idx] != '-') &&
  396.                        (line_buf[line_idx] != '.'))
  397.                         if ((line_buf[line_idx] == '\0'))
  398.                         {
  399.                                 *i1 = iarr[0];
  400.                                 *i2 = iarr[1];
  401.                                 *i3 = iarr[2];
  402.                                 *i4 = iarr[3];
  403.                                 *f1 = rarr[0];
  404.                                 *f2 = rarr[1];
  405.                                 *f3 = rarr[2];
  406.                                 *f4 = rarr[3];
  407.                                 *f5 = rarr[4];
  408.                                 *f6 = rarr[5];
  409.                                 return;
  410.                         }
  411.  
  412.                 /* read a double from line */
  413.                 rarr[i] = atof (&line_buf[line_idx]);
  414.  
  415.                 /* traverse numerical field to next ' ' or ',' */
  416.                 line_idx--;
  417.                 while ((line_buf[++line_idx] != ' ') && (line_buf[line_idx] != '        ') &&
  418.                        (line_buf[line_idx] != ',') && (line_buf[line_idx] != '\0'))
  419.                 {
  420.                         /* test for non-numerical characters */
  421.                         if (((line_buf[line_idx] < '0') || (line_buf[line_idx] > '9')) &&
  422.                             (line_buf[line_idx] != '.') && (line_buf[line_idx] != '+') &&
  423.                             (line_buf[line_idx] != '-') && (line_buf[line_idx] != 'E') &&
  424.                             (line_buf[line_idx] != 'e'))
  425.                         {
  426.                                 fprintf (
  427.                                     output_fp,
  428.                                     "\n  COMMAND DATA CARD \"%s\" ERROR:"
  429.                                     "\n  NON-NUMERICAL CHARACTER '%c' IN FLOAT FIELD AT CHAR. "
  430.                                     "%d\n",
  431.                                     gm,
  432.                                     line_buf[line_idx],
  433.                                     (line_idx + 1));
  434.                                 stop (-1);
  435.                         }
  436.  
  437.                 } /* while( (line_buff[++line_idx] ... */
  438.  
  439.                 /* Return on end of line */
  440.                 if (line_buf[line_idx] == '\0')
  441.                 {
  442.                         *i1 = iarr[0];
  443.                         *i2 = iarr[1];
  444.                         *i3 = iarr[2];
  445.                         *i4 = iarr[3];
  446.                         *f1 = rarr[0];
  447.                         *f2 = rarr[1];
  448.                         *f3 = rarr[2];
  449.                         *f4 = rarr[3];
  450.                         *f5 = rarr[4];
  451.                         *f6 = rarr[5];
  452.                         return;
  453.                 }
  454.  
  455.         } /* for( i = 0; i < nflt; i++ ) */
  456.  
  457.         *i1 = iarr[0];
  458.         *i2 = iarr[1];
  459.         *i3 = iarr[2];
  460.         *i4 = iarr[3];
  461.         *f1 = rarr[0];
  462.         *f2 = rarr[1];
  463.         *f3 = rarr[2];
  464.         *f4 = rarr[3];
  465.         *f5 = rarr[4];
  466.         *f6 = rarr[5];
  467.  
  468.         return;
  469. }
  470.  
  471. /*-----------------------------------------------------------------------*/
  472.