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. /* pointers to input/output files */
  35. extern FILE *input_fp, *output_fp, *plot_fp;
  36.  
  37. /* common  /plot/ */
  38. extern plot_t plot;
  39.  
  40. /*-------------------------------------------------------------------*/
  41.  
  42. /* arc generates segment geometry data for an arc of ns segments */
  43. void arc (int itg, int ns, double rada, double ang1, double ang2, double rad)
  44. {
  45.         int ist, i, mreq;
  46.         double ang, dang, xs1, xs2, zs1, zs2;
  47.  
  48.         ist = data.n;
  49.         data.n += ns;
  50.         data.np = data.n;
  51.         data.mp = data.m;
  52.         data.ipsym = 0;
  53.  
  54.         if (ns < 1)
  55.                 return;
  56.  
  57.         if (fabsl (ang2 - ang1) < 360.00001)
  58.         {
  59.                 /* Reallocate tags buffer */
  60.                 mem_realloc ((void *) &data.itag, (data.n + data.m) * sizeof (int));
  61.  
  62.                 /* Reallocate wire buffers */
  63.                 mreq = data.n * sizeof (double);
  64.                 mem_realloc ((void *) &data.x1, mreq);
  65.                 mem_realloc ((void *) &data.y1, mreq);
  66.                 mem_realloc ((void *) &data.z1, mreq);
  67.                 mem_realloc ((void *) &data.x2, mreq);
  68.                 mem_realloc ((void *) &data.y2, mreq);
  69.                 mem_realloc ((void *) &data.z2, mreq);
  70.                 mem_realloc ((void *) &data.bi, mreq);
  71.  
  72.                 ang = ang1 * TA;
  73.                 dang = (ang2 - ang1) * TA / ns;
  74.                 xs1 = rada * cosl (ang);
  75.                 zs1 = rada * sinl (ang);
  76.  
  77.                 for (i = ist; i < data.n; i++)
  78.                 {
  79.                         ang += dang;
  80.                         xs2 = rada * cosl (ang);
  81.                         zs2 = rada * sinl (ang);
  82.                         data.x1[i] = xs1;
  83.  
  84.                         data.y1[i] = 0.;
  85.                         data.z1[i] = zs1;
  86.                         data.x2[i] = xs2;
  87.                         data.y2[i] = 0.;
  88.                         data.z2[i] = zs2;
  89.                         xs1 = xs2;
  90.                         zs1 = zs2;
  91.                         data.bi[i] = rad;
  92.                         data.itag[i] = itg;
  93.  
  94.                 } /* for( i = ist; i < data.n; i++ ) */
  95.  
  96.         } /* if( fabsl( ang2- ang1) < 360.00001) */
  97.         else
  98.         {
  99.                 fprintf (output_fp, "\n  ERROR -- ARC ANGLE EXCEEDS 360 DEGREES");
  100.                 stop (-1);
  101.         }
  102.  
  103.         return;
  104. }
  105.  
  106. /*-----------------------------------------------------------------------*/
  107.  
  108. /* connect sets up segment connection data in arrays icon1 and */
  109. /* icon2 by searching for segment ends that are in contact. */
  110. void conect (int ignd)
  111. {
  112.         int i, iz, ic, j, jx, ix, ixx, iseg, iend, jend, jump, ipf;
  113.         // int nsflg;
  114.         double sep = 0., xi1, yi1, zi1, xi2, yi2, zi2;
  115.         double slen, xa, ya, za, xs, ys, zs;
  116.  
  117.         segj.maxcon = 1;
  118.  
  119.         if (ignd != 0)
  120.         {
  121.                 fprintf (output_fp, "\n\n   GROUND PLANE SPECIFIED.\n");
  122.  
  123.                 if (ignd > 0)
  124.                         fprintf (
  125.                             output_fp,
  126.                             "\n     WHERE WIRE ENDS TOUCH GROUND, CURRENT WILL"
  127.                             " BE INTERPOLATED TO IMAGE IN GROUND PLANE.\n");
  128.  
  129.                 if (data.ipsym == 2)
  130.                 {
  131.                         data.np = 2 * data.np;
  132.                         data.mp = 2 * data.mp;
  133.                 }
  134.  
  135.                 if (abs (data.ipsym) > 2)
  136.                 {
  137.                         data.np = data.n;
  138.                         data.mp = data.m;
  139.                 }
  140.  
  141.                 /*** possibly should be error condition?? **/
  142.                 if (data.np > data.n)
  143.                 {
  144.                         fprintf (output_fp, "\n ERROR: NP > N IN CONECT()");
  145.                         stop (-1);
  146.                 }
  147.  
  148.                 if ((data.np == data.n) && (data.mp == data.m))
  149.                         data.ipsym = 0;
  150.  
  151.         } /* if( ignd != 0) */
  152.  
  153.         if (data.n != 0)
  154.         {
  155.                 /* Allocate memory to connections */
  156.                 mem_realloc ((void *) &data.icon1, (data.n + data.m) * sizeof (int));
  157.                 mem_realloc ((void *) &data.icon2, (data.n + data.m) * sizeof (int));
  158.  
  159.                 for (i = 0; i < data.n; i++)
  160.                 {
  161.                         data.icon1[i] = data.icon2[i] = 0;
  162.                         iz = i + 1;
  163.                         xi1 = data.x1[i];
  164.                         yi1 = data.y1[i];
  165.                         zi1 = data.z1[i];
  166.                         xi2 = data.x2[i];
  167.                         yi2 = data.y2[i];
  168.                         zi2 = data.z2[i];
  169.                         slen = sqrtl (
  170.                                    (xi2 - xi1) * (xi2 - xi1) + (yi2 - yi1) * (yi2 - yi1) +
  171.                                    (zi2 - zi1) * (zi2 - zi1)) *
  172.                                SMIN;
  173.  
  174.                         /* determine connection data for end 1 of segment. */
  175.                         jump = FALSE;
  176.                         if (ignd > 0)
  177.                         {
  178.                                 if (zi1 <= -slen)
  179.                                 {
  180.                                         fprintf (
  181.                                             output_fp,
  182.                                             "\n  GEOMETRY DATA ERROR -- SEGMENT"
  183.                                             " %d EXTENDS BELOW GROUND",
  184.                                             iz);
  185.                                         stop (-1);
  186.                                 }
  187.  
  188.                                 if (zi1 <= slen)
  189.                                 {
  190.                                         data.icon1[i] = iz;
  191.                                         data.z1[i] = 0.;
  192.                                         jump = TRUE;
  193.  
  194.                                 } /* if( zi1 <= slen) */
  195.  
  196.                         } /* if( ignd > 0) */
  197.  
  198.                         if (!jump)
  199.                         {
  200.                                 ic = i;
  201.                                 for (j = 1; j < data.n; j++)
  202.                                 {
  203.                                         ic++;
  204.                                         if (ic >= data.n)
  205.                                                 ic = 0;
  206.  
  207.                                         sep = fabsl (xi1 - data.x1[ic]) +
  208.                                               fabsl (yi1 - data.y1[ic]) +
  209.                                               fabsl (zi1 - data.z1[ic]);
  210.                                         if (sep <= slen)
  211.                                         {
  212.                                                 data.icon1[i] = -(ic + 1);
  213.                                                 break;
  214.                                         }
  215.  
  216.                                         sep = fabsl (xi1 - data.x2[ic]) +
  217.                                               fabsl (yi1 - data.y2[ic]) +
  218.                                               fabsl (zi1 - data.z2[ic]);
  219.                                         if (sep <= slen)
  220.                                         {
  221.                                                 data.icon1[i] = (ic + 1);
  222.                                                 break;
  223.                                         }
  224.  
  225.                                 } /* for( j = 1; j < data.n; j++) */
  226.  
  227.                         } /* if( ! jump ) */
  228.  
  229.                         /* determine connection data for end 2 of segment. */
  230.                         if ((ignd > 0) || jump)
  231.                         {
  232.                                 if (zi2 <= -slen)
  233.                                 {
  234.                                         fprintf (
  235.                                             output_fp,
  236.                                             "\n  GEOMETRY DATA ERROR -- SEGMENT"
  237.                                             " %d EXTENDS BELOW GROUND",
  238.                                             iz);
  239.                                         stop (-1);
  240.                                 }
  241.  
  242.                                 if (zi2 <= slen)
  243.                                 {
  244.                                         if (data.icon1[i] == iz)
  245.                                         {
  246.                                                 fprintf (
  247.                                                     output_fp,
  248.                                                     "\n  GEOMETRY DATA ERROR -- SEGMENT"
  249.                                                     " %d LIES IN GROUND PLANE",
  250.                                                     iz);
  251.                                                 stop (-1);
  252.                                         }
  253.  
  254.                                         data.icon2[i] = iz;
  255.                                         data.z2[i] = 0.;
  256.                                         continue;
  257.  
  258.                                 } /* if( zi2 <= slen) */
  259.  
  260.                         } /* if( ignd > 0) */
  261.  
  262.                         ic = i;
  263.                         for (j = 1; j < data.n; j++)
  264.                         {
  265.                                 ic++;
  266.                                 if (ic >= data.n)
  267.                                         ic = 0;
  268.  
  269.                                 sep = fabsl (xi2 - data.x1[ic]) + fabsl (yi2 - data.y1[ic]) +
  270.                                       fabsl (zi2 - data.z1[ic]);
  271.                                 if (sep <= slen)
  272.                                 {
  273.                                         data.icon2[i] = (ic + 1);
  274.                                         break;
  275.                                 }
  276.  
  277.                                 sep = fabsl (xi2 - data.x2[ic]) + fabsl (yi2 - data.y2[ic]) +
  278.                                       fabsl (zi2 - data.z2[ic]);
  279.                                 if (sep <= slen)
  280.                                 {
  281.                                         data.icon2[i] = -(ic + 1);
  282.                                         break;
  283.                                 }
  284.  
  285.                         } /* for( j = 1; j < data.n; j++ ) */
  286.  
  287.                 } /* for( i = 0; i < data.n; i++ ) */
  288.  
  289.                 /* find wire-surface connections for new patches */
  290.                 if (data.m != 0)
  291.                 {
  292.                         ix = -1;
  293.                         i = 0;
  294.                         while (++i <= data.m)
  295.                         {
  296.                                 ix++;
  297.                                 xs = data.px[ix];
  298.                                 ys = data.py[ix];
  299.                                 zs = data.pz[ix];
  300.  
  301.                                 for (iseg = 0; iseg < data.n; iseg++)
  302.                                 {
  303.                                         xi1 = data.x1[iseg];
  304.                                         yi1 = data.y1[iseg];
  305.                                         zi1 = data.z1[iseg];
  306.                                         xi2 = data.x2[iseg];
  307.                                         yi2 = data.y2[iseg];
  308.                                         zi2 = data.z2[iseg];
  309.  
  310.                                         /* for first end of segment */
  311.                                         slen = (fabsl (xi2 - xi1) + fabsl (yi2 - yi1) +
  312.                                                 fabsl (zi2 - zi1)) *
  313.                                                SMIN;
  314.                                         sep = fabsl (xi1 - xs) + fabsl (yi1 - ys) +
  315.                                               fabsl (zi1 - zs);
  316.  
  317.                                         /* connection - divide patch into 4 patches at present
  318.                                          * array loc. */
  319.                                         if (sep <= slen)
  320.                                         {
  321.                                                 data.icon1[iseg] = PCHCON + i;
  322.                                                 ic = 0;
  323.                                                 subph (i, ic);
  324.                                                 break;
  325.                                         }
  326.  
  327.                                         sep = fabsl (xi2 - xs) + fabsl (yi2 - ys) +
  328.                                               fabsl (zi2 - zs);
  329.                                         if (sep <= slen)
  330.                                         {
  331.                                                 data.icon2[iseg] = PCHCON + i;
  332.                                                 ic = 0;
  333.                                                 subph (i, ic);
  334.                                                 break;
  335.                                         }
  336.  
  337.                                 } /* for( iseg = 0; iseg < data.n; iseg++ ) */
  338.  
  339.                         } /* while( ++i <= data.m ) */
  340.  
  341.                 } /* if( data.m != 0) */
  342.  
  343.         } /* if( data.n != 0) */
  344.  
  345.         fprintf (
  346.             output_fp,
  347.             "\n\n"
  348.             "   TOTAL SEGMENTS USED=%5d     "
  349.             "NO. SEG. IN A SYMMETRIC CELL=%5d"
  350.             "     SYMMETRY FLAG=  %d\n",
  351.             data.n,
  352.             data.np,
  353.             data.ipsym);
  354.  
  355.         if (data.m > 0)
  356.                 fprintf (
  357.                     output_fp,
  358.                     "\n"
  359.                     "       TOTAL PATCHES USED: %d   PATCHES"
  360.                     " IN A SYMMETRIC CELL: %d",
  361.                     data.m,
  362.                     data.mp);
  363.  
  364.         iseg = (data.n + data.m) / (data.np + data.mp);
  365.         if (iseg != 1)
  366.         {
  367.                 /*** may be error condition?? ***/
  368.                 if (data.ipsym == 0)
  369.                 {
  370.                         fprintf (output_fp, "\n  ERROR: IPSYM=0 IN CONECT()");
  371.                         stop (-1);
  372.                 }
  373.  
  374.                 if (data.ipsym < 0)
  375.                         fprintf (
  376.                             output_fp,
  377.                             "\n  STRUCTURE HAS %d FOLD ROTATIONAL SYMMETRY\n",
  378.                             iseg);
  379.                 else
  380.                 {
  381.                         ic = iseg / 2;
  382.                         if (iseg == 8)
  383.                                 ic = 3;
  384.                         fprintf (output_fp, "\n  STRUCTURE HAS %d PLANES OF SYMMETRY\n", ic);
  385.                 } /* if( data.ipsym < 0 ) */
  386.  
  387.         } /* if( iseg == 1) */
  388.  
  389.         if (data.n == 0)
  390.                 return;
  391.  
  392.         /* Allocate to connection buffers */
  393.         mem_realloc ((void *) &segj.jco, segj.maxcon * sizeof (int));
  394.  
  395.         /* adjust connected seg. ends to exactly coincide.  print junctions */
  396.         /* of 3 or more seg.  also find old seg. connecting to new seg. */
  397.         iseg = 0;
  398.         ipf = FALSE;
  399.         for (j = 0; j < data.n; j++)
  400.         {
  401.                 jx = j + 1;
  402.                 iend = -1;
  403.                 jend = -1;
  404.                 ix = data.icon1[j];
  405.                 ic = 1;
  406.                 segj.jco[0] = -jx;
  407.                 xa = data.x1[j];
  408.                 ya = data.y1[j];
  409.                 za = data.z1[j];
  410.  
  411.                 while (TRUE)
  412.                 {
  413.                         if ((ix != 0) && (ix != (j + 1)) && (ix <= PCHCON))
  414.                         {
  415.                                 // nsflg=0;
  416.  
  417.                                 do
  418.                                 {
  419.                                         if (ix == 0)
  420.                                         {
  421.                                                 fprintf (
  422.                                                     output_fp,
  423.                                                     "\n  CONNECT - SEGMENT CONNECTION ERROR "
  424.                                                     "FOR SEGMENT: %d",
  425.                                                     ix);
  426.                                                 stop (-1);
  427.                                         }
  428.  
  429.                                         if (ix < 0)
  430.                                                 ix = -ix;
  431.                                         else
  432.                                                 jend = -jend;
  433.  
  434.                                         jump = FALSE;
  435.  
  436.                                         if (ix == jx)
  437.                                                 break;
  438.  
  439.                                         if (ix < jx)
  440.                                         {
  441.                                                 jump = TRUE;
  442.                                                 break;
  443.                                         }
  444.  
  445.                                         /* Record max. no. of connections */
  446.                                         ic++;
  447.                                         if (ic >= segj.maxcon)
  448.                                         {
  449.                                                 segj.maxcon = ic + 1;
  450.                                                 mem_realloc (
  451.                                                     (void *) &segj.jco,
  452.                                                     segj.maxcon * sizeof (int));
  453.                                         }
  454.                                         segj.jco[ic - 1] = ix * jend;
  455.  
  456.                                         /*if( ix > 0){
  457.                                               //nsflg=1;
  458.                                         }*/
  459.                                         ixx = ix - 1;
  460.                                         if (jend != 1)
  461.                                         {
  462.                                                 xa = xa + data.x1[ixx];
  463.                                                 ya = ya + data.y1[ixx];
  464.                                                 za = za + data.z1[ixx];
  465.                                                 ix = data.icon1[ixx];
  466.                                                 continue;
  467.                                         }
  468.  
  469.                                         xa = xa + data.x2[ixx];
  470.                                         ya = ya + data.y2[ixx];
  471.                                         za = za + data.z2[ixx];
  472.                                         ix = data.icon2[ixx];
  473.  
  474.                                 } /* do */
  475.                                 while (ix != 0);
  476.  
  477.                                 if (jump && (iend == 1))
  478.                                         break;
  479.                                 else if (jump)
  480.                                 {
  481.                                         iend = 1;
  482.                                         jend = 1;
  483.                                         ix = data.icon2[j];
  484.                                         ic = 1;
  485.                                         segj.jco[0] = jx;
  486.                                         xa = data.x2[j];
  487.                                         ya = data.y2[j];
  488.                                         za = data.z2[j];
  489.                                         continue;
  490.                                 }
  491.  
  492.                                 sep = (double) ic;
  493.                                 xa = xa / sep;
  494.                                 ya = ya / sep;
  495.                                 za = za / sep;
  496.  
  497.                                 for (i = 0; i < ic; i++)
  498.                                 {
  499.                                         ix = segj.jco[i];
  500.                                         if (ix <= 0)
  501.                                         {
  502.                                                 ix = -ix;
  503.                                                 ixx = ix - 1;
  504.                                                 data.x1[ixx] = xa;
  505.                                                 data.y1[ixx] = ya;
  506.                                                 data.z1[ixx] = za;
  507.                                                 continue;
  508.                                         }
  509.  
  510.                                         ixx = ix - 1;
  511.                                         data.x2[ixx] = xa;
  512.                                         data.y2[ixx] = ya;
  513.                                         data.z2[ixx] = za;
  514.  
  515.                                 } /* for( i = 0; i < ic; i++ ) */
  516.  
  517.                                 if (!ipf)
  518.                                 {
  519.                                         fprintf (
  520.                                             output_fp,
  521.                                             "\n\n"
  522.                                             "         - MULTIPLE WIRE JUNCTIONS -\n"
  523.                                             " JUNCTION    SEGMENTS  (- FOR END 1, + FOR END "
  524.                                             "2)");
  525.  
  526.                                         ipf = TRUE;
  527.                                 }
  528.                                 if (ic >= 3)
  529.                                 {
  530.                                         iseg++;
  531.                                         fprintf (output_fp, "\n %5d     ", iseg);
  532.  
  533.                                         for (i = 1; i <= ic; i++)
  534.                                         {
  535.                                                 fprintf (output_fp, " %4d", segj.jco[i - 1]);
  536.                                                 if (!(i % 20))
  537.                                                         fprintf (
  538.                                                             output_fp, "\n              ");
  539.                                         }
  540.                                 } /* if( ic >= 3) */
  541.  
  542.                         } /*if( (ix != 0) && (ix != j) && (ix <= PCHCON) ) */
  543.  
  544.                         if (iend == 1)
  545.                                 break;
  546.  
  547.                         iend = 1;
  548.                         jend = 1;
  549.                         ix = data.icon2[j];
  550.                         ic = 1;
  551.                         segj.jco[0] = jx;
  552.                         xa = data.x2[j];
  553.                         ya = data.y2[j];
  554.                         za = data.z2[j];
  555.  
  556.                 } /* while( TRUE ) */
  557.  
  558.         } /* for( j = 0; j < data.n; j++ ) */
  559.           /* print NONE if nothing printed */
  560.         if (ipf == 1 && iseg == 0)
  561.         {
  562.                 fprintf (output_fp, "\n  NONE");
  563.         }
  564.  
  565.         mem_realloc ((void *) &segj.ax, segj.maxcon * sizeof (double));
  566.         mem_realloc ((void *) &segj.bx, segj.maxcon * sizeof (double));
  567.         mem_realloc ((void *) &segj.cx, segj.maxcon * sizeof (double));
  568.  
  569.         return;
  570. }
  571.  
  572. /*-----------------------------------------------------------------------*/
  573.  
  574. /* datagn is the main routine for input of geometry data. */
  575. void datagn (void)
  576. {
  577.         char gm[3];
  578.         char ifx[2] = {'*', 'X'}, ify[2] = {'*', 'Y'}, ifz[2] = {'*', 'Z'};
  579.         char ipt[4] = {'P', 'R', 'T', 'Q'};
  580.  
  581.         /* input card mnemonic list */
  582.         /* "XT" stands for "exit", added for testing */
  583. #define GM_NUM 12
  584.         char *atst[GM_NUM] = {
  585.             "GW", "GX", "GR", "GS", "GE", "GM", "SP", "SM", "GA", "SC", "GH", "GF"};
  586.  
  587.         int nwire, isct, iphd, i1, i2, itg, iy, iz, mreq;
  588.         int ix, i, ns, gm_num; /* geometry card id as a number */
  589.         double rad, xs1, xs2, ys1, ys2, zs1, zs2, x4 = 0, y4 = 0, z4 = 0;
  590.         double x3 = 0, y3 = 0, z3 = 0, xw1, xw2, yw1, yw2, zw1, zw2;
  591.         double dummy;
  592.  
  593.         data.ipsym = 0;
  594.         nwire = 0;
  595.         data.n = 0;
  596.         data.np = 0;
  597.         data.m = 0;
  598.         data.mp = 0;
  599.         isct = 0;
  600.         iphd = FALSE;
  601.  
  602.         /* read geometry data card and branch to */
  603.         /* section for operation requested */
  604.         do
  605.         {
  606.                 readgm (gm, &itg, &ns, &xw1, &yw1, &zw1, &xw2, &yw2, &zw2, &rad);
  607.  
  608.                 /* identify card id mnemonic */
  609.                 for (gm_num = 0; gm_num < GM_NUM; gm_num++)
  610.                         if (strncmp (gm, atst[gm_num], 2) == 0)
  611.                                 break;
  612.  
  613.                 if (iphd == FALSE)
  614.                 {
  615.                         fprintf (
  616.                             output_fp,
  617.                             "\n\n\n\n"
  618.                             "                                 "
  619.                             "- - - STRUCTURE SPECIFICATION - - -\n\n"
  620.                             "                                     "
  621.                             "COORDINATES MUST BE INPUT IN\n"
  622.                             "                                     "
  623.                             "METERS OR BE SCALED TO METERS\n"
  624.                             "                                     "
  625.                             "BEFORE STRUCTURE INPUT IS ENDED\n");
  626.  
  627.                         fprintf (
  628.                             output_fp,
  629.                             "\n\n"
  630.                             "  WIRE                                                          "
  631.                             "                     NO. OF    FIRST  LAST     TAG\n"
  632.                             "  NO.        X1         Y1         Z1          X2         Y2     "
  633.                             "    Z2      RADIUS   SEG.     SEG.   SEG.     NO.");
  634.  
  635.                         iphd = 1;
  636.                 }
  637.  
  638.                 if (gm_num != 10)
  639.                         isct = 0;
  640.  
  641.                 switch (gm_num)
  642.                 {
  643.                 case 0: /* "gw" card, generate segment data for straight wire. */
  644.  
  645.                         nwire++;
  646.                         i1 = data.n + 1;
  647.                         i2 = data.n + ns;
  648.  
  649.                         fprintf (
  650.                             output_fp,
  651.                             "\n"
  652.                             " %5d %10.5f %10.5f %10.5f  %10.5f"
  653.                             " %10.5f %10.5f %10.5f  %5d    %5d %5d   %5d",
  654.                             nwire,
  655.                             xw1,
  656.                             yw1,
  657.                             zw1,
  658.                             xw2,
  659.                             yw2,
  660.                             zw2,
  661.                             rad,
  662.                             ns,
  663.                             i1,
  664.                             i2,
  665.                             itg);
  666.  
  667.                         if (rad != 0)
  668.                         {
  669.                                 xs1 = 1.;
  670.                                 ys1 = 1.;
  671.                         }
  672.                         else
  673.                         {
  674.                                 readgm (
  675.                                     gm,
  676.                                     &ix,
  677.                                     &iy,
  678.                                     &xs1,
  679.                                     &ys1,
  680.                                     &zs1,
  681.                                     &dummy,
  682.                                     &dummy,
  683.                                     &dummy,
  684.                                     &dummy);
  685.  
  686.                                 if (strcmp (gm, "GC") != 0)
  687.                                 {
  688.                                         fprintf (output_fp, "\n  GEOMETRY DATA CARD ERROR");
  689.                                         stop (-1);
  690.                                 }
  691.  
  692.                                 fprintf (
  693.                                     output_fp,
  694.                                     "\n  ABOVE WIRE IS TAPERED.  SEGMENT LENGTH RATIO: %9.5f\n"
  695.                                     "                                 "
  696.                                     "RADIUS FROM: %9.5f TO: %9.5f",
  697.                                     xs1,
  698.                                     ys1,
  699.                                     zs1);
  700.  
  701.                                 if ((ys1 == 0) || (zs1 == 0))
  702.                                 {
  703.                                         fprintf (output_fp, "\n  GEOMETRY DATA CARD ERROR");
  704.                                         stop (-1);
  705.                                 }
  706.  
  707.                                 rad = ys1;
  708.                                 ys1 = powl ((zs1 / ys1), (1. / (ns - 1.)));
  709.                         }
  710.  
  711.                         wire (xw1, yw1, zw1, xw2, yw2, zw2, rad, xs1, ys1, ns, itg);
  712.  
  713.                         continue;
  714.  
  715.                         /* reflect structure along x,y, or z */
  716.                         /* axes or rotate to form cylinder.  */
  717.                 case 1: /* "gx" card */
  718.  
  719.                         iy = ns / 10;
  720.                         iz = ns - iy * 10;
  721.                         ix = iy / 10;
  722.                         iy = iy - ix * 10;
  723.  
  724.                         if (ix != 0)
  725.                                 ix = 1;
  726.                         if (iy != 0)
  727.                                 iy = 1;
  728.                         if (iz != 0)
  729.                                 iz = 1;
  730.  
  731.                         fprintf (
  732.                             output_fp,
  733.                             "\n  STRUCTURE REFLECTED ALONG THE AXES %c %c %c"
  734.                             " - TAGS INCREMENTED BY %d\n",
  735.                             ifx[ix],
  736.                             ify[iy],
  737.                             ifz[iz],
  738.                             itg);
  739.  
  740.                         reflc (ix, iy, iz, itg, ns);
  741.  
  742.                         continue;
  743.  
  744.                 case 2: /* "gr" card */
  745.  
  746.                         fprintf (
  747.                             output_fp,
  748.                             "\n  STRUCTURE ROTATED ABOUT Z-AXIS %d TIMES"
  749.                             " - LABELS INCREMENTED BY %d\n",
  750.                             ns,
  751.                             itg);
  752.  
  753.                         ix = -1;
  754.                         iz = 0;
  755.                         reflc (ix, iy, iz, itg, ns);
  756.  
  757.                         continue;
  758.  
  759.                 case 3: /* "gs" card, scale structure dimensions by factor xw1. */
  760.  
  761.                         if (data.n > 0)
  762.                         {
  763.                                 for (i = 0; i < data.n; i++)
  764.                                 {
  765.                                         data.x1[i] = data.x1[i] * xw1;
  766.                                         data.y1[i] = data.y1[i] * xw1;
  767.                                         data.z1[i] = data.z1[i] * xw1;
  768.                                         data.x2[i] = data.x2[i] * xw1;
  769.                                         data.y2[i] = data.y2[i] * xw1;
  770.                                         data.z2[i] = data.z2[i] * xw1;
  771.                                         data.bi[i] = data.bi[i] * xw1;
  772.                                 }
  773.                         } /* if( data.n >= n2) */
  774.  
  775.                         if (data.m > 0)
  776.                         {
  777.                                 yw1 = xw1 * xw1;
  778.                                 for (i = 0; i < data.m; i++)
  779.                                 {
  780.                                         data.px[i] = data.px[i] * xw1;
  781.                                         data.py[i] = data.py[i] * xw1;
  782.                                         data.pz[i] = data.pz[i] * xw1;
  783.                                         data.pbi[i] = data.pbi[i] * yw1;
  784.                                 }
  785.                         } /* if( data.m >= m2) */
  786.  
  787.                         fprintf (output_fp, "\n      STRUCTURE SCALED BY FACTOR %10.5f", xw1);
  788.  
  789.                         continue;
  790.  
  791.                 case 4: /* "ge" card, terminate structure geometry input. */
  792.  
  793.                         if (ns != 0)
  794.                         {
  795.                                 plot.iplp1 = 1;
  796.                                 plot.iplp2 = 1;
  797.                         }
  798.  
  799.                         conect (itg);
  800.  
  801.                         if (data.n != 0)
  802.                         {
  803.                                 /* Allocate wire buffers */
  804.                                 mreq = data.n * sizeof (double);
  805.                                 mem_realloc ((void *) &data.si, mreq);
  806.                                 mem_realloc ((void *) &data.sab, mreq);
  807.                                 mem_realloc ((void *) &data.cab, mreq);
  808.                                 mem_realloc ((void *) &data.salp, mreq);
  809.                                 mem_realloc ((void *) &data.x, mreq);
  810.                                 mem_realloc ((void *) &data.y, mreq);
  811.                                 mem_realloc ((void *) &data.z, mreq);
  812.  
  813.                                 fprintf (
  814.                                     output_fp,
  815.                                     "\n\n\n\n\n"
  816.                                     "                                "
  817.                                     " - - - - SEGMENTATION DATA - - - -\n\n"
  818.                                     "                                       "
  819.                                     " COORDINATES IN METERS\n\n"
  820.                                     "                        "
  821.                                     " I+ AND I- INDICATE THE SEGMENTS BEFORE AND AFTER I\n\n");
  822.  
  823.                                 fprintf (
  824.                                     output_fp,
  825.                                     "\n"
  826.                                     "  SEG.   COORDINATES OF SEG. CENTER     SEG.     "
  827.                                     "ORIENTATION ANGLES    WIRE    CONNECTION DATA   TAG\n"
  828.                                     "  NO.       X         Y         Z       LENGTH   "
  829.                                     "  ALPHA     BETA      RADIUS    I-   I    I+    NO.");
  830.  
  831.                                 for (i = 0; i < data.n; i++)
  832.                                 {
  833.                                         xw1 = data.x2[i] - data.x1[i];
  834.                                         yw1 = data.y2[i] - data.y1[i];
  835.                                         zw1 = data.z2[i] - data.z1[i];
  836.                                         data.x[i] = (data.x1[i] + data.x2[i]) / 2.;
  837.                                         data.y[i] = (data.y1[i] + data.y2[i]) / 2.;
  838.                                         data.z[i] = (data.z1[i] + data.z2[i]) / 2.;
  839.                                         xw2 = xw1 * xw1 + yw1 * yw1 + zw1 * zw1;
  840.                                         yw2 = sqrtl (xw2);
  841.                                         yw2 = (xw2 / yw2 + yw2) * .5;
  842.                                         data.si[i] = yw2;
  843.                                         data.cab[i] = xw1 / yw2;
  844.                                         data.sab[i] = yw1 / yw2;
  845.                                         xw2 = zw1 / yw2;
  846.  
  847.                                         if (xw2 > 1.)
  848.                                                 xw2 = 1.;
  849.                                         if (xw2 < -1.)
  850.                                                 xw2 = -1.;
  851.  
  852.                                         data.salp[i] = xw2;
  853.                                         xw2 = asinl (xw2) * TD;
  854.                                         yw2 = atan2l (yw1, xw1) * TD;
  855.  
  856.                                         fprintf (
  857.                                             output_fp,
  858.                                             "\n"
  859.                                             " %5d%10.5f%10.5f%10.5f%10.5f"
  860.                                             " %10.5f%10.5f%10.5f%7d%7d%7d%7d",
  861.                                             i + 1,
  862.                                             data.x[i],
  863.                                             data.y[i],
  864.                                             data.z[i],
  865.                                             data.si[i],
  866.                                             xw2,
  867.                                             yw2,
  868.                                             data.bi[i],
  869.                                             data.icon1[i],
  870.                                             i + 1,
  871.                                             data.icon2[i],
  872.                                             data.itag[i]);
  873.  
  874.                                         if (plot.iplp1 == 1)
  875.                                                 fprintf (
  876.                                                     plot_fp,
  877.                                                     "%12.5E %12.5E %12.5E "
  878.                                                     "%12.5E %12.5E %12.5E %12.5E%7d%7d%7d\n",
  879.                                                     data.x[i],
  880.                                                     data.y[i],
  881.                                                     data.z[i],
  882.                                                     data.si[i],
  883.                                                     xw2,
  884.                                                     yw2,
  885.                                                     data.bi[i],
  886.                                                     data.icon1[i],
  887.                                                     i + 1,
  888.                                                     data.icon2[i]);
  889.  
  890.                                         if ((data.si[i] <= 1.e-20) || (data.bi[i] <= 0.))
  891.                                         {
  892.                                                 fprintf (output_fp, "\n SEGMENT DATA ERROR");
  893.                                                 stop (-1);
  894.                                         }
  895.  
  896.                                 } /* for( i = 0; i < data.n; i++ ) */
  897.  
  898.                         } /* if( data.n != 0) */
  899.  
  900.                         fprintf (output_fp, "\n\n\n");
  901.                         if (data.m != 0)
  902.                         {
  903.                                 fprintf (
  904.                                     output_fp,
  905.                                     "\n\n\n"
  906.                                     "                                   "
  907.                                     " - - - - - SURFACE PATCH DATA - - - -\n"
  908.                                     "                                            "
  909.                                     " COORDINATES IN METERS\n\n"
  910.                                     " PATCH      COORD. OF PATCH CENTER           UNIT NORMAL "
  911.                                     "VECTOR      "
  912.                                     " PATCH           COMPONENTS OF UNIT TANGENT VECTORS\n"
  913.                                     "  No:       X          Y          Z          X        Y  "
  914.                                     "      Z      "
  915.                                     " AREA         X1       Y1       Z1        X2       Y2    "
  916.                                     "  Z2");
  917.  
  918.                                 for (i = 0; i < data.m; i++)
  919.                                 {
  920.                                         xw1 = (data.t1y[i] * data.t2z[i] -
  921.                                                data.t1z[i] * data.t2y[i]) *
  922.                                               data.psalp[i];
  923.                                         yw1 = (data.t1z[i] * data.t2x[i] -
  924.                                                data.t1x[i] * data.t2z[i]) *
  925.                                               data.psalp[i];
  926.                                         zw1 = (data.t1x[i] * data.t2y[i] -
  927.                                                data.t1y[i] * data.t2x[i]) *
  928.                                               data.psalp[i];
  929.  
  930.                                         fprintf (
  931.                                             output_fp,
  932.                                             "\n"
  933.                                             " %4d %10.5f %10.5f %10.5f  %8.5f %8.5f %8.5f"
  934.                                             " %10.5f  %8.5f %8.5f %8.5f  %8.5f %8.5f %8.5f",
  935.                                             i + 1,
  936.                                             data.px[i],
  937.                                             data.py[i],
  938.                                             data.pz[i],
  939.                                             xw1,
  940.                                             yw1,
  941.                                             zw1,
  942.                                             data.pbi[i],
  943.                                             data.t1x[i],
  944.                                             data.t1y[i],
  945.                                             data.t1z[i],
  946.                                             data.t2x[i],
  947.                                             data.t2y[i],
  948.                                             data.t2z[i]);
  949.  
  950.                                 } /* for( i = 0; i < data.m; i++ ) */
  951.  
  952.                         } /* if( data.m == 0) */
  953.  
  954.                         data.npm = data.n + data.m;
  955.                         data.np2m = data.n + 2 * data.m;
  956.                         data.np3m = data.n + 3 * data.m;
  957.  
  958.                         return;
  959.  
  960.                         /* "gm" card, move structure or reproduce */
  961.                         /* original structure in new positions.   */
  962.                 case 5:
  963.  
  964.                         fprintf (
  965.                             output_fp,
  966.                             "\n     THE STRUCTURE HAS BEEN MOVED, MOVE DATA CARD IS:\n"
  967.                             "   %3d %5d %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f",
  968.                             itg,
  969.                             ns,
  970.                             xw1,
  971.                             yw1,
  972.                             zw1,
  973.                             xw2,
  974.                             yw2,
  975.                             zw2,
  976.                             rad);
  977.  
  978.                         xw1 = xw1 * TA;
  979.                         yw1 = yw1 * TA;
  980.                         zw1 = zw1 * TA;
  981.  
  982.                         move (xw1, yw1, zw1, xw2, yw2, zw2, (int) (rad + .5), ns, itg);
  983.                         continue;
  984.  
  985.                 case 6: /* "sp" card, generate single new patch */
  986.  
  987.                         i1 = data.m + 1;
  988.                         ns++;
  989.  
  990.                         if (itg != 0)
  991.                         {
  992.                                 fprintf (output_fp, "\n  PATCH DATA ERROR");
  993.                                 stop (-1);
  994.                         }
  995.  
  996.                         fprintf (
  997.                             output_fp,
  998.                             "\n"
  999.                             " %5d%c %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f",
  1000.                             i1,
  1001.                             ipt[ns - 1],
  1002.                             xw1,
  1003.                             yw1,
  1004.                             zw1,
  1005.                             xw2,
  1006.                             yw2,
  1007.                             zw2);
  1008.  
  1009.                         if ((ns == 2) || (ns == 4))
  1010.                                 isct = 1;
  1011.  
  1012.                         if (ns > 1)
  1013.                         {
  1014.                                 readgm (gm, &ix, &iy, &x3, &y3, &z3, &x4, &y4, &z4, &dummy);
  1015.  
  1016.                                 if ((ns == 2) || (itg > 0))
  1017.                                 {
  1018.                                         x4 = xw1 + x3 - xw2;
  1019.                                         y4 = yw1 + y3 - yw2;
  1020.                                         z4 = zw1 + z3 - zw2;
  1021.                                 }
  1022.  
  1023.                                 fprintf (
  1024.                                     output_fp,
  1025.                                     "\n"
  1026.                                     "      %11.5f %11.5f %11.5f %11.5f %11.5f %11.5f",
  1027.                                     x3,
  1028.                                     y3,
  1029.                                     z3,
  1030.                                     x4,
  1031.                                     y4,
  1032.                                     z4);
  1033.  
  1034.                                 if (strcmp (gm, "SC") != 0)
  1035.                                 {
  1036.                                         fprintf (output_fp, "\n  PATCH DATA ERROR");
  1037.                                         stop (-1);
  1038.                                 }
  1039.  
  1040.                         } /* if( ns > 1) */
  1041.                         else
  1042.                         {
  1043.                                 xw2 = xw2 * TA;
  1044.                                 yw2 = yw2 * TA;
  1045.                         }
  1046.  
  1047.                         patch (itg, ns, xw1, yw1, zw1, xw2, yw2, zw2, x3, y3, z3, x4, y4, z4);
  1048.  
  1049.                         continue;
  1050.  
  1051.                 case 7: /* "sm" card, generate multiple-patch surface */
  1052.  
  1053.                         i1 = data.m + 1;
  1054.                         fprintf (
  1055.                             output_fp,
  1056.                             "\n"
  1057.                             " %5d%c %10.5f %11.5f %11.5f %11.5f %11.5f %11.5f"
  1058.                             "     SURFACE - %d BY %d PATCHES",
  1059.                             i1,
  1060.                             ipt[1],
  1061.                             xw1,
  1062.                             yw1,
  1063.                             zw1,
  1064.                             xw2,
  1065.                             yw2,
  1066.                             zw2,
  1067.                             itg,
  1068.                             ns);
  1069.  
  1070.                         if ((itg < 1) || (ns < 1))
  1071.                         {
  1072.                                 fprintf (output_fp, "\n  PATCH DATA ERROR");
  1073.                                 stop (-1);
  1074.                         }
  1075.  
  1076.                         readgm (gm, &ix, &iy, &x3, &y3, &z3, &x4, &y4, &z4, &dummy);
  1077.  
  1078.                         if ((ns == 2) || (itg > 0))
  1079.                         {
  1080.                                 x4 = xw1 + x3 - xw2;
  1081.                                 y4 = yw1 + y3 - yw2;
  1082.                                 z4 = zw1 + z3 - zw2;
  1083.                         }
  1084.  
  1085.                         fprintf (
  1086.                             output_fp,
  1087.                             "\n"
  1088.                             "      %11.5f %11.5f %11.5f %11.5f %11.5f %11.5f",
  1089.                             x3,
  1090.                             y3,
  1091.                             z3,
  1092.                             x4,
  1093.                             y4,
  1094.                             z4);
  1095.  
  1096.                         if (strcmp (gm, "SC") != 0)
  1097.                         {
  1098.                                 fprintf (output_fp, "\n  PATCH DATA ERROR");
  1099.                                 stop (-1);
  1100.                         }
  1101.  
  1102.                         patch (itg, ns, xw1, yw1, zw1, xw2, yw2, zw2, x3, y3, z3, x4, y4, z4);
  1103.  
  1104.                         continue;
  1105.  
  1106.                 case 8: /* "ga" card, generate segment data for wire arc */
  1107.  
  1108.                         nwire++;
  1109.                         i1 = data.n + 1;
  1110.                         i2 = data.n + ns;
  1111.  
  1112.                         fprintf (
  1113.                             output_fp,
  1114.                             "\n"
  1115.                             " %5d  ARC RADIUS: %9.5f  FROM: %8.3f TO: %8.3f DEGREES"
  1116.                             "       %11.5f %5d %5d %5d %4d",
  1117.                             nwire,
  1118.                             xw1,
  1119.                             yw1,
  1120.                             zw1,
  1121.                             xw2,
  1122.                             ns,
  1123.                             i1,
  1124.                             i2,
  1125.                             itg);
  1126.  
  1127.                         arc (itg, ns, xw1, yw1, zw1, xw2);
  1128.  
  1129.                         continue;
  1130.  
  1131.                 case 9: /* "sc" card */
  1132.  
  1133.                         if (isct == 0)
  1134.                         {
  1135.                                 fprintf (output_fp, "\n  PATCH DATA ERROR");
  1136.                                 stop (-1);
  1137.                         }
  1138.  
  1139.                         i1 = data.m + 1;
  1140.                         ns++;
  1141.  
  1142.                         if ((itg != 0) || ((ns != 2) && (ns != 4)))
  1143.                         {
  1144.                                 fprintf (output_fp, "\n  PATCH DATA ERROR");
  1145.                                 stop (-1);
  1146.                         }
  1147.  
  1148.                         xs1 = x4;
  1149.                         ys1 = y4;
  1150.                         zs1 = z4;
  1151.                         xs2 = x3;
  1152.                         ys2 = y3;
  1153.                         zs2 = z3;
  1154.                         x3 = xw1;
  1155.                         y3 = yw1;
  1156.                         z3 = zw1;
  1157.  
  1158.                         if (ns == 4)
  1159.                         {
  1160.                                 x4 = xw2;
  1161.                                 y4 = yw2;
  1162.                                 z4 = zw2;
  1163.                         }
  1164.  
  1165.                         xw1 = xs1;
  1166.                         yw1 = ys1;
  1167.                         zw1 = zs1;
  1168.                         xw2 = xs2;
  1169.                         yw2 = ys2;
  1170.                         zw2 = zs2;
  1171.  
  1172.                         if (ns != 4)
  1173.                         {
  1174.                                 x4 = xw1 + x3 - xw2;
  1175.                                 y4 = yw1 + y3 - yw2;
  1176.                                 z4 = zw1 + z3 - zw2;
  1177.                         }
  1178.  
  1179.                         fprintf (
  1180.                             output_fp,
  1181.                             "\n"
  1182.                             " %5d%c %10.5f %11.5f %11.5f %11.5f %11.5f %11.5f",
  1183.                             i1,
  1184.                             ipt[ns - 1],
  1185.                             xw1,
  1186.                             yw1,
  1187.                             zw1,
  1188.                             xw2,
  1189.                             yw2,
  1190.                             zw2);
  1191.  
  1192.                         fprintf (
  1193.                             output_fp,
  1194.                             "\n"
  1195.                             "      %11.5f %11.5f %11.5f  %11.5f %11.5f %11.5f",
  1196.                             x3,
  1197.                             y3,
  1198.                             z3,
  1199.                             x4,
  1200.                             y4,
  1201.                             z4);
  1202.  
  1203.                         patch (itg, ns, xw1, yw1, zw1, xw2, yw2, zw2, x3, y3, z3, x4, y4, z4);
  1204.  
  1205.                         continue;
  1206.  
  1207.                 case 10: /* "gh" card, generate helix */
  1208.  
  1209.                         nwire++;
  1210.                         i1 = data.n + 1;
  1211.                         i2 = data.n + ns;
  1212.  
  1213.                         fprintf (
  1214.                             output_fp,
  1215.                             "\n"
  1216.                             " %5d HELIX STRUCTURE - SPACING OF TURNS: %8.3f AXIAL"
  1217.                             " LENGTH: %8.3f  %8.3f %5d %5d %5d %4d\n      "
  1218.                             " RADIUS X1:%8.3f Y1:%8.3f X2:%8.3f Y2:%8.3f ",
  1219.                             nwire,
  1220.                             xw1,
  1221.                             yw1,
  1222.                             rad,
  1223.                             ns,
  1224.                             i1,
  1225.                             i2,
  1226.                             itg,
  1227.                             zw1,
  1228.                             xw2,
  1229.                             yw2,
  1230.                             zw2);
  1231.  
  1232.                         helix (xw1, yw1, zw1, xw2, yw2, zw2, rad, ns, itg);
  1233.  
  1234.                         continue;
  1235.  
  1236.                 case 11: /* "gf" card, not supported */
  1237.                         abort_on_error (-5);
  1238.                         continue;
  1239.                 default: /* error message */
  1240.  
  1241.                         fprintf (output_fp, "\n  GEOMETRY DATA CARD ERROR");
  1242.                         fprintf (
  1243.                             output_fp,
  1244.                             "\n"
  1245.                             " %2s %3d %5d %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f",
  1246.                             gm,
  1247.                             itg,
  1248.                             ns,
  1249.                             xw1,
  1250.                             yw1,
  1251.                             zw1,
  1252.                             xw2,
  1253.                             yw2,
  1254.                             zw2,
  1255.                             rad);
  1256.  
  1257.                         stop (-1);
  1258.  
  1259.                 } /* switch( gm_num ) */
  1260.  
  1261.         } /* do */
  1262.         while (TRUE);
  1263.  
  1264.         return;
  1265. }
  1266.  
  1267. /*-----------------------------------------------------------------------*/
  1268.  
  1269. /* subroutine helix generates segment geometry */
  1270. /* data for a helix of ns segments */
  1271. void helix (
  1272.     double s,
  1273.     double hl,
  1274.     double a1,
  1275.     double b1,
  1276.     double a2,
  1277.     double b2,
  1278.     double rad,
  1279.     int ns,
  1280.     int itg)
  1281. {
  1282.         int ist, i, mreq;
  1283.         // double turns;
  1284.         double zinc, copy, sangle, hdia, turn, pitch, hmaj, hmin;
  1285.  
  1286.         ist = data.n;
  1287.         data.n += ns;
  1288.         data.np = data.n;
  1289.         data.mp = data.m;
  1290.         data.ipsym = 0;
  1291.  
  1292.         if (ns < 1)
  1293.                 return;
  1294.  
  1295.         // turns= fabsl( hl/ s);
  1296.         zinc = fabsl (hl / ns);
  1297.  
  1298.         /* Reallocate tags buffer */
  1299.         mem_realloc ((void *) &data.itag, (data.n + data.m) * sizeof (int)); /*????*/
  1300.  
  1301.         /* Reallocate wire buffers */
  1302.         mreq = data.n * sizeof (double);
  1303.         mem_realloc ((void *) &data.x1, mreq);
  1304.         mem_realloc ((void *) &data.y1, mreq);
  1305.         mem_realloc ((void *) &data.z1, mreq);
  1306.         mem_realloc ((void *) &data.x2, mreq);
  1307.         mem_realloc ((void *) &data.y2, mreq);
  1308.         mem_realloc ((void *) &data.z2, mreq);
  1309.         mem_realloc ((void *) &data.bi, mreq);
  1310.  
  1311.         data.z1[ist] = 0.;
  1312.         for (i = ist; i < data.n; i++)
  1313.         {
  1314.                 data.bi[i] = rad;
  1315.                 data.itag[i] = itg;
  1316.  
  1317.                 if (i != ist)
  1318.                         data.z1[i] = data.z1[i - 1] + zinc;
  1319.  
  1320.                 data.z2[i] = data.z1[i] + zinc;
  1321.  
  1322.                 if (a2 == a1)
  1323.                 {
  1324.                         if (b1 == 0.)
  1325.                                 b1 = a1;
  1326.  
  1327.                         data.x1[i] = a1 * cosl (2. * PI * data.z1[i] / s);
  1328.                         data.y1[i] = b1 * sinl (2. * PI * data.z1[i] / s);
  1329.                         data.x2[i] = a1 * cosl (2. * PI * data.z2[i] / s);
  1330.                         data.y2[i] = b1 * sinl (2. * PI * data.z2[i] / s);
  1331.                 }
  1332.                 else
  1333.                 {
  1334.                         if (b2 == 0.)
  1335.                                 b2 = a2;
  1336.  
  1337.                         data.x1[i] = (a1 + (a2 - a1) * data.z1[i] / fabsl (hl)) *
  1338.                                      cosl (2. * PI * data.z1[i] / s);
  1339.                         data.y1[i] = (b1 + (b2 - b1) * data.z1[i] / fabsl (hl)) *
  1340.                                      sinl (2. * PI * data.z1[i] / s);
  1341.                         data.x2[i] = (a1 + (a2 - a1) * data.z2[i] / fabsl (hl)) *
  1342.                                      cosl (2. * PI * data.z2[i] / s);
  1343.                         data.y2[i] = (b1 + (b2 - b1) * data.z2[i] / fabsl (hl)) *
  1344.                                      sinl (2. * PI * data.z2[i] / s);
  1345.  
  1346.                 } /* if( a2 == a1) */
  1347.  
  1348.                 if (hl > 0.)
  1349.                         continue;
  1350.  
  1351.                 copy = data.x1[i];
  1352.                 data.x1[i] = data.y1[i];
  1353.                 data.y1[i] = copy;
  1354.                 copy = data.x2[i];
  1355.                 data.x2[i] = data.y2[i];
  1356.                 data.y2[i] = copy;
  1357.  
  1358.         } /* for( i = ist; i < data.n; i++ ) */
  1359.  
  1360.         if (a2 != a1)
  1361.         {
  1362.                 sangle = atanl (a2 / (fabsl (hl) + (fabsl (hl) * a1) / (a2 - a1)));
  1363.                 fprintf (output_fp, "\n       THE CONE ANGLE OF THE SPIRAL IS %10.4f", sangle);
  1364.                 return;
  1365.         }
  1366.  
  1367.         if (a1 == b1)
  1368.         {
  1369.                 hdia = 2. * a1;
  1370.                 turn = hdia * PI;
  1371.                 pitch = atanl (s / (PI * hdia));
  1372.                 turn = turn / cosl (pitch);
  1373.                 pitch = 180. * pitch / PI;
  1374.         }
  1375.         else
  1376.         {
  1377.                 if (a1 >= b1)
  1378.                 {
  1379.                         hmaj = 2. * a1;
  1380.                         hmin = 2. * b1;
  1381.                 }
  1382.                 else
  1383.                 {
  1384.                         hmaj = 2. * b1;
  1385.                         hmin = 2. * a1;
  1386.                 }
  1387.  
  1388.                 hdia = sqrtl ((hmaj * hmaj + hmin * hmin) / 2 * hmaj);
  1389.                 turn = 2. * PI * hdia;
  1390.                 pitch = (180. / PI) * atanl (s / (PI * hdia));
  1391.  
  1392.         } /* if( a1 == b1) */
  1393.  
  1394.         fprintf (
  1395.             output_fp,
  1396.             "\n"
  1397.             "       THE PITCH ANGLE IS: %.4f    THE LENGTH OF WIRE/TURN IS: %.4f",
  1398.             pitch,
  1399.             turn);
  1400.  
  1401.         return;
  1402. }
  1403.  
  1404. /*-----------------------------------------------------------------------*/
  1405.  
  1406. /* isegno returns the segment number of the mth segment having the */
  1407. /* tag number itagi.  if itagi=0 segment number m is returned. */
  1408. int isegno (int itagi, int mx)
  1409. {
  1410.         int icnt, i, iseg;
  1411.  
  1412.         if (mx <= 0)
  1413.         {
  1414.                 fprintf (
  1415.                     output_fp,
  1416.                     "\n  CHECK DATA, PARAMETER SPECIFYING SEGMENT"
  1417.                     " POSITION IN A GROUP OF EQUAL TAGS MUST NOT BE ZERO");
  1418.                 stop (-1);
  1419.         }
  1420.  
  1421.         icnt = 0;
  1422.         if (itagi == 0)
  1423.         {
  1424.                 iseg = mx;
  1425.                 return (iseg);
  1426.         }
  1427.  
  1428.         if (data.n > 0)
  1429.         {
  1430.                 for (i = 0; i < data.n; i++)
  1431.                 {
  1432.                         if (data.itag[i] != itagi)
  1433.                                 continue;
  1434.  
  1435.                         icnt++;
  1436.                         if (icnt == mx)
  1437.                         {
  1438.                                 iseg = i + 1;
  1439.                                 return (iseg);
  1440.                         }
  1441.  
  1442.                 } /* for( i = 0; i < data.n; i++ ) */
  1443.  
  1444.         } /* if( data.n > 0) */
  1445.  
  1446.         fprintf (
  1447.             output_fp,
  1448.             "\n\n"
  1449.             "  NO SEGMENT HAS AN ITAG OF %d",
  1450.             itagi);
  1451.         stop (-1);
  1452.  
  1453.         return (0);
  1454. }
  1455.  
  1456. /*-----------------------------------------------------------------------*/
  1457.  
  1458. /* subroutine move moves the structure with respect to its */
  1459. /* coordinate system or reproduces structure in new positions. */
  1460. /* structure is rotated about x,y,z axes by rox,roy,roz */
  1461. /* respectively, then shifted by xs,ys,zs */
  1462. void move (
  1463.     double rox,
  1464.     double roy,
  1465.     double roz,
  1466.     double xs,
  1467.     double ys,
  1468.     double zs,
  1469.     int its,
  1470.     int nrpt,
  1471.     int itgi)
  1472. {
  1473.         int nrp, ix, i1, k, ir, i, ii, mreq;
  1474.         double sps, cps, sth, cth, sph, cph, xx, xy;
  1475.         double xz, yx, yy, yz, zx, zy, zz, xi, yi, zi;
  1476.  
  1477.         if (fabsl (rox) + fabsl (roy) > 1.0e-10)
  1478.                 data.ipsym = data.ipsym * 3;
  1479.  
  1480.         sps = sinl (rox);
  1481.         cps = cosl (rox);
  1482.         sth = sinl (roy);
  1483.         cth = cosl (roy);
  1484.         sph = sinl (roz);
  1485.         cph = cosl (roz);
  1486.         xx = cph * cth;
  1487.         xy = cph * sth * sps - sph * cps;
  1488.         xz = cph * sth * cps + sph * sps;
  1489.         yx = sph * cth;
  1490.         yy = sph * sth * sps + cph * cps;
  1491.         yz = sph * sth * cps - cph * sps;
  1492.         zx = -sth;
  1493.         zy = cth * sps;
  1494.         zz = cth * cps;
  1495.  
  1496.         if (nrpt == 0)
  1497.                 nrp = 1;
  1498.         else
  1499.                 nrp = nrpt;
  1500.  
  1501.         ix = 1;
  1502.         if (data.n > 0)
  1503.         {
  1504.                 i1 = isegno (its, 1);
  1505.                 if (i1 < 1)
  1506.                         i1 = 1;
  1507.  
  1508.                 ix = i1;
  1509.                 if (nrpt == 0)
  1510.                         k = i1 - 1;
  1511.                 else
  1512.                 {
  1513.                         k = data.n;
  1514.                         /* Reallocate tags buffer */
  1515.                         mreq = data.n + data.m + (data.n + 1 - i1) * nrpt;
  1516.                         mem_realloc ((void *) &data.itag, mreq * sizeof (int));
  1517.  
  1518.                         /* Reallocate wire buffers */
  1519.                         mreq = (data.n + (data.n + 1 - i1) * nrpt) * sizeof (double);
  1520.                         mem_realloc ((void *) &data.x1, mreq);
  1521.                         mem_realloc ((void *) &data.y1, mreq);
  1522.                         mem_realloc ((void *) &data.z1, mreq);
  1523.                         mem_realloc ((void *) &data.x2, mreq);
  1524.                         mem_realloc ((void *) &data.y2, mreq);
  1525.                         mem_realloc ((void *) &data.z2, mreq);
  1526.                         mem_realloc ((void *) &data.bi, mreq);
  1527.                 }
  1528.  
  1529.                 for (ir = 0; ir < nrp; ir++)
  1530.                 {
  1531.                         for (i = i1 - 1; i < data.n; i++)
  1532.                         {
  1533.                                 xi = data.x1[i];
  1534.                                 yi = data.y1[i];
  1535.                                 zi = data.z1[i];
  1536.                                 data.x1[k] = xi * xx + yi * xy + zi * xz + xs;
  1537.                                 data.y1[k] = xi * yx + yi * yy + zi * yz + ys;
  1538.                                 data.z1[k] = xi * zx + yi * zy + zi * zz + zs;
  1539.                                 xi = data.x2[i];
  1540.                                 yi = data.y2[i];
  1541.                                 zi = data.z2[i];
  1542.                                 data.x2[k] = xi * xx + yi * xy + zi * xz + xs;
  1543.                                 data.y2[k] = xi * yx + yi * yy + zi * yz + ys;
  1544.                                 data.z2[k] = xi * zx + yi * zy + zi * zz + zs;
  1545.                                 data.bi[k] = data.bi[i];
  1546.                                 data.itag[k] = data.itag[i];
  1547.                                 if (data.itag[i] != 0)
  1548.                                         data.itag[k] = data.itag[i] + itgi;
  1549.  
  1550.                                 k++;
  1551.  
  1552.                         } /* for( i = i1; i < data.n; i++ ) */
  1553.  
  1554.                         i1 = data.n + 1;
  1555.                         data.n = k;
  1556.  
  1557.                 } /* for( ir = 0; ir < nrp; ir++ ) */
  1558.  
  1559.         } /* if( data.n >= n2) */
  1560.  
  1561.         if (data.m > 0)
  1562.         {
  1563.                 i1 = 0;
  1564.                 if (nrpt == 0)
  1565.                         k = 0;
  1566.                 else
  1567.                         k = data.m;
  1568.  
  1569.                 /* Reallocate patch buffers */
  1570.                 mreq = data.m * (1 + nrpt) * sizeof (double);
  1571.                 mem_realloc ((void *) &data.px, mreq);
  1572.                 mem_realloc ((void *) &data.py, mreq);
  1573.                 mem_realloc ((void *) &data.pz, mreq);
  1574.                 mem_realloc ((void *) &data.t1x, mreq);
  1575.                 mem_realloc ((void *) &data.t1y, mreq);
  1576.                 mem_realloc ((void *) &data.t1z, mreq);
  1577.                 mem_realloc ((void *) &data.t2x, mreq);
  1578.                 mem_realloc ((void *) &data.t2y, mreq);
  1579.                 mem_realloc ((void *) &data.t2z, mreq);
  1580.                 mem_realloc ((void *) &data.pbi, mreq);
  1581.                 mem_realloc ((void *) &data.psalp, mreq);
  1582.  
  1583.                 for (ii = 0; ii < nrp; ii++)
  1584.                 {
  1585.                         for (i = i1; i < data.m; i++)
  1586.                         {
  1587.                                 xi = data.px[i];
  1588.                                 yi = data.py[i];
  1589.                                 zi = data.pz[i];
  1590.                                 data.px[k] = xi * xx + yi * xy + zi * xz + xs;
  1591.                                 data.py[k] = xi * yx + yi * yy + zi * yz + ys;
  1592.                                 data.pz[k] = xi * zx + yi * zy + zi * zz + zs;
  1593.                                 xi = data.t1x[i];
  1594.                                 yi = data.t1y[i];
  1595.                                 zi = data.t1z[i];
  1596.                                 data.t1x[k] = xi * xx + yi * xy + zi * xz;
  1597.                                 data.t1y[k] = xi * yx + yi * yy + zi * yz;
  1598.                                 data.t1z[k] = xi * zx + yi * zy + zi * zz;
  1599.                                 xi = data.t2x[i];
  1600.                                 yi = data.t2y[i];
  1601.                                 zi = data.t2z[i];
  1602.                                 data.t2x[k] = xi * xx + yi * xy + zi * xz;
  1603.                                 data.t2y[k] = xi * yx + yi * yy + zi * yz;
  1604.                                 data.t2z[k] = xi * zx + yi * zy + zi * zz;
  1605.                                 data.psalp[k] = data.psalp[i];
  1606.                                 data.pbi[k] = data.pbi[i];
  1607.                                 k++;
  1608.  
  1609.                         } /* for( i = i1; i < data.m; i++ ) */
  1610.  
  1611.                         i1 = data.m;
  1612.                         data.m = k;
  1613.  
  1614.                 } /* for( ii = 0; ii < nrp; ii++ ) */
  1615.  
  1616.         } /* if( data.m >= m2) */
  1617.  
  1618.         if ((nrpt == 0) && (ix == 1))
  1619.                 return;
  1620.  
  1621.         data.np = data.n;
  1622.         data.mp = data.m;
  1623.         data.ipsym = 0;
  1624.  
  1625.         return;
  1626. }
  1627.  
  1628. /*-----------------------------------------------------------------------*/
  1629.  
  1630. /* patch generates and modifies patch geometry data */
  1631. void patch (
  1632.     int nx,
  1633.     int ny,
  1634.     double ax1,
  1635.     double ay1,
  1636.     double az1,
  1637.     double ax2,
  1638.     double ay2,
  1639.     double az2,
  1640.     double ax3,
  1641.     double ay3,
  1642.     double az3,
  1643.     double ax4,
  1644.     double ay4,
  1645.     double az4)
  1646. {
  1647.         int mi, ntp, iy, ix, mreq;
  1648.         double s1x = 0., s1y = 0., s1z = 0., s2x = 0., s2y = 0., s2z = 0., xst = 0.;
  1649.         double znv, xnv, ynv, xa, xn2, yn2, zn2, salpn, xs, ys, zs, xt, yt, zt;
  1650.  
  1651.         /* new patches.  for nx=0, ny=1,2,3,4 patch is (respectively) */
  1652.         /* arbitrary, rectagular, triangular, or quadrilateral. */
  1653.         /* for nx and ny  > 0 a rectangular surface is produced with */
  1654.         /* nx by ny rectangular patches. */
  1655.  
  1656.         data.m++;
  1657.         mi = data.m - 1;
  1658.  
  1659.         /* Reallocate patch buffers */
  1660.         mreq = data.m * sizeof (double);
  1661.         mem_realloc ((void *) &data.px, mreq);
  1662.         mem_realloc ((void *) &data.py, mreq);
  1663.         mem_realloc ((void *) &data.pz, mreq);
  1664.         mem_realloc ((void *) &data.t1x, mreq);
  1665.         mem_realloc ((void *) &data.t1y, mreq);
  1666.         mem_realloc ((void *) &data.t1z, mreq);
  1667.         mem_realloc ((void *) &data.t2x, mreq);
  1668.         mem_realloc ((void *) &data.t2y, mreq);
  1669.         mem_realloc ((void *) &data.t2z, mreq);
  1670.         mem_realloc ((void *) &data.pbi, mreq);
  1671.         mem_realloc ((void *) &data.psalp, mreq);
  1672.  
  1673.         if (nx > 0)
  1674.                 ntp = 2;
  1675.         else
  1676.                 ntp = ny;
  1677.  
  1678.         if (ntp <= 1)
  1679.         {
  1680.                 data.px[mi] = ax1;
  1681.                 data.py[mi] = ay1;
  1682.                 data.pz[mi] = az1;
  1683.                 data.pbi[mi] = az2;
  1684.                 znv = cosl (ax2);
  1685.                 xnv = znv * cosl (ay2);
  1686.                 ynv = znv * sinl (ay2);
  1687.                 znv = sinl (ax2);
  1688.                 xa = sqrtl (xnv * xnv + ynv * ynv);
  1689.  
  1690.                 if (xa >= 1.0e-6)
  1691.                 {
  1692.                         data.t1x[mi] = -ynv / xa;
  1693.                         data.t1y[mi] = xnv / xa;
  1694.                         data.t1z[mi] = 0.;
  1695.                 }
  1696.                 else
  1697.                 {
  1698.                         data.t1x[mi] = 1.;
  1699.                         data.t1y[mi] = 0.;
  1700.                         data.t1z[mi] = 0.;
  1701.                 }
  1702.  
  1703.         } /* if( ntp <= 1) */
  1704.         else
  1705.         {
  1706.                 s1x = ax2 - ax1;
  1707.                 s1y = ay2 - ay1;
  1708.                 s1z = az2 - az1;
  1709.                 s2x = ax3 - ax2;
  1710.                 s2y = ay3 - ay2;
  1711.                 s2z = az3 - az2;
  1712.  
  1713.                 if (nx != 0)
  1714.                 {
  1715.                         s1x = s1x / nx;
  1716.                         s1y = s1y / nx;
  1717.                         s1z = s1z / nx;
  1718.                         s2x = s2x / ny;
  1719.                         s2y = s2y / ny;
  1720.                         s2z = s2z / ny;
  1721.                 }
  1722.  
  1723.                 xnv = s1y * s2z - s1z * s2y;
  1724.                 ynv = s1z * s2x - s1x * s2z;
  1725.                 znv = s1x * s2y - s1y * s2x;
  1726.                 xa = sqrtl (xnv * xnv + ynv * ynv + znv * znv);
  1727.                 xnv = xnv / xa;
  1728.                 ynv = ynv / xa;
  1729.                 znv = znv / xa;
  1730.                 xst = sqrtl (s1x * s1x + s1y * s1y + s1z * s1z);
  1731.                 data.t1x[mi] = s1x / xst;
  1732.                 data.t1y[mi] = s1y / xst;
  1733.                 data.t1z[mi] = s1z / xst;
  1734.  
  1735.                 if (ntp <= 2)
  1736.                 {
  1737.                         data.px[mi] = ax1 + .5 * (s1x + s2x);
  1738.                         data.py[mi] = ay1 + .5 * (s1y + s2y);
  1739.                         data.pz[mi] = az1 + .5 * (s1z + s2z);
  1740.                         data.pbi[mi] = xa;
  1741.                 }
  1742.                 else
  1743.                 {
  1744.                         if (ntp != 4)
  1745.                         {
  1746.                                 data.px[mi] = (ax1 + ax2 + ax3) / 3.;
  1747.                                 data.py[mi] = (ay1 + ay2 + ay3) / 3.;
  1748.                                 data.pz[mi] = (az1 + az2 + az3) / 3.;
  1749.                                 data.pbi[mi] = .5 * xa;
  1750.                         }
  1751.                         else
  1752.                         {
  1753.                                 s1x = ax3 - ax1;
  1754.                                 s1y = ay3 - ay1;
  1755.                                 s1z = az3 - az1;
  1756.                                 s2x = ax4 - ax1;
  1757.                                 s2y = ay4 - ay1;
  1758.                                 s2z = az4 - az1;
  1759.                                 xn2 = s1y * s2z - s1z * s2y;
  1760.                                 yn2 = s1z * s2x - s1x * s2z;
  1761.                                 zn2 = s1x * s2y - s1y * s2x;
  1762.                                 xst = sqrtl (xn2 * xn2 + yn2 * yn2 + zn2 * zn2);
  1763.                                 salpn = 1. / (3. * (xa + xst));
  1764.                                 data.px[mi] =
  1765.                                     (xa * (ax1 + ax2 + ax3) + xst * (ax1 + ax3 + ax4)) * salpn;
  1766.                                 data.py[mi] =
  1767.                                     (xa * (ay1 + ay2 + ay3) + xst * (ay1 + ay3 + ay4)) * salpn;
  1768.                                 data.pz[mi] =
  1769.                                     (xa * (az1 + az2 + az3) + xst * (az1 + az3 + az4)) * salpn;
  1770.                                 data.pbi[mi] = .5 * (xa + xst);
  1771.                                 s1x = (xnv * xn2 + ynv * yn2 + znv * zn2) / xst;
  1772.  
  1773.                                 if (s1x <= 0.9998)
  1774.                                 {
  1775.                                         fprintf (
  1776.                                             output_fp,
  1777.                                             "\n  ERROR -- CORNERS OF QUADRILATERAL"
  1778.                                             " PATCH DO NOT LIE IN A PLANE");
  1779.                                         stop (-1);
  1780.                                 }
  1781.  
  1782.                         } /* if( ntp != 4) */
  1783.  
  1784.                 } /* if( ntp <= 2) */
  1785.  
  1786.         } /* if( ntp <= 1) */
  1787.  
  1788.         data.t2x[mi] = ynv * data.t1z[mi] - znv * data.t1y[mi];
  1789.         data.t2y[mi] = znv * data.t1x[mi] - xnv * data.t1z[mi];
  1790.         data.t2z[mi] = xnv * data.t1y[mi] - ynv * data.t1x[mi];
  1791.         data.psalp[mi] = 1.;
  1792.  
  1793.         if (nx != 0)
  1794.         {
  1795.                 data.m += nx * ny - 1;
  1796.  
  1797.                 /* Reallocate patch buffers */
  1798.                 mreq = data.m * sizeof (double);
  1799.                 mem_realloc ((void *) &data.px, mreq);
  1800.                 mem_realloc ((void *) &data.py, mreq);
  1801.                 mem_realloc ((void *) &data.pz, mreq);
  1802.                 mem_realloc ((void *) &data.t1x, mreq);
  1803.                 mem_realloc ((void *) &data.t1y, mreq);
  1804.                 mem_realloc ((void *) &data.t1z, mreq);
  1805.                 mem_realloc ((void *) &data.t2x, mreq);
  1806.                 mem_realloc ((void *) &data.t2y, mreq);
  1807.                 mem_realloc ((void *) &data.t2z, mreq);
  1808.                 mem_realloc ((void *) &data.pbi, mreq);
  1809.                 mem_realloc ((void *) &data.psalp, mreq);
  1810.  
  1811.                 xn2 = data.px[mi] - s1x - s2x;
  1812.                 yn2 = data.py[mi] - s1y - s2y;
  1813.                 zn2 = data.pz[mi] - s1z - s2z;
  1814.                 xs = data.t1x[mi];
  1815.                 ys = data.t1y[mi];
  1816.                 zs = data.t1z[mi];
  1817.                 xt = data.t2x[mi];
  1818.                 yt = data.t2y[mi];
  1819.                 zt = data.t2z[mi];
  1820.  
  1821.                 for (iy = 0; iy < ny; iy++)
  1822.                 {
  1823.                         xn2 += s2x;
  1824.                         yn2 += s2y;
  1825.                         zn2 += s2z;
  1826.  
  1827.                         for (ix = 1; ix <= nx; ix++)
  1828.                         {
  1829.                                 xst = (double) ix;
  1830.                                 data.px[mi] = xn2 + xst * s1x;
  1831.                                 data.py[mi] = yn2 + xst * s1y;
  1832.                                 data.pz[mi] = zn2 + xst * s1z;
  1833.                                 data.pbi[mi] = xa;
  1834.                                 data.psalp[mi] = 1.;
  1835.                                 data.t1x[mi] = xs;
  1836.                                 data.t1y[mi] = ys;
  1837.                                 data.t1z[mi] = zs;
  1838.                                 data.t2x[mi] = xt;
  1839.                                 data.t2y[mi] = yt;
  1840.                                 data.t2z[mi] = zt;
  1841.                                 mi++;
  1842.                         } /* for( ix = 0; ix < nx; ix++ ) */
  1843.  
  1844.                 } /* for( iy = 0; iy < ny; iy++ ) */
  1845.  
  1846.         } /* if( nx != 0) */
  1847.  
  1848.         data.ipsym = 0;
  1849.         data.np = data.n;
  1850.         data.mp = data.m;
  1851.  
  1852.         return;
  1853. }
  1854.  
  1855. /*-----------------------------------------------------------------------*/
  1856.  
  1857. /*** this function was an 'entry point' (part of) 'patch()' ***/
  1858. void subph (int nx, int ny)
  1859. {
  1860.         int mia, ix, iy, mi, mreq;
  1861.         double xs, ys, zs, xa, xst, s1x, s1y, s1z, s2x, s2y, s2z, saln, xt, yt;
  1862.  
  1863.         /* Reallocate patch buffers */
  1864.         if (ny == 0)
  1865.                 data.m += 3;
  1866.         else
  1867.                 data.m += 4;
  1868.  
  1869.         mreq = data.m * sizeof (double);
  1870.         mem_realloc ((void *) &data.px, mreq);
  1871.         mem_realloc ((void *) &data.py, mreq);
  1872.         mem_realloc ((void *) &data.pz, mreq);
  1873.         mem_realloc ((void *) &data.t1x, mreq);
  1874.         mem_realloc ((void *) &data.t1y, mreq);
  1875.         mem_realloc ((void *) &data.t1z, mreq);
  1876.         mem_realloc ((void *) &data.t2x, mreq);
  1877.         mem_realloc ((void *) &data.t2y, mreq);
  1878.         mem_realloc ((void *) &data.t2z, mreq);
  1879.         mem_realloc ((void *) &data.pbi, mreq);
  1880.         mem_realloc ((void *) &data.psalp, mreq);
  1881.         mem_realloc ((void *) &data.icon1, (data.n + data.m) * sizeof (int));
  1882.         mem_realloc ((void *) &data.icon2, (data.n + data.m) * sizeof (int));
  1883.  
  1884.         /* Shift patches to make room for new ones */
  1885.         if ((ny == 0) && (nx != data.m))
  1886.         {
  1887.                 for (iy = data.m - 1; iy > nx + 2; iy--)
  1888.                 {
  1889.                         ix = iy - 3;
  1890.                         data.px[iy] = data.px[ix];
  1891.                         data.py[iy] = data.py[ix];
  1892.                         data.pz[iy] = data.pz[ix];
  1893.                         data.pbi[iy] = data.pbi[ix];
  1894.                         data.psalp[iy] = data.psalp[ix];
  1895.                         data.t1x[iy] = data.t1x[ix];
  1896.                         data.t1y[iy] = data.t1y[ix];
  1897.                         data.t1z[iy] = data.t1z[ix];
  1898.                         data.t2x[iy] = data.t2x[ix];
  1899.                         data.t2y[iy] = data.t2y[ix];
  1900.                         data.t2z[iy] = data.t2z[ix];
  1901.                 }
  1902.  
  1903.         } /* if( (ny == 0) || (nx != m) ) */
  1904.  
  1905.         /* divide patch for connection */
  1906.         mi = nx - 1;
  1907.         xs = data.px[mi];
  1908.         ys = data.py[mi];
  1909.         zs = data.pz[mi];
  1910.         xa = data.pbi[mi] / 4.;
  1911.         xst = sqrtl (xa) / 2.;
  1912.         s1x = data.t1x[mi];
  1913.         s1y = data.t1y[mi];
  1914.         s1z = data.t1z[mi];
  1915.         s2x = data.t2x[mi];
  1916.         s2y = data.t2y[mi];
  1917.         s2z = data.t2z[mi];
  1918.         saln = data.psalp[mi];
  1919.         xt = xst;
  1920.         yt = xst;
  1921.  
  1922.         if (ny == 0)
  1923.                 mia = mi;
  1924.         else
  1925.         {
  1926.                 data.mp++;
  1927.                 mia = data.m - 1;
  1928.         }
  1929.  
  1930.         for (ix = 1; ix <= 4; ix++)
  1931.         {
  1932.                 data.px[mia] = xs + xt * s1x + yt * s2x;
  1933.                 data.py[mia] = ys + xt * s1y + yt * s2y;
  1934.                 data.pz[mia] = zs + xt * s1z + yt * s2z;
  1935.                 data.pbi[mia] = xa;
  1936.                 data.t1x[mia] = s1x;
  1937.                 data.t1y[mia] = s1y;
  1938.                 data.t1z[mia] = s1z;
  1939.                 data.t2x[mia] = s2x;
  1940.                 data.t2y[mia] = s2y;
  1941.                 data.t2z[mia] = s2z;
  1942.                 data.psalp[mia] = saln;
  1943.  
  1944.                 if (ix == 2)
  1945.                         yt = -yt;
  1946.  
  1947.                 if ((ix == 1) || (ix == 3))
  1948.                         xt = -xt;
  1949.  
  1950.                 mia++;
  1951.         }
  1952.  
  1953.         if (nx <= data.mp)
  1954.                 data.mp += 3;
  1955.  
  1956.         if (ny > 0)
  1957.                 data.pz[mi] = 10000.;
  1958.  
  1959.         return;
  1960. }
  1961.  
  1962. /*-----------------------------------------------------------------------*/
  1963.  
  1964. void readgm (
  1965.     char *gm,
  1966.     int *i1,
  1967.     int *i2,
  1968.     double *x1,
  1969.     double *y1,
  1970.     double *z1,
  1971.     double *x2,
  1972.     double *y2,
  1973.     double *z2,
  1974.     double *rad)
  1975. {
  1976.         char line_buf[134];
  1977.         int nlin, i, line_idx;
  1978.         int nint = 2, nflt = 7;
  1979.         int iarr[2] = {0, 0};
  1980.         double rarr[7] = {0., 0., 0., 0., 0., 0., 0.};
  1981.  
  1982.         /* read a line from input file */
  1983.         load_line (line_buf, input_fp);
  1984.         {
  1985.                 gm[0] = 'E';
  1986.                 gm[1] = 'N';
  1987.                 gm[2] = '\0';
  1988.                 nlin = 2;
  1989.         }
  1990.  
  1991.         /* get line length */
  1992.         nlin = strlen (line_buf);
  1993.  
  1994.         /* abort if card's mnemonic too short or missing */
  1995.         if (nlin < 2)
  1996.         {
  1997.                 fprintf (
  1998.                     output_fp,
  1999.                     "\n  GEOMETRY DATA CARD ERROR:"
  2000.                     "\n  CARD'S MNEMONIC CODE TOO SHORT OR MISSING.");
  2001.                 stop (-1);
  2002.         }
  2003.  
  2004.         /* extract card's mnemonic code */
  2005.         /* was strncpy (gm, line_buf, 2) ; */
  2006.         gm[0] = line_buf[0];
  2007.         gm[1] = line_buf[1];
  2008.         gm[2] = '\0';
  2009.  
  2010.         /* Exit if "XT" command read (for testing) */
  2011.         if (strcmp (gm, "XT") == 0)
  2012.         {
  2013.                 fprintf (stderr, "\nnec2c: Exiting after an \"XT\" command in readgm()\n");
  2014.                 fprintf (
  2015.                     output_fp, "\n\n  nec2c: Exiting after an \"XT\" command in readgm()");
  2016.                 stop (0);
  2017.         }
  2018.  
  2019.         /* Return if only mnemonic on card */
  2020.         if (nlin == 2)
  2021.         {
  2022.                 *i1 = *i2 = 0;
  2023.                 *x1 = *y1 = *z1 = *x2 = *y2 = *z2 = *rad = 0.;
  2024.                 return;
  2025.         }
  2026.  
  2027.         /* read integers from line */
  2028.         line_idx = 1;
  2029.         for (i = 0; i < nint; i++)
  2030.         {
  2031.                 /* Find first numerical character */
  2032.                 while (((line_buf[++line_idx] < '0') || (line_buf[line_idx] > '9')) &&
  2033.                        (line_buf[line_idx] != '+') && (line_buf[line_idx] != '-'))
  2034.                         if ((line_buf[line_idx] == '\0'))
  2035.                         {
  2036.                                 *i1 = iarr[0];
  2037.                                 *i2 = iarr[1];
  2038.                                 *x1 = rarr[0];
  2039.                                 *y1 = rarr[1];
  2040.                                 *z1 = rarr[2];
  2041.                                 *x2 = rarr[3];
  2042.                                 *y2 = rarr[4];
  2043.                                 *z2 = rarr[5];
  2044.                                 *rad = rarr[6];
  2045.                                 return;
  2046.                         }
  2047.  
  2048.                 /* read an integer from line */
  2049.                 iarr[i] = atoi (&line_buf[line_idx]);
  2050.  
  2051.                 /* traverse numerical field to next ' ' or ',' or '\0' */
  2052.                 line_idx--;
  2053.                 while ((line_buf[++line_idx] != ' ') && (line_buf[line_idx] != '        ') &&
  2054.                        (line_buf[line_idx] != ',') && (line_buf[line_idx] != '\0'))
  2055.                 {
  2056.                         /* test for non-numerical characters */
  2057.                         if (((line_buf[line_idx] < '0') || (line_buf[line_idx] > '9')) &&
  2058.                             (line_buf[line_idx] != '+') && (line_buf[line_idx] != '-'))
  2059.                         {
  2060.                                 fprintf (
  2061.                                     output_fp,
  2062.                                     "\n  GEOMETRY DATA CARD \"%s\" ERROR:"
  2063.                                     "\n  NON-NUMERICAL CHARACTER '%c' IN INTEGER FIELD AT "
  2064.                                     "CHAR. %d\n",
  2065.                                     gm,
  2066.                                     line_buf[line_idx],
  2067.                                     (line_idx + 1));
  2068.                                 stop (-1);
  2069.                         }
  2070.  
  2071.                 } /* while( (line_buff[++line_idx] ... */
  2072.  
  2073.                 /* Return on end of line */
  2074.                 if (line_buf[line_idx] == '\0')
  2075.                 {
  2076.                         *i1 = iarr[0];
  2077.                         *i2 = iarr[1];
  2078.                         *x1 = rarr[0];
  2079.                         *y1 = rarr[1];
  2080.                         *z1 = rarr[2];
  2081.                         *x2 = rarr[3];
  2082.                         *y2 = rarr[4];
  2083.                         *z2 = rarr[5];
  2084.                         *rad = rarr[6];
  2085.                         return;
  2086.                 }
  2087.  
  2088.         } /* for( i = 0; i < nint; i++ ) */
  2089.  
  2090.         /* read doubles from line */
  2091.         for (i = 0; i < nflt; i++)
  2092.         {
  2093.                 /* Find first numerical character */
  2094.                 while (((line_buf[++line_idx] < '0') || (line_buf[line_idx] > '9')) &&
  2095.                        (line_buf[line_idx] != '+') && (line_buf[line_idx] != '-') &&
  2096.                        (line_buf[line_idx] != '.'))
  2097.                         if ((line_buf[line_idx] == '\0'))
  2098.                         {
  2099.                                 *i1 = iarr[0];
  2100.                                 *i2 = iarr[1];
  2101.                                 *x1 = rarr[0];
  2102.                                 *y1 = rarr[1];
  2103.                                 *z1 = rarr[2];
  2104.                                 *x2 = rarr[3];
  2105.                                 *y2 = rarr[4];
  2106.                                 *z2 = rarr[5];
  2107.                                 *rad = rarr[6];
  2108.                                 return;
  2109.                         }
  2110.  
  2111.                 /* read a double from line */
  2112.                 rarr[i] = atof (&line_buf[line_idx]);
  2113.  
  2114.                 /* traverse numerical field to next ' ' or ',' or '\0' */
  2115.                 line_idx--;
  2116.                 while ((line_buf[++line_idx] != ' ') && (line_buf[line_idx] != '        ') &&
  2117.                        (line_buf[line_idx] != ',') && (line_buf[line_idx] != '\0'))
  2118.                 {
  2119.                         /* test for non-numerical characters */
  2120.                         if (((line_buf[line_idx] < '0') || (line_buf[line_idx] > '9')) &&
  2121.                             (line_buf[line_idx] != '.') && (line_buf[line_idx] != '+') &&
  2122.                             (line_buf[line_idx] != '-') && (line_buf[line_idx] != 'E') &&
  2123.                             (line_buf[line_idx] != 'e'))
  2124.                         {
  2125.                                 fprintf (
  2126.                                     output_fp,
  2127.                                     "\n  GEOMETRY DATA CARD \"%s\" ERROR:"
  2128.                                     "\n  NON-NUMERICAL CHARACTER '%c' IN FLOAT FIELD AT CHAR. "
  2129.                                     "%d.\n",
  2130.                                     gm,
  2131.                                     line_buf[line_idx],
  2132.                                     (line_idx + 1));
  2133.                                 stop (-1);
  2134.                         }
  2135.  
  2136.                 } /* while( (line_buff[++line_idx] ... */
  2137.  
  2138.                 /* Return on end of line */
  2139.                 if (line_buf[line_idx] == '\0')
  2140.                 {
  2141.                         *i1 = iarr[0];
  2142.                         *i2 = iarr[1];
  2143.                         *x1 = rarr[0];
  2144.                         *y1 = rarr[1];
  2145.                         *z1 = rarr[2];
  2146.                         *x2 = rarr[3];
  2147.                         *y2 = rarr[4];
  2148.                         *z2 = rarr[5];
  2149.                         *rad = rarr[6];
  2150.                         return;
  2151.                 }
  2152.  
  2153.         } /* for( i = 0; i < nflt; i++ ) */
  2154.  
  2155.         *i1 = iarr[0];
  2156.         *i2 = iarr[1];
  2157.         *x1 = rarr[0];
  2158.         *y1 = rarr[1];
  2159.         *z1 = rarr[2];
  2160.         *x2 = rarr[3];
  2161.         *y2 = rarr[4];
  2162.         *z2 = rarr[5];
  2163.         *rad = rarr[6];
  2164.  
  2165.         return;
  2166. }
  2167.  
  2168. /*-----------------------------------------------------------------------*/
  2169.  
  2170. /* reflc reflects partial structure along x,y, or z axes or rotates */
  2171. /* structure to complete a symmetric structure. */
  2172. void reflc (int ix, int iy, int iz, int itx, int nop)
  2173. {
  2174.         int iti, i, nx, itagi, k, mreq;
  2175.         double e1, e2, fnop, sam, cs, ss, xk, yk;
  2176.  
  2177.         data.np = data.n;
  2178.         data.mp = data.m;
  2179.         data.ipsym = 0;
  2180.         iti = itx;
  2181.  
  2182.         if (ix >= 0)
  2183.         {
  2184.                 if (nop == 0)
  2185.                         return;
  2186.  
  2187.                 data.ipsym = 1;
  2188.  
  2189.                 /* reflect along z axis */
  2190.                 if (iz != 0)
  2191.                 {
  2192.                         data.ipsym = 2;
  2193.  
  2194.                         if (data.n > 0)
  2195.                         {
  2196.                                 /* Reallocate tags buffer */
  2197.                                 mem_realloc (
  2198.                                     (void *) &data.itag, (2 * data.n + data.m) * sizeof (int));
  2199.  
  2200.                                 /* Reallocate wire buffers */
  2201.                                 mreq = 2 * data.n * sizeof (double);
  2202.                                 mem_realloc ((void *) &data.x1, mreq);
  2203.                                 mem_realloc ((void *) &data.y1, mreq);
  2204.                                 mem_realloc ((void *) &data.z1, mreq);
  2205.                                 mem_realloc ((void *) &data.x2, mreq);
  2206.                                 mem_realloc ((void *) &data.y2, mreq);
  2207.                                 mem_realloc ((void *) &data.z2, mreq);
  2208.                                 mem_realloc ((void *) &data.bi, mreq);
  2209.  
  2210.                                 for (i = 0; i < data.n; i++)
  2211.                                 {
  2212.                                         nx = i + data.n;
  2213.                                         e1 = data.z1[i];
  2214.                                         e2 = data.z2[i];
  2215.  
  2216.                                         if ((fabsl (e1) + fabsl (e2) <= 1.0e-5) ||
  2217.                                             (e1 * e2 < -1.0e-6))
  2218.                                         {
  2219.                                                 fprintf (
  2220.                                                     output_fp,
  2221.                                                     "\n  GEOMETRY DATA ERROR--SEGMENT %d"
  2222.                                                     " LIES IN PLANE OF SYMMETRY",
  2223.                                                     i + 1);
  2224.                                                 stop (-1);
  2225.                                         }
  2226.  
  2227.                                         data.x1[nx] = data.x1[i];
  2228.                                         data.y1[nx] = data.y1[i];
  2229.                                         data.z1[nx] = -e1;
  2230.                                         data.x2[nx] = data.x2[i];
  2231.                                         data.y2[nx] = data.y2[i];
  2232.                                         data.z2[nx] = -e2;
  2233.                                         itagi = data.itag[i];
  2234.  
  2235.                                         if (itagi == 0)
  2236.                                                 data.itag[nx] = 0;
  2237.                                         if (itagi != 0)
  2238.                                                 data.itag[nx] = itagi + iti;
  2239.  
  2240.                                         data.bi[nx] = data.bi[i];
  2241.  
  2242.                                 } /* for( i = 0; i < data.n; i++ ) */
  2243.  
  2244.                                 data.n = data.n * 2;
  2245.                                 iti = iti * 2;
  2246.  
  2247.                         } /* if( data.n > 0) */
  2248.  
  2249.                         if (data.m > 0)
  2250.                         {
  2251.                                 /* Reallocate patch buffers */
  2252.                                 mreq = 2 * data.m * sizeof (double);
  2253.                                 mem_realloc ((void *) &data.px, mreq);
  2254.                                 mem_realloc ((void *) &data.py, mreq);
  2255.                                 mem_realloc ((void *) &data.pz, mreq);
  2256.                                 mem_realloc ((void *) &data.t1x, mreq);
  2257.                                 mem_realloc ((void *) &data.t1y, mreq);
  2258.                                 mem_realloc ((void *) &data.t1z, mreq);
  2259.                                 mem_realloc ((void *) &data.t2x, mreq);
  2260.                                 mem_realloc ((void *) &data.t2y, mreq);
  2261.                                 mem_realloc ((void *) &data.t2z, mreq);
  2262.                                 mem_realloc ((void *) &data.pbi, mreq);
  2263.                                 mem_realloc ((void *) &data.psalp, mreq);
  2264.  
  2265.                                 for (i = 0; i < data.m; i++)
  2266.                                 {
  2267.                                         nx = i + data.m;
  2268.                                         if (fabsl (data.pz[i]) <= 1.0e-10)
  2269.                                         {
  2270.                                                 fprintf (
  2271.                                                     output_fp,
  2272.                                                     "\n  GEOMETRY DATA ERROR--PATCH %d"
  2273.                                                     " LIES IN PLANE OF SYMMETRY",
  2274.                                                     i + 1);
  2275.                                                 stop (-1);
  2276.                                         }
  2277.  
  2278.                                         data.px[nx] = data.px[i];
  2279.                                         data.py[nx] = data.py[i];
  2280.                                         data.pz[nx] = -data.pz[i];
  2281.                                         data.t1x[nx] = data.t1x[i];
  2282.                                         data.t1y[nx] = data.t1y[i];
  2283.                                         data.t1z[nx] = -data.t1z[i];
  2284.                                         data.t2x[nx] = data.t2x[i];
  2285.                                         data.t2y[nx] = data.t2y[i];
  2286.                                         data.t2z[nx] = -data.t2z[i];
  2287.                                         data.psalp[nx] = -data.psalp[i];
  2288.                                         data.pbi[nx] = data.pbi[i];
  2289.                                 }
  2290.  
  2291.                                 data.m = data.m * 2;
  2292.  
  2293.                         } /* if( data.m >= m2) */
  2294.  
  2295.                 } /* if( iz != 0) */
  2296.  
  2297.                 /* reflect along y axis */
  2298.                 if (iy != 0)
  2299.                 {
  2300.                         if (data.n > 0)
  2301.                         {
  2302.                                 /* Reallocate tags buffer */
  2303.                                 mem_realloc (
  2304.                                     (void *) &data.itag,
  2305.                                     (2 * data.n + data.m) * sizeof (int)); /*????*/
  2306.  
  2307.                                 /* Reallocate wire buffers */
  2308.                                 mreq = 2 * data.n * sizeof (double);
  2309.                                 mem_realloc ((void *) &data.x1, mreq);
  2310.                                 mem_realloc ((void *) &data.y1, mreq);
  2311.                                 mem_realloc ((void *) &data.z1, mreq);
  2312.                                 mem_realloc ((void *) &data.x2, mreq);
  2313.                                 mem_realloc ((void *) &data.y2, mreq);
  2314.                                 mem_realloc ((void *) &data.z2, mreq);
  2315.                                 mem_realloc ((void *) &data.bi, mreq);
  2316.  
  2317.                                 for (i = 0; i < data.n; i++)
  2318.                                 {
  2319.                                         nx = i + data.n;
  2320.                                         e1 = data.y1[i];
  2321.                                         e2 = data.y2[i];
  2322.  
  2323.                                         if ((fabsl (e1) + fabsl (e2) <= 1.0e-5) ||
  2324.                                             (e1 * e2 < -1.0e-6))
  2325.                                         {
  2326.                                                 fprintf (
  2327.                                                     output_fp,
  2328.                                                     "\n  GEOMETRY DATA ERROR--SEGMENT %d"
  2329.                                                     " LIES IN PLANE OF SYMMETRY",
  2330.                                                     i + 1);
  2331.                                                 stop (-1);
  2332.                                         }
  2333.  
  2334.                                         data.x1[nx] = data.x1[i];
  2335.                                         data.y1[nx] = -e1;
  2336.                                         data.z1[nx] = data.z1[i];
  2337.                                         data.x2[nx] = data.x2[i];
  2338.                                         data.y2[nx] = -e2;
  2339.                                         data.z2[nx] = data.z2[i];
  2340.                                         itagi = data.itag[i];
  2341.  
  2342.                                         if (itagi == 0)
  2343.                                                 data.itag[nx] = 0;
  2344.                                         if (itagi != 0)
  2345.                                                 data.itag[nx] = itagi + iti;
  2346.  
  2347.                                         data.bi[nx] = data.bi[i];
  2348.  
  2349.                                 } /* for( i = n2-1; i < data.n; i++ ) */
  2350.  
  2351.                                 data.n = data.n * 2;
  2352.                                 iti = iti * 2;
  2353.  
  2354.                         } /* if( data.n >= n2) */
  2355.  
  2356.                         if (data.m > 0)
  2357.                         {
  2358.                                 /* Reallocate patch buffers */
  2359.                                 mreq = 2 * data.m * sizeof (double);
  2360.                                 mem_realloc ((void *) &data.px, mreq);
  2361.                                 mem_realloc ((void *) &data.py, mreq);
  2362.                                 mem_realloc ((void *) &data.pz, mreq);
  2363.                                 mem_realloc ((void *) &data.t1x, mreq);
  2364.                                 mem_realloc ((void *) &data.t1y, mreq);
  2365.                                 mem_realloc ((void *) &data.t1z, mreq);
  2366.                                 mem_realloc ((void *) &data.t2x, mreq);
  2367.                                 mem_realloc ((void *) &data.t2y, mreq);
  2368.                                 mem_realloc ((void *) &data.t2z, mreq);
  2369.                                 mem_realloc ((void *) &data.pbi, mreq);
  2370.                                 mem_realloc ((void *) &data.psalp, mreq);
  2371.  
  2372.                                 for (i = 0; i < data.m; i++)
  2373.                                 {
  2374.                                         nx = i + data.m;
  2375.                                         if (fabsl (data.py[i]) <= 1.0e-10)
  2376.                                         {
  2377.                                                 fprintf (
  2378.                                                     output_fp,
  2379.                                                     "\n  GEOMETRY DATA ERROR--PATCH %d"
  2380.                                                     " LIES IN PLANE OF SYMMETRY",
  2381.                                                     i + 1);
  2382.                                                 stop (-1);
  2383.                                         }
  2384.  
  2385.                                         data.px[nx] = data.px[i];
  2386.                                         data.py[nx] = -data.py[i];
  2387.                                         data.pz[nx] = data.pz[i];
  2388.                                         data.t1x[nx] = data.t1x[i];
  2389.                                         data.t1y[nx] = -data.t1y[i];
  2390.                                         data.t1z[nx] = data.t1z[i];
  2391.                                         data.t2x[nx] = data.t2x[i];
  2392.                                         data.t2y[nx] = -data.t2y[i];
  2393.                                         data.t2z[nx] = data.t2z[i];
  2394.                                         data.psalp[nx] = -data.psalp[i];
  2395.                                         data.pbi[nx] = data.pbi[i];
  2396.  
  2397.                                 } /* for( i = m2; i <= data.m; i++ ) */
  2398.  
  2399.                                 data.m = data.m * 2;
  2400.  
  2401.                         } /* if( data.m >= m2) */
  2402.  
  2403.                 } /* if( iy != 0) */
  2404.  
  2405.                 /* reflect along x axis */
  2406.                 if (ix == 0)
  2407.                         return;
  2408.  
  2409.                 if (data.n > 0)
  2410.                 {
  2411.                         /* Reallocate tags buffer */
  2412.                         mem_realloc (
  2413.                             (void *) &data.itag,
  2414.                             (2 * data.n + data.m) * sizeof (int)); /*????*/
  2415.  
  2416.                         /* Reallocate wire buffers */
  2417.                         mreq = 2 * data.n * sizeof (double);
  2418.                         mem_realloc ((void *) &data.x1, mreq);
  2419.                         mem_realloc ((void *) &data.y1, mreq);
  2420.                         mem_realloc ((void *) &data.z1, mreq);
  2421.                         mem_realloc ((void *) &data.x2, mreq);
  2422.                         mem_realloc ((void *) &data.y2, mreq);
  2423.                         mem_realloc ((void *) &data.z2, mreq);
  2424.                         mem_realloc ((void *) &data.bi, mreq);
  2425.  
  2426.                         for (i = 0; i < data.n; i++)
  2427.                         {
  2428.                                 nx = i + data.n;
  2429.                                 e1 = data.x1[i];
  2430.                                 e2 = data.x2[i];
  2431.  
  2432.                                 if ((fabsl (e1) + fabsl (e2) <= 1.0e-5) || (e1 * e2 < -1.0e-6))
  2433.                                 {
  2434.                                         fprintf (
  2435.                                             output_fp,
  2436.                                             "\n  GEOMETRY DATA ERROR--SEGMENT %d"
  2437.                                             " LIES IN PLANE OF SYMMETRY",
  2438.                                             i + 1);
  2439.                                         stop (-1);
  2440.                                 }
  2441.  
  2442.                                 data.x1[nx] = -e1;
  2443.                                 data.y1[nx] = data.y1[i];
  2444.                                 data.z1[nx] = data.z1[i];
  2445.                                 data.x2[nx] = -e2;
  2446.                                 data.y2[nx] = data.y2[i];
  2447.                                 data.z2[nx] = data.z2[i];
  2448.                                 itagi = data.itag[i];
  2449.  
  2450.                                 if (itagi == 0)
  2451.                                         data.itag[nx] = 0;
  2452.                                 if (itagi != 0)
  2453.                                         data.itag[nx] = itagi + iti;
  2454.  
  2455.                                 data.bi[nx] = data.bi[i];
  2456.                         }
  2457.  
  2458.                         data.n = data.n * 2;
  2459.  
  2460.                 } /* if( data.n > 0) */
  2461.  
  2462.                 if (data.m == 0)
  2463.                         return;
  2464.  
  2465.                 /* Reallocate patch buffers */
  2466.                 mreq = 2 * data.m * sizeof (double);
  2467.                 mem_realloc ((void *) &data.px, mreq);
  2468.                 mem_realloc ((void *) &data.py, mreq);
  2469.                 mem_realloc ((void *) &data.pz, mreq);
  2470.                 mem_realloc ((void *) &data.t1x, mreq);
  2471.                 mem_realloc ((void *) &data.t1y, mreq);
  2472.                 mem_realloc ((void *) &data.t1z, mreq);
  2473.                 mem_realloc ((void *) &data.t2x, mreq);
  2474.                 mem_realloc ((void *) &data.t2y, mreq);
  2475.                 mem_realloc ((void *) &data.t2z, mreq);
  2476.                 mem_realloc ((void *) &data.pbi, mreq);
  2477.                 mem_realloc ((void *) &data.psalp, mreq);
  2478.  
  2479.                 for (i = 0; i < data.m; i++)
  2480.                 {
  2481.                         nx = i + data.m;
  2482.                         if (fabsl (data.px[i]) <= 1.0e-10)
  2483.                         {
  2484.                                 fprintf (
  2485.                                     output_fp,
  2486.                                     "\n  GEOMETRY DATA ERROR--PATCH %d"
  2487.                                     " LIES IN PLANE OF SYMMETRY",
  2488.                                     i + 1);
  2489.                                 stop (-1);
  2490.                         }
  2491.  
  2492.                         data.px[nx] = -data.px[i];
  2493.                         data.py[nx] = data.py[i];
  2494.                         data.pz[nx] = data.pz[i];
  2495.                         data.t1x[nx] = -data.t1x[i];
  2496.                         data.t1y[nx] = data.t1y[i];
  2497.                         data.t1z[nx] = data.t1z[i];
  2498.                         data.t2x[nx] = -data.t2x[i];
  2499.                         data.t2y[nx] = data.t2y[i];
  2500.                         data.t2z[nx] = data.t2z[i];
  2501.                         data.psalp[nx] = -data.psalp[i];
  2502.                         data.pbi[nx] = data.pbi[i];
  2503.                 }
  2504.  
  2505.                 data.m = data.m * 2;
  2506.                 return;
  2507.  
  2508.         } /* if( ix >= 0) */
  2509.  
  2510.         /* reproduce structure with rotation to form cylindrical structure */
  2511.         fnop = (double) nop;
  2512.         data.ipsym = -1;
  2513.         sam = TP / fnop;
  2514.         cs = cosl (sam);
  2515.         ss = sinl (sam);
  2516.  
  2517.         if (data.n > 0)
  2518.         {
  2519.                 data.n *= nop;
  2520.                 nx = data.np;
  2521.  
  2522.                 /* Reallocate tags buffer */
  2523.                 mem_realloc ((void *) &data.itag, (data.n + data.m) * sizeof (int)); /*????*/
  2524.  
  2525.                 /* Reallocate wire buffers */
  2526.                 mreq = data.n * sizeof (double);
  2527.                 mem_realloc ((void *) &data.x1, mreq);
  2528.                 mem_realloc ((void *) &data.y1, mreq);
  2529.                 mem_realloc ((void *) &data.z1, mreq);
  2530.                 mem_realloc ((void *) &data.x2, mreq);
  2531.                 mem_realloc ((void *) &data.y2, mreq);
  2532.                 mem_realloc ((void *) &data.z2, mreq);
  2533.                 mem_realloc ((void *) &data.bi, mreq);
  2534.  
  2535.                 for (i = nx; i < data.n; i++)
  2536.                 {
  2537.                         k = i - data.np;
  2538.                         xk = data.x1[k];
  2539.                         yk = data.y1[k];
  2540.                         data.x1[i] = xk * cs - yk * ss;
  2541.                         data.y1[i] = xk * ss + yk * cs;
  2542.                         data.z1[i] = data.z1[k];
  2543.                         xk = data.x2[k];
  2544.                         yk = data.y2[k];
  2545.                         data.x2[i] = xk * cs - yk * ss;
  2546.                         data.y2[i] = xk * ss + yk * cs;
  2547.                         data.z2[i] = data.z2[k];
  2548.                         data.bi[i] = data.bi[k];
  2549.                         itagi = data.itag[k];
  2550.  
  2551.                         if (itagi == 0)
  2552.                                 data.itag[i] = 0;
  2553.                         if (itagi != 0)
  2554.                                 data.itag[i] = itagi + iti;
  2555.                 }
  2556.  
  2557.         } /* if( data.n >= n2) */
  2558.  
  2559.         if (data.m == 0)
  2560.                 return;
  2561.  
  2562.         data.m *= nop;
  2563.         nx = data.mp;
  2564.  
  2565.         /* Reallocate patch buffers */
  2566.         mreq = data.m * sizeof (double);
  2567.         mem_realloc ((void *) &data.px, mreq);
  2568.         mem_realloc ((void *) &data.py, mreq);
  2569.         mem_realloc ((void *) &data.pz, mreq);
  2570.         mem_realloc ((void *) &data.t1x, mreq);
  2571.         mem_realloc ((void *) &data.t1y, mreq);
  2572.         mem_realloc ((void *) &data.t1z, mreq);
  2573.         mem_realloc ((void *) &data.t2x, mreq);
  2574.         mem_realloc ((void *) &data.t2y, mreq);
  2575.         mem_realloc ((void *) &data.t2z, mreq);
  2576.         mem_realloc ((void *) &data.pbi, mreq);
  2577.         mem_realloc ((void *) &data.psalp, mreq);
  2578.  
  2579.         for (i = nx; i < data.m; i++)
  2580.         {
  2581.                 k = i - data.mp;
  2582.                 xk = data.px[k];
  2583.                 yk = data.py[k];
  2584.                 data.px[i] = xk * cs - yk * ss;
  2585.                 data.py[i] = xk * ss + yk * cs;
  2586.                 data.pz[i] = data.pz[k];
  2587.                 xk = data.t1x[k];
  2588.                 yk = data.t1y[k];
  2589.                 data.t1x[i] = xk * cs - yk * ss;
  2590.                 data.t1y[i] = xk * ss + yk * cs;
  2591.                 data.t1z[i] = data.t1z[k];
  2592.                 xk = data.t2x[k];
  2593.                 yk = data.t2y[k];
  2594.                 data.t2x[i] = xk * cs - yk * ss;
  2595.                 data.t2y[i] = xk * ss + yk * cs;
  2596.                 data.t2z[i] = data.t2z[k];
  2597.                 data.psalp[i] = data.psalp[k];
  2598.                 data.pbi[i] = data.pbi[k];
  2599.  
  2600.         } /* for( i = nx; i < data.m; i++ ) */
  2601.  
  2602.         return;
  2603. }
  2604.  
  2605. /*-----------------------------------------------------------------------*/
  2606.  
  2607. /* subroutine wire generates segment geometry */
  2608. /* data for a straight wire of ns segments. */
  2609. void wire (
  2610.     double xw1,
  2611.     double yw1,
  2612.     double zw1,
  2613.     double xw2,
  2614.     double yw2,
  2615.     double zw2,
  2616.     double rad,
  2617.     double rdel,
  2618.     double rrad,
  2619.     int ns,
  2620.     int itg)
  2621. {
  2622.         int ist, i, mreq;
  2623.         double xd, yd, zd, delz, rd, fns, radz;
  2624.         double xs1, ys1, zs1, xs2, ys2, zs2;
  2625.  
  2626.         ist = data.n;
  2627.         data.n = data.n + ns;
  2628.         data.np = data.n;
  2629.         data.mp = data.m;
  2630.         data.ipsym = 0;
  2631.  
  2632.         if (ns < 1)
  2633.                 return;
  2634.  
  2635.         /* Reallocate tags buffer */
  2636.         mem_realloc ((void *) &data.itag, (data.n + data.m) * sizeof (int)); /*????*/
  2637.  
  2638.         /* Reallocate wire buffers */
  2639.         mreq = data.n * sizeof (double);
  2640.         mem_realloc ((void *) &data.x1, mreq);
  2641.         mem_realloc ((void *) &data.y1, mreq);
  2642.         mem_realloc ((void *) &data.z1, mreq);
  2643.         mem_realloc ((void *) &data.x2, mreq);
  2644.         mem_realloc ((void *) &data.y2, mreq);
  2645.         mem_realloc ((void *) &data.z2, mreq);
  2646.         mem_realloc ((void *) &data.bi, mreq);
  2647.  
  2648.         xd = xw2 - xw1;
  2649.         yd = yw2 - yw1;
  2650.         zd = zw2 - zw1;
  2651.  
  2652.         if (fabsl (rdel - 1.) >= 1.0e-6)
  2653.         {
  2654.                 delz = sqrtl (xd * xd + yd * yd + zd * zd);
  2655.                 xd = xd / delz;
  2656.                 yd = yd / delz;
  2657.                 zd = zd / delz;
  2658.                 delz = delz * (1. - rdel) / (1. - powl (rdel, ns));
  2659.                 rd = rdel;
  2660.         }
  2661.         else
  2662.         {
  2663.                 fns = ns;
  2664.                 xd = xd / fns;
  2665.                 yd = yd / fns;
  2666.                 zd = zd / fns;
  2667.                 delz = 1.;
  2668.                 rd = 1.;
  2669.         }
  2670.  
  2671.         radz = rad;
  2672.         xs1 = xw1;
  2673.         ys1 = yw1;
  2674.         zs1 = zw1;
  2675.  
  2676.         for (i = ist; i < data.n; i++)
  2677.         {
  2678.                 data.itag[i] = itg;
  2679.                 xs2 = xs1 + xd * delz;
  2680.                 ys2 = ys1 + yd * delz;
  2681.                 zs2 = zs1 + zd * delz;
  2682.                 data.x1[i] = xs1;
  2683.                 data.y1[i] = ys1;
  2684.                 data.z1[i] = zs1;
  2685.                 data.x2[i] = xs2;
  2686.                 data.y2[i] = ys2;
  2687.                 data.z2[i] = zs2;
  2688.                 data.bi[i] = radz;
  2689.                 delz = delz * rd;
  2690.                 radz = radz * rrad;
  2691.                 xs1 = xs2;
  2692.                 ys1 = ys2;
  2693.                 zs1 = zs2;
  2694.         }
  2695.  
  2696.         data.x2[data.n - 1] = xw2;
  2697.         data.y2[data.n - 1] = yw2;
  2698.         data.z2[data.n - 1] = zw2;
  2699.  
  2700.         return;
  2701. }
  2702.  
  2703. /*-----------------------------------------------------------------------*/
  2704.