home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / libs / gle / util / fitz / fit.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-29  |  65.5 KB  |  2,712 lines
  1. /* fit.f -- translated by f2c (version of 26 February 1990  17:38:00).
  2.    You must link the resulting object file with the libraries:
  3.     -lF77 -lI77 -lm -lc   (in that order)
  4. */
  5.  
  6. #include "f2c.h"
  7.  
  8. /* Common Block Declarations */
  9.  
  10. struct {
  11.     integer nit;
  12. } idlc_;
  13.  
  14. #define idlc_1 idlc_
  15.  
  16. struct {
  17.     integer itpv;
  18. } idpi_;
  19.  
  20. #define idpi_1 idpi_
  21.  
  22. /* Subroutine */ int idbvip_(md, ncp, ndp, xd, yd, zd, nip, xi, yi, zi, iwk, 
  23.     wk)
  24. integer *md, *ncp, *ndp;
  25. real *xd, *yd, *zd;
  26. integer *nip;
  27. real *xi, *yi, *zi;
  28. integer *iwk;
  29. real *wk;
  30. {
  31.  
  32.     /* System generated locals */
  33.     integer i_1, i_2;
  34.  
  35.     /* Local variables */
  36.     static integer jwit, jwit0;
  37.     extern /* Subroutine */ int err2090_();
  38.     static integer jwipc, jwipl, ncppv, ndppv, jwiwk, nippv, jwipt, jwiwl, 
  39.         jwiwp, nl;
  40.     extern /* Subroutine */ int idcldp_();
  41.     static integer nt;
  42.     extern /* Subroutine */ int idtang_(), idlctn_(), idpdrv_(), idptip_();
  43.     static integer md0, iip, ncp0, ndp0, nip0;
  44.  
  45. /* this subroutine performs bivariate interpolation when the pro- */
  46. /* jections of the data points in the x-y plane are irregularly */
  47. /* distributed in the plane. */
  48. /* the input parameters are */
  49. /*     md  = mode of computation (must be 1, 2, or 3), */
  50. /*         = 1 for new ncp and/or new xd-yd, */
  51. /*         = 2 for old ncp, old xd-yd, new xi-yi, */
  52. /*         = 3 for old ncp, old xd-yd, old xi-yi, */
  53. /*     ncp = number of additional data points used for esti- */
  54. /*           mating partial derivatives at each data point */
  55. /*           (must be 2 or greater, but smaller than ndp), */
  56. /*     ndp = number of data points (must be 4 or greater), */
  57. /*     xd  = array of dimension ndp containing the x */
  58. /*           coordinates of the data points, */
  59. /*     yd  = array of dimension ndp containing the y */
  60. /*           coordinates of the data points, */
  61. /*     zd  = array of dimension ndp containing the z */
  62. /*           coordinates of the data points, */
  63. /*     nip = number of output points at which interpolation */
  64. /*           is to be performed (must be 1 or greater), */
  65. /*     xi  = array of dimension nip containing the x */
  66. /*           coordinates of the output points, */
  67. /*     yi  = array of dimension nip containing the y */
  68. /*           coordinates of the output points. */
  69. /* the output parameter is */
  70. /*     zi  = array of dimension nip where interpolated z */
  71. /*           values are to be stored. */
  72. /* the other parameters are */
  73. /*     iwk = integer array of dimension */
  74. /*              max0(31,27+ncp)*ndp+nip */
  75. /*           used internally as a work area, */
  76. /*     wk  = array of dimension 8*ndp used internally as a */
  77. /*           work area. */
  78. /* the very first call to this subroutine and the call with a new */
  79. /* ncp value, a new ndp value, and/or new contents of the xd and */
  80. /* yd arrays must be made with md=1.  the call with md=2 must be */
  81. /* preceded by another call with the same ncp and ndp values and */
  82. /* with the same contents of the xd and yd arrays.  the call with */
  83. /* md=3 must be preceded by another call with the same ncp, ndp, */
  84. /* and nip values and with the same contents of the xd, yd, xi, */
  85. /* and yi arrays.  between the call with md=2 or md=3 and its */
  86. /* preceding call, the iwk and wk arrays must not be disturbed. */
  87. /* use of a value between 3 and 5 (inclusive) for ncp is recom- */
  88. /* mended unless there are evidences that dictate otherwise. */
  89. /* the lun constant in the data initialization statement is the */
  90. /* logical unit number of the standard output unit and is, */
  91. /* therefore, system dependent. */
  92. /* this subroutine calls the idcldp, idlctn, idpdrv, idptip, and */
  93. /* idtang subroutines. */
  94. /* declaration statements */
  95.     /* Parameter adjustments */
  96.     --wk;
  97.     --iwk;
  98.     --zi;
  99.     --yi;
  100.     --xi;
  101.     --zd;
  102.     --yd;
  103.     --xd;
  104.  
  105.     /* Function Body */
  106. /* setting of some input parameters to local variables. */
  107. /* (for md=1,2,3) */
  108. /* L10: */
  109.     md0 = *md;
  110.     ncp0 = *ncp;
  111.     ndp0 = *ndp;
  112.     nip0 = *nip;
  113. /* error check.  (for md=1,2,3) */
  114. /* L20: */
  115.     if (md0 < 1 || md0 > 3) {
  116.     goto L90;
  117.     }
  118.     if (ncp0 < 2 || ncp0 >= ndp0) {
  119.     goto L90;
  120.     }
  121.     if (ndp0 < 4) {
  122.     goto L90;
  123.     }
  124.     if (nip0 < 1) {
  125.     goto L90;
  126.     }
  127.     if (md0 >= 2) {
  128.     goto L21;
  129.     }
  130.     iwk[1] = ncp0;
  131.     iwk[2] = ndp0;
  132.     goto L22;
  133. L21:
  134.     ncppv = iwk[1];
  135.     ndppv = iwk[2];
  136.     if (ncp0 != ncppv) {
  137.     goto L90;
  138.     }
  139.     if (ndp0 != ndppv) {
  140.     goto L90;
  141.     }
  142. L22:
  143.     if (md0 >= 3) {
  144.     goto L23;
  145.     }
  146.     iwk[3] = *nip;
  147.     goto L30;
  148. L23:
  149.     nippv = iwk[3];
  150.     if (nip0 != nippv) {
  151.     goto L90;
  152.     }
  153. /* allocation of storage areas in the iwk array.  (for md=1,2,3) */
  154. L30:
  155.     jwipt = 16;
  156.     jwiwl = ndp0 * 6 + 1;
  157.     jwiwk = jwiwl;
  158.     jwipl = ndp0 * 24 + 1;
  159.     jwiwp = ndp0 * 30 + 1;
  160.     jwipc = ndp0 * 27 + 1;
  161. /* Computing MAX */
  162.     i_1 = 31, i_2 = ncp0 + 27;
  163.     jwit0 = max(i_1,i_2) * ndp0;
  164. /* triangulates the x-y plane.  (for md=1) */
  165. /* L40: */
  166.     if (md0 > 1) {
  167.     goto L50;
  168.     }
  169.     idtang_(&ndp0, &xd[1], &yd[1], &nt, &iwk[jwipt], &nl, &iwk[jwipl], &iwk[
  170.         jwiwl], &iwk[jwiwp], &wk[1]);
  171.     iwk[5] = nt;
  172.     iwk[6] = nl;
  173.     if (nt == 0) {
  174.     return 0;
  175.     }
  176. /* determines ncp points closest to each data point.  (for md=1) */
  177. L50:
  178.     if (md0 > 1) {
  179.     goto L60;
  180.     }
  181.     idcldp_(&ndp0, &xd[1], &yd[1], &ncp0, &iwk[jwipc]);
  182.     if (iwk[jwipc] == 0) {
  183.     return 0;
  184.     }
  185. /* locates all points at which interpolation is to be performed. */
  186. /* (for md=1,2) */
  187. L60:
  188.     if (md0 == 3) {
  189.     goto L70;
  190.     }
  191.     idlc_1.nit = 0;
  192.     jwit = jwit0;
  193.     i_1 = nip0;
  194.     for (iip = 1; iip <= i_1; ++iip) {
  195.     ++jwit;
  196.     idlctn_(&ndp0, &xd[1], &yd[1], &nt, &iwk[jwipt], &nl, &iwk[jwipl], &
  197.         xi[iip], &yi[iip], &iwk[jwit], &iwk[jwiwk], &wk[1]);
  198. /* L61: */
  199.     }
  200. /* estimates partial derivatives at all data points. */
  201. /* (for md=1,2,3) */
  202. L70:
  203.     idpdrv_(&ndp0, &xd[1], &yd[1], &zd[1], &ncp0, &iwk[jwipc], &wk[1]);
  204. /* interpolates the zi values.  (for md=1,2,3) */
  205. /* L80: */
  206.     idpi_1.itpv = 0;
  207.     jwit = jwit0;
  208.     i_1 = nip0;
  209.     for (iip = 1; iip <= i_1; ++iip) {
  210.     ++jwit;
  211.     idptip_(&xd[1], &yd[1], &zd[1], &nt, &iwk[jwipt], &nl, &iwk[jwipl], &
  212.         wk[1], &iwk[jwit], &xi[iip], &yi[iip], &zi[iip]);
  213. /* L81: */
  214.     }
  215.     return 0;
  216. /* error exit */
  217. L90:
  218.     err2090_();
  219.     return 0;
  220. /* format statement for error message */
  221. /* L2090: */
  222. } /* idbvip_ */
  223.  
  224. /* Subroutine */ int idcldp_(ndp, xd, yd, ncp, ipc)
  225. integer *ndp;
  226. real *xd, *yd;
  227. integer *ncp, *ipc;
  228. {
  229.     /* Initialized data */
  230.  
  231.     static integer ncpmx = 25;
  232.  
  233.     /* System generated locals */
  234.     integer i_1, i_2, i_3;
  235.     real r_1, r_2;
  236.  
  237.     /* Local variables */
  238.     static real dsqi;
  239.     static integer ip2mn, ip3mn;
  240.     extern /* Subroutine */ int err2090_(), err2091_(), err2092_();
  241.     static integer nclpt;
  242.     static real dsqmn;
  243.     static integer j1;
  244.     static real dsqmx;
  245.     static integer j3, j4, j2;
  246.     static real x1, y1;
  247.     static integer ip1, ip2, ip3;
  248.     static real dx12, dy12, dx13, dy13;
  249.     static integer jmx, ipc0[25], ncp0, ndp0;
  250.     static real dsq0[25];
  251.  
  252. /* this subroutine selects several data points that are closest */
  253. /* to each of the data point. */
  254. /* the input parameters are */
  255. /*     ndp = number of data points, */
  256. /*     xd,yd = arrays of dimension ndp containing the x and y */
  257. /*           coordinates of the data points, */
  258. /*     ncp = number of data points closest to each data */
  259. /*           points. */
  260. /* the output parameter is */
  261. /*     ipc = integer array of dimension ncp*ndp, where the */
  262. /*           point numbers of ncp data points closest to */
  263. /*           each of the ndp data points are to be stored. */
  264. /* this subroutine arbitrarily sets a restriction that ncp must */
  265. /* not exceed 25. */
  266. /* the lun constant in the data initialization statement is the */
  267. /* logical unit number of the standard output unit and is, */
  268. /* therefore, system dependent. */
  269. /* declaration statements */
  270.     /* Parameter adjustments */
  271.     --ipc;
  272.     --yd;
  273.     --xd;
  274.  
  275.     /* Function Body */
  276. /* statement function */
  277. /* preliminary processing */
  278. /* L10: */
  279.     ndp0 = *ndp;
  280.     ncp0 = *ncp;
  281.     if (ndp0 < 2) {
  282.     goto L90;
  283.     }
  284.     if (ncp0 < 1 || ncp0 > ncpmx || ncp0 >= ndp0) {
  285.     goto L90;
  286.     }
  287. /* calculation */
  288. /* L20: */
  289.     i_1 = ndp0;
  290.     for (ip1 = 1; ip1 <= i_1; ++ip1) {
  291. /* - selects ncp points. */
  292.     x1 = xd[ip1];
  293.     y1 = yd[ip1];
  294.     j1 = 0;
  295.     dsqmx = (float)0.;
  296.     i_2 = ndp0;
  297.     for (ip2 = 1; ip2 <= i_2; ++ip2) {
  298.         if (ip2 == ip1) {
  299.         goto L22;
  300.         }
  301. /* Computing 2nd power */
  302.         r_1 = xd[ip2] - x1;
  303. /* Computing 2nd power */
  304.         r_2 = yd[ip2] - y1;
  305.         dsqi = r_1 * r_1 + r_2 * r_2;
  306.         ++j1;
  307.         dsq0[j1 - 1] = dsqi;
  308.         ipc0[j1 - 1] = ip2;
  309.         if (dsqi <= dsqmx) {
  310.         goto L21;
  311.         }
  312.         dsqmx = dsqi;
  313.         jmx = j1;
  314. L21:
  315.         if (j1 >= ncp0) {
  316.         goto L23;
  317.         }
  318. L22:
  319.     ;}
  320. L23:
  321.     ip2mn = ip2 + 1;
  322.     if (ip2mn > ndp0) {
  323.         goto L30;
  324.     }
  325.     i_2 = ndp0;
  326.     for (ip2 = ip2mn; ip2 <= i_2; ++ip2) {
  327.         if (ip2 == ip1) {
  328.         goto L25;
  329.         }
  330. /* Computing 2nd power */
  331.         r_1 = xd[ip2] - x1;
  332. /* Computing 2nd power */
  333.         r_2 = yd[ip2] - y1;
  334.         dsqi = r_1 * r_1 + r_2 * r_2;
  335.         if (dsqi >= dsqmx) {
  336.         goto L25;
  337.         }
  338.         dsq0[jmx - 1] = dsqi;
  339.         ipc0[jmx - 1] = ip2;
  340.         dsqmx = (float)0.;
  341.         i_3 = ncp0;
  342.         for (j1 = 1; j1 <= i_3; ++j1) {
  343.         if (dsq0[j1 - 1] <= dsqmx) {
  344.             goto L24;
  345.         }
  346.         dsqmx = dsq0[j1 - 1];
  347.         jmx = j1;
  348. L24:
  349.         ;}
  350. L25:
  351.     ;}
  352. /* - checks if all the ncp+1 points are collinear. */
  353. L30:
  354.     ip2 = ipc0[0];
  355.     dx12 = xd[ip2] - x1;
  356.     dy12 = yd[ip2] - y1;
  357.     i_2 = ncp0;
  358.     for (j3 = 2; j3 <= i_2; ++j3) {
  359.         ip3 = ipc0[j3 - 1];
  360.         dx13 = xd[ip3] - x1;
  361.         dy13 = yd[ip3] - y1;
  362.         if (dy13 * dx12 - dx13 * dy12 != (float)0.) {
  363.         goto L50;
  364.         }
  365. /* L31: */
  366.     }
  367. /* - searches for the closest noncollinear point. */
  368. /* L40: */
  369.     nclpt = 0;
  370.     i_2 = ndp0;
  371.     for (ip3 = 1; ip3 <= i_2; ++ip3) {
  372.         if (ip3 == ip1) {
  373.         goto L43;
  374.         }
  375.         i_3 = ncp0;
  376.         for (j4 = 1; j4 <= i_3; ++j4) {
  377.         if (ip3 == ipc0[j4 - 1]) {
  378.             goto L43;
  379.         }
  380. /* L41: */
  381.         }
  382.         dx13 = xd[ip3] - x1;
  383.         dy13 = yd[ip3] - y1;
  384.         if (dy13 * dx12 - dx13 * dy12 == (float)0.) {
  385.         goto L43;
  386.         }
  387. /* Computing 2nd power */
  388.         r_1 = xd[ip3] - x1;
  389. /* Computing 2nd power */
  390.         r_2 = yd[ip3] - y1;
  391.         dsqi = r_1 * r_1 + r_2 * r_2;
  392.         if (nclpt == 0) {
  393.         goto L42;
  394.         }
  395.         if (dsqi >= dsqmn) {
  396.         goto L43;
  397.         }
  398. L42:
  399.         nclpt = 1;
  400.         dsqmn = dsqi;
  401.         ip3mn = ip3;
  402. L43:
  403.     ;}
  404.     if (nclpt == 0) {
  405.         goto L91;
  406.     }
  407.     dsqmx = dsqmn;
  408.     ipc0[jmx - 1] = ip3mn;
  409. /* - replaces the local array for the output array. */
  410. L50:
  411.     j1 = (ip1 - 1) * ncp0;
  412.     i_2 = ncp0;
  413.     for (j2 = 1; j2 <= i_2; ++j2) {
  414.         ++j1;
  415.         ipc[j1] = ipc0[j2 - 1];
  416. /* L51: */
  417.     }
  418. /* L59: */
  419.     }
  420.     return 0;
  421. /* error exit */
  422. L90:
  423.     err2090_();
  424.     goto L92;
  425. L91:
  426.     err2091_();
  427. L92:
  428.     err2092_();
  429.     ipc[1] = 0;
  430.     return 0;
  431. /* format statements for error messages */
  432. /* L2090: */
  433. /* L2091: */
  434. /* L2092: */
  435. } /* idcldp_ */
  436.  
  437. /* Subroutine */ int idgrid_(xd, yd, nt, ipt, nl, ipl, nxi, nyi, xi, yi, ngp, 
  438.     igp)
  439. real *xd, *yd;
  440. integer *nt, *ipt, *nl, *ipl, *nxi, *nyi;
  441. real *xi, *yi;
  442. integer *ngp, *igp;
  443. {
  444.     /* System generated locals */
  445.     integer i_1, i_2, i_3, i_4;
  446.     real r_1, r_2;
  447.  
  448.     /* Local variables */
  449.     static integer il0t3, insd, it0t3;
  450.     static real ximn, yimn, ximx, yimx;
  451.     static integer jigp0, jigp1, jngp0, jngp1, l, ilp1t3, iximn, iximx;
  452.     static real x1, y1, x2, y2, x3, y3;
  453.     static integer jigp1i, il0, nl0, ip1, ip2, it0, nxinyi, ip3, nt0, ixi, 
  454.         iyi;
  455.     static real yii, xii;
  456.     static integer izi;
  457.     static real xmn, ymn, xmx, ymx;
  458.     static integer ngp0, ngp1, ilp1, nxi0, nyi0;
  459.  
  460. /* this subroutine organizes grid points for surface fitting by */
  461. /* sorting them in ascending order of triangle numbers and of the */
  462. /* border line segment number. */
  463. /* the input parameters are */
  464. /*     xd,yd = arrays of dimension ndp containing the x and y */
  465. /*           coordinates of the data points, where ndp is the */
  466. /*           number of the data points, */
  467. /*     nt  = number of triangles, */
  468. /*     ipt = integer array of dimension 3*nt containing the */
  469. /*           point numbers of the vertexes of the triangles, */
  470. /*     nl  = number of border line segments, */
  471. /*     ipl = integer array of dimension 3*nl containing the */
  472. /*           point numbers of the end points of the border */
  473. /*           line segments and their respective triangle */
  474. /*           numbers, */
  475. /*     nxi = number of grid points in the x coordinate, */
  476. /*     nyi = number of grid points in the y coordinate, */
  477. /*     xi,yi = arrays of dimension nxi and nyi containing */
  478. /*           the x and y coordinates of the grid points, */
  479. /*           respectively. */
  480. /* the output parameters are */
  481. /*     ngp = integer array of dimension 2*(nt+2*nl) where the */
  482. /*           number of grid points that belong to each of the */
  483. /*           triangles or of the border line segments are to */
  484. /*           be stored, */
  485. /*     igp = integer array of dimension nxi*nyi where the */
  486. /*           grid point numbers are to be stored in ascending */
  487. /*           order of the triangle number and the border line */
  488. /*           segment number. */
  489. /* declaration statements */
  490. /* statement functions */
  491. /* preliminary processing */
  492.     /* Parameter adjustments */
  493.     --igp;
  494.     --ngp;
  495.     --yi;
  496.     --xi;
  497.     --ipl;
  498.     --ipt;
  499.     --yd;
  500.     --xd;
  501.  
  502.     /* Function Body */
  503.     nt0 = *nt;
  504.     nl0 = *nl;
  505.     nxi0 = *nxi;
  506.     nyi0 = *nyi;
  507.     nxinyi = nxi0 * nyi0;
  508. /* Computing MIN */
  509.     r_1 = xi[1], r_2 = xi[nxi0];
  510.     ximn = dmin(r_1,r_2);
  511. /* Computing MAX */
  512.     r_1 = xi[1], r_2 = xi[nxi0];
  513.     ximx = dmax(r_1,r_2);
  514. /* Computing MIN */
  515.     r_1 = yi[1], r_2 = yi[nyi0];
  516.     yimn = dmin(r_1,r_2);
  517. /* Computing MAX */
  518.     r_1 = yi[1], r_2 = yi[nyi0];
  519.     yimx = dmax(r_1,r_2);
  520. /* determines grid points inside the data area. */
  521.     jngp0 = 0;
  522.     jngp1 = (nt0 + (nl0 << 1) << 1) + 1;
  523.     jigp0 = 0;
  524.     jigp1 = nxinyi + 1;
  525.     i_1 = nt0;
  526.     for (it0 = 1; it0 <= i_1; ++it0) {
  527.     ngp0 = 0;
  528.     ngp1 = 0;
  529.     it0t3 = it0 * 3;
  530.     ip1 = ipt[it0t3 - 2];
  531.     ip2 = ipt[it0t3 - 1];
  532.     ip3 = ipt[it0t3];
  533.     x1 = xd[ip1];
  534.     y1 = yd[ip1];
  535.     x2 = xd[ip2];
  536.     y2 = yd[ip2];
  537.     x3 = xd[ip3];
  538.     y3 = yd[ip3];
  539. /* Computing MIN */
  540.     r_1 = min(x1,x2);
  541.     xmn = dmin(r_1,x3);
  542. /* Computing MAX */
  543.     r_1 = max(x1,x2);
  544.     xmx = dmax(r_1,x3);
  545. /* Computing MIN */
  546.     r_1 = min(y1,y2);
  547.     ymn = dmin(r_1,y3);
  548. /* Computing MAX */
  549.     r_1 = max(y1,y2);
  550.     ymx = dmax(r_1,y3);
  551.     insd = 0;
  552.     i_2 = nxi0;
  553.     for (ixi = 1; ixi <= i_2; ++ixi) {
  554.         if (xi[ixi] >= xmn && xi[ixi] <= xmx) {
  555.         goto L10;
  556.         }
  557.         if (insd == 0) {
  558.         goto L20;
  559.         }
  560.         iximx = ixi - 1;
  561.         goto L30;
  562. L10:
  563.         if (insd == 1) {
  564.         goto L20;
  565.         }
  566.         insd = 1;
  567.         iximn = ixi;
  568. L20:
  569.     ;}
  570.     if (insd == 0) {
  571.         goto L150;
  572.     }
  573.     iximx = nxi0;
  574. L30:
  575.     i_2 = nyi0;
  576.     for (iyi = 1; iyi <= i_2; ++iyi) {
  577.         yii = yi[iyi];
  578.         if (yii < ymn || yii > ymx) {
  579.         goto L140;
  580.         }
  581.         i_3 = iximx;
  582.         for (ixi = iximn; ixi <= i_3; ++ixi) {
  583.         xii = xi[ixi];
  584.         l = 0;
  585.         if ((r_1 = (x1 - xii) * (y2 - yii) - (y1 - yii) * (x2 - xii)) 
  586.             < (float)0.) {
  587.             goto L130;
  588.         } else if (r_1 == 0) {
  589.             goto L40;
  590.         } else {
  591.             goto L50;
  592.         }
  593. L40:
  594.         l = 1;
  595. L50:
  596.         if ((r_1 = (x2 - xii) * (y3 - yii) - (y2 - yii) * (x3 - xii)) 
  597.             < (float)0.) {
  598.             goto L130;
  599.         } else if (r_1 == 0) {
  600.             goto L60;
  601.         } else {
  602.             goto L70;
  603.         }
  604. L60:
  605.         l = 1;
  606. L70:
  607.         if ((r_1 = (x3 - xii) * (y1 - yii) - (y3 - yii) * (x1 - xii)) 
  608.             < (float)0.) {
  609.             goto L130;
  610.         } else if (r_1 == 0) {
  611.             goto L80;
  612.         } else {
  613.             goto L90;
  614.         }
  615. L80:
  616.         l = 1;
  617. L90:
  618.         izi = nxi0 * (iyi - 1) + ixi;
  619.         if (l == 1) {
  620.             goto L100;
  621.         }
  622.         ++ngp0;
  623.         ++jigp0;
  624.         igp[jigp0] = izi;
  625.         goto L130;
  626. L100:
  627.         if (jigp1 > nxinyi) {
  628.             goto L120;
  629.         }
  630.         i_4 = nxinyi;
  631.         for (jigp1i = jigp1; jigp1i <= i_4; ++jigp1i) {
  632.             if (izi == igp[jigp1i]) {
  633.             goto L130;
  634.             }
  635. /* L110: */
  636.         }
  637. L120:
  638.         ++ngp1;
  639.         --jigp1;
  640.         igp[jigp1] = izi;
  641. L130:
  642.         ;}
  643. L140:
  644.     ;}
  645. L150:
  646.     ++jngp0;
  647.     ngp[jngp0] = ngp0;
  648.     --jngp1;
  649.     ngp[jngp1] = ngp1;
  650. /* L160: */
  651.     }
  652. /* determines grid points outside the data area. */
  653. /* - in semi-infinite rectangular area. */
  654.     i_1 = nl0;
  655.     for (il0 = 1; il0 <= i_1; ++il0) {
  656.     ngp0 = 0;
  657.     ngp1 = 0;
  658.     il0t3 = il0 * 3;
  659.     ip1 = ipl[il0t3 - 2];
  660.     ip2 = ipl[il0t3 - 1];
  661.     x1 = xd[ip1];
  662.     y1 = yd[ip1];
  663.     x2 = xd[ip2];
  664.     y2 = yd[ip2];
  665.     xmn = ximn;
  666.     xmx = ximx;
  667.     ymn = yimn;
  668.     ymx = yimx;
  669.     if (y2 >= y1) {
  670.         xmn = dmin(x1,x2);
  671.     }
  672.     if (y2 <= y1) {
  673.         xmx = dmax(x1,x2);
  674.     }
  675.     if (x2 <= x1) {
  676.         ymn = dmin(y1,y2);
  677.     }
  678.     if (x2 >= x1) {
  679.         ymx = dmax(y1,y2);
  680.     }
  681.     insd = 0;
  682.     i_2 = nxi0;
  683.     for (ixi = 1; ixi <= i_2; ++ixi) {
  684.         if (xi[ixi] >= xmn && xi[ixi] <= xmx) {
  685.         goto L170;
  686.         }
  687.         if (insd == 0) {
  688.         goto L180;
  689.         }
  690.         iximx = ixi - 1;
  691.         goto L190;
  692. L170:
  693.         if (insd == 1) {
  694.         goto L180;
  695.         }
  696.         insd = 1;
  697.         iximn = ixi;
  698. L180:
  699.     ;}
  700.     if (insd == 0) {
  701.         goto L310;
  702.     }
  703.     iximx = nxi0;
  704. L190:
  705.     i_2 = nyi0;
  706.     for (iyi = 1; iyi <= i_2; ++iyi) {
  707.         yii = yi[iyi];
  708.         if (yii < ymn || yii > ymx) {
  709.         goto L300;
  710.         }
  711.         i_3 = iximx;
  712.         for (ixi = iximn; ixi <= i_3; ++ixi) {
  713.         xii = xi[ixi];
  714.         l = 0;
  715.         if ((r_1 = (x1 - xii) * (y2 - yii) - (y1 - yii) * (x2 - xii)) 
  716.             < (float)0.) {
  717.             goto L210;
  718.         } else if (r_1 == 0) {
  719.             goto L200;
  720.         } else {
  721.             goto L290;
  722.         }
  723. L200:
  724.         l = 1;
  725. L210:
  726.         if ((r_1 = (x2 - x1) * (xii - x1) + (y2 - y1) * (yii - y1)) < 
  727.             (float)0.) {
  728.             goto L290;
  729.         } else if (r_1 == 0) {
  730.             goto L220;
  731.         } else {
  732.             goto L230;
  733.         }
  734. L220:
  735.         l = 1;
  736. L230:
  737.         if ((r_1 = (x1 - x2) * (xii - x2) + (y1 - y2) * (yii - y2)) < 
  738.             (float)0.) {
  739.             goto L290;
  740.         } else if (r_1 == 0) {
  741.             goto L240;
  742.         } else {
  743.             goto L250;
  744.         }
  745. L240:
  746.         l = 1;
  747. L250:
  748.         izi = nxi0 * (iyi - 1) + ixi;
  749.         if (l == 1) {
  750.             goto L260;
  751.         }
  752.         ++ngp0;
  753.         ++jigp0;
  754.         igp[jigp0] = izi;
  755.         goto L290;
  756. L260:
  757.         if (jigp1 > nxinyi) {
  758.             goto L280;
  759.         }
  760.         i_4 = nxinyi;
  761.         for (jigp1i = jigp1; jigp1i <= i_4; ++jigp1i) {
  762.             if (izi == igp[jigp1i]) {
  763.             goto L290;
  764.             }
  765. /* L270: */
  766.         }
  767. L280:
  768.         ++ngp1;
  769.         --jigp1;
  770.         igp[jigp1] = izi;
  771. L290:
  772.         ;}
  773. L300:
  774.     ;}
  775. L310:
  776.     ++jngp0;
  777.     ngp[jngp0] = ngp0;
  778.     --jngp1;
  779.     ngp[jngp1] = ngp1;
  780. /* - in semi-infinite triangular area. */
  781.     ngp0 = 0;
  782.     ngp1 = 0;
  783.     ilp1 = il0 % nl0 + 1;
  784.     ilp1t3 = ilp1 * 3;
  785.     ip3 = ipl[ilp1t3 - 1];
  786.     x3 = xd[ip3];
  787.     y3 = yd[ip3];
  788.     xmn = ximn;
  789.     xmx = ximx;
  790.     ymn = yimn;
  791.     ymx = yimx;
  792.     if (y3 >= y2 && y2 >= y1) {
  793.         xmn = x2;
  794.     }
  795.     if (y3 <= y2 && y2 <= y1) {
  796.         xmx = x2;
  797.     }
  798.     if (x3 <= x2 && x2 <= x1) {
  799.         ymn = y2;
  800.     }
  801.     if (x3 >= x2 && x2 >= x1) {
  802.         ymx = y2;
  803.     }
  804.     insd = 0;
  805.     i_2 = nxi0;
  806.     for (ixi = 1; ixi <= i_2; ++ixi) {
  807.         if (xi[ixi] >= xmn && xi[ixi] <= xmx) {
  808.         goto L320;
  809.         }
  810.         if (insd == 0) {
  811.         goto L330;
  812.         }
  813.         iximx = ixi - 1;
  814.         goto L340;
  815. L320:
  816.         if (insd == 1) {
  817.         goto L330;
  818.         }
  819.         insd = 1;
  820.         iximn = ixi;
  821. L330:
  822.     ;}
  823.     if (insd == 0) {
  824.         goto L440;
  825.     }
  826.     iximx = nxi0;
  827. L340:
  828.     i_2 = nyi0;
  829.     for (iyi = 1; iyi <= i_2; ++iyi) {
  830.         yii = yi[iyi];
  831.         if (yii < ymn || yii > ymx) {
  832.         goto L430;
  833.         }
  834.         i_3 = iximx;
  835.         for (ixi = iximn; ixi <= i_3; ++ixi) {
  836.         xii = xi[ixi];
  837.         l = 0;
  838.         if ((r_1 = (x1 - x2) * (xii - x2) + (y1 - y2) * (yii - y2)) < 
  839.             (float)0.) {
  840.             goto L360;
  841.         } else if (r_1 == 0) {
  842.             goto L350;
  843.         } else {
  844.             goto L420;
  845.         }
  846. L350:
  847.         l = 1;
  848. L360:
  849.         if ((r_1 = (x3 - x2) * (xii - x2) + (y3 - y2) * (yii - y2)) < 
  850.             (float)0.) {
  851.             goto L380;
  852.         } else if (r_1 == 0) {
  853.             goto L370;
  854.         } else {
  855.             goto L420;
  856.         }
  857. L370:
  858.         l = 1;
  859. L380:
  860.         izi = nxi0 * (iyi - 1) + ixi;
  861.         if (l == 1) {
  862.             goto L390;
  863.         }
  864.         ++ngp0;
  865.         ++jigp0;
  866.         igp[jigp0] = izi;
  867.         goto L420;
  868. L390:
  869.         if (jigp1 > nxinyi) {
  870.             goto L410;
  871.         }
  872.         i_4 = nxinyi;
  873.         for (jigp1i = jigp1; jigp1i <= i_4; ++jigp1i) {
  874.             if (izi == igp[jigp1i]) {
  875.             goto L420;
  876.             }
  877. /* L400: */
  878.         }
  879. L410:
  880.         ++ngp1;
  881.         --jigp1;
  882.         igp[jigp1] = izi;
  883. L420:
  884.         ;}
  885. L430:
  886.     ;}
  887. L440:
  888.     ++jngp0;
  889.     ngp[jngp0] = ngp0;
  890.     --jngp1;
  891.     ngp[jngp1] = ngp1;
  892. /* L450: */
  893.     }
  894.     return 0;
  895. } /* idgrid_ */
  896.  
  897. /* Subroutine */ int idlctn_(ndp, xd, yd, nt, ipt, nl, ipl, xii, yii, iti, 
  898.     iwk, wk)
  899. integer *ndp;
  900. real *xd, *yd;
  901. integer *nt, *ipt, *nl, *ipl;
  902. real *xii, *yii;
  903. integer *iti, *iwk;
  904. real *wk;
  905. {
  906.     /* System generated locals */
  907.     integer i_1;
  908.     real r_1, r_2;
  909.  
  910.     /* Local variables */
  911.     static integer idsc[9], il1t3, itsc, it0t3, jiwk, ntsc[9], ntsci, i1, i2, 
  912.         i3, itipv;
  913.     static real x0, y0, x1, y1, x2, y2, x3, y3, xi, yi;
  914.     static integer il1, il2, nl0, ip1, ip2, it0, ip3, nt0;
  915.     static real xs1, xs2, ys1, ys2;
  916.     static integer idp, isc, ntl, jwk;
  917.     static real xmn, ymn, xmx, ymx;
  918.     static integer ndp0;
  919.  
  920. /* this subroutine locates a point, i.e., determines to what tri- */
  921. /* angle a given point (xii,yii) belongs.  when the given point */
  922. /* does not lie inside the data area, this subroutine determines */
  923. /* the border line segment when the point lies in an outside */
  924. /* rectangular area, and two border line segments when the point */
  925. /* lies in an outside triangular area. */
  926. /* the input parameters are */
  927. /*     ndp = number of data points, */
  928. /*     xd,yd = arrays of dimension ndp containing the x and y */
  929. /*           coordinates of the data points, */
  930. /*     nt  = number of triangles, */
  931. /*     ipt = integer array of dimension 3*nt containing the */
  932. /*           point numbers of the vertexes of the triangles, */
  933. /*     nl  = number of border line segments, */
  934. /*     ipl = integer array of dimension 3*nl containing the */
  935. /*           point numbers of the end points of the border */
  936. /*           line segments and their respective triangle */
  937. /*           numbers, */
  938. /*     xii,yii = x and y coordinates of the point to be */
  939. /*           located. */
  940. /* the output parameter is */
  941. /*     iti = triangle number, when the point is inside the */
  942. /*           data area, or */
  943. /*           two border line segment numbers, il1 and il2, */
  944. /*           coded to il1*(nt+nl)+il2, when the point is */
  945. /*           outside the data area. */
  946. /* the other parameters are */
  947. /*     iwk = integer array of dimension 18*ndp used inter- */
  948. /*           nally as a work area, */
  949. /*     wk  = array of dimension 8*ndp used internally as a */
  950. /*           work area. */
  951. /* declaration statements */
  952. /* statement functions */
  953. /* preliminary processing */
  954.     /* Parameter adjustments */
  955.     --wk;
  956.     --iwk;
  957.     --ipl;
  958.     --ipt;
  959.     --yd;
  960.     --xd;
  961.  
  962.     /* Function Body */
  963.     ndp0 = *ndp;
  964.     nt0 = *nt;
  965.     nl0 = *nl;
  966.     ntl = nt0 + nl0;
  967.     x0 = *xii;
  968.     y0 = *yii;
  969. /* processing for a new set of data points */
  970.     if (idlc_1.nit != 0) {
  971.     goto L80;
  972.     }
  973.     idlc_1.nit = 1;
  974. /* - divides the x-y plane into nine rectangular sections. */
  975.     xmn = xd[1];
  976.     xmx = xmn;
  977.     ymn = yd[1];
  978.     ymx = ymn;
  979.     i_1 = ndp0;
  980.     for (idp = 2; idp <= i_1; ++idp) {
  981.     xi = xd[idp];
  982.     yi = yd[idp];
  983.     xmn = dmin(xi,xmn);
  984.     xmx = dmax(xi,xmx);
  985.     ymn = dmin(yi,ymn);
  986.     ymx = dmax(yi,ymx);
  987. /* L10: */
  988.     }
  989.     xs1 = (xmn + xmn + xmx) / (float)3.;
  990.     xs2 = (xmn + xmx + xmx) / (float)3.;
  991.     ys1 = (ymn + ymn + ymx) / (float)3.;
  992.     ys2 = (ymn + ymx + ymx) / (float)3.;
  993. /* - determines and stores in the iwk array triangle numbers of */
  994. /* - the triangles associated with each of the nine sections. */
  995.     for (isc = 1; isc <= 9; ++isc) {
  996.     ntsc[isc - 1] = 0;
  997.     idsc[isc - 1] = 0;
  998. /* L20: */
  999.     }
  1000.     it0t3 = 0;
  1001.     jwk = 0;
  1002.     i_1 = nt0;
  1003.     for (it0 = 1; it0 <= i_1; ++it0) {
  1004.     it0t3 += 3;
  1005.     i1 = ipt[it0t3 - 2];
  1006.     i2 = ipt[it0t3 - 1];
  1007.     i3 = ipt[it0t3];
  1008. /* Computing MIN */
  1009.     r_1 = xd[i1], r_2 = xd[i2], r_1 = min(r_1,r_2), r_2 = xd[i3];
  1010.     xmn = dmin(r_1,r_2);
  1011. /* Computing MAX */
  1012.     r_1 = xd[i1], r_2 = xd[i2], r_1 = max(r_1,r_2), r_2 = xd[i3];
  1013.     xmx = dmax(r_1,r_2);
  1014. /* Computing MIN */
  1015.     r_1 = yd[i1], r_2 = yd[i2], r_1 = min(r_1,r_2), r_2 = yd[i3];
  1016.     ymn = dmin(r_1,r_2);
  1017. /* Computing MAX */
  1018.     r_1 = yd[i1], r_2 = yd[i2], r_1 = max(r_1,r_2), r_2 = yd[i3];
  1019.     ymx = dmax(r_1,r_2);
  1020.     if (ymn > ys1) {
  1021.         goto L30;
  1022.     }
  1023.     if (xmn <= xs1) {
  1024.         idsc[0] = 1;
  1025.     }
  1026.     if (xmx >= xs1 && xmn <= xs2) {
  1027.         idsc[1] = 1;
  1028.     }
  1029.     if (xmx >= xs2) {
  1030.         idsc[2] = 1;
  1031.     }
  1032. L30:
  1033.     if (ymx < ys1 || ymn > ys2) {
  1034.         goto L40;
  1035.     }
  1036.     if (xmn <= xs1) {
  1037.         idsc[3] = 1;
  1038.     }
  1039.     if (xmx >= xs1 && xmn <= xs2) {
  1040.         idsc[4] = 1;
  1041.     }
  1042.     if (xmx >= xs2) {
  1043.         idsc[5] = 1;
  1044.     }
  1045. L40:
  1046.     if (ymx < ys2) {
  1047.         goto L50;
  1048.     }
  1049.     if (xmn <= xs1) {
  1050.         idsc[6] = 1;
  1051.     }
  1052.     if (xmx >= xs1 && xmn <= xs2) {
  1053.         idsc[7] = 1;
  1054.     }
  1055.     if (xmx >= xs2) {
  1056.         idsc[8] = 1;
  1057.     }
  1058. L50:
  1059.     for (isc = 1; isc <= 9; ++isc) {
  1060.         if (idsc[isc - 1] == 0) {
  1061.         goto L60;
  1062.         }
  1063.         jiwk = ntsc[isc - 1] * 9 + isc;
  1064.         iwk[jiwk] = it0;
  1065.         ++ntsc[isc - 1];
  1066.         idsc[isc - 1] = 0;
  1067. L60:
  1068.     ;}
  1069. /* - stores in the wk array the minimum and maximum of the x and */
  1070. /* - y coordinate values for each of the triangle. */
  1071.     jwk += 4;
  1072.     wk[jwk - 3] = xmn;
  1073.     wk[jwk - 2] = xmx;
  1074.     wk[jwk - 1] = ymn;
  1075.     wk[jwk] = ymx;
  1076. /* L70: */
  1077.     }
  1078.     goto L110;
  1079. /* checks if in the same triangle as previous. */
  1080. L80:
  1081.     it0 = itipv;
  1082.     if (it0 > nt0) {
  1083.     goto L90;
  1084.     }
  1085.     it0t3 = it0 * 3;
  1086.     ip1 = ipt[it0t3 - 2];
  1087.     x1 = xd[ip1];
  1088.     y1 = yd[ip1];
  1089.     ip2 = ipt[it0t3 - 1];
  1090.     x2 = xd[ip2];
  1091.     y2 = yd[ip2];
  1092.     if ((x1 - x0) * (y2 - y0) - (y1 - y0) * (x2 - x0) < (float)0.) {
  1093.     goto L110;
  1094.     }
  1095.     ip3 = ipt[it0t3];
  1096.     x3 = xd[ip3];
  1097.     y3 = yd[ip3];
  1098.     if ((x2 - x0) * (y3 - y0) - (y2 - y0) * (x3 - x0) < (float)0.) {
  1099.     goto L110;
  1100.     }
  1101.     if ((x3 - x0) * (y1 - y0) - (y3 - y0) * (x1 - x0) < (float)0.) {
  1102.     goto L110;
  1103.     }
  1104.     goto L170;
  1105. /* checks if on the same border line segment. */
  1106. L90:
  1107.     il1 = it0 / ntl;
  1108.     il2 = it0 - il1 * ntl;
  1109.     il1t3 = il1 * 3;
  1110.     ip1 = ipl[il1t3 - 2];
  1111.     x1 = xd[ip1];
  1112.     y1 = yd[ip1];
  1113.     ip2 = ipl[il1t3 - 1];
  1114.     x2 = xd[ip2];
  1115.     y2 = yd[ip2];
  1116.     if (il2 != il1) {
  1117.     goto L100;
  1118.     }
  1119.     if ((x1 - x2) * (x0 - x2) + (y1 - y2) * (y0 - y2) < (float)0.) {
  1120.     goto L110;
  1121.     }
  1122.     if ((x2 - x1) * (x0 - x1) + (y2 - y1) * (y0 - y1) < (float)0.) {
  1123.     goto L110;
  1124.     }
  1125.     if ((x1 - x0) * (y2 - y0) - (y1 - y0) * (x2 - x0) > (float)0.) {
  1126.     goto L110;
  1127.     }
  1128.     goto L170;
  1129. /* checks if between the same two border line segments. */
  1130. L100:
  1131.     if ((x1 - x2) * (x0 - x2) + (y1 - y2) * (y0 - y2) > (float)0.) {
  1132.     goto L110;
  1133.     }
  1134.     ip3 = ipl[il2 * 3 - 1];
  1135.     x3 = xd[ip3];
  1136.     y3 = yd[ip3];
  1137.     if ((x3 - x2) * (x0 - x2) + (y3 - y2) * (y0 - y2) <= (float)0.) {
  1138.     goto L170;
  1139.     }
  1140. /* locates inside the data area. */
  1141. /* - determines the section in which the point in question lies. */
  1142. L110:
  1143.     isc = 1;
  1144.     if (x0 >= xs1) {
  1145.     ++isc;
  1146.     }
  1147.     if (x0 >= xs2) {
  1148.     ++isc;
  1149.     }
  1150.     if (y0 >= ys1) {
  1151.     isc += 3;
  1152.     }
  1153.     if (y0 >= ys2) {
  1154.     isc += 3;
  1155.     }
  1156. /* - searches through the triangles associated with the section. */
  1157.     ntsci = ntsc[isc - 1];
  1158.     if (ntsci <= 0) {
  1159.     goto L130;
  1160.     }
  1161.     jiwk = isc - 9;
  1162.     i_1 = ntsci;
  1163.     for (itsc = 1; itsc <= i_1; ++itsc) {
  1164.     jiwk += 9;
  1165.     it0 = iwk[jiwk];
  1166.     jwk = it0 << 2;
  1167.     if (x0 < wk[jwk - 3]) {
  1168.         goto L120;
  1169.     }
  1170.     if (x0 > wk[jwk - 2]) {
  1171.         goto L120;
  1172.     }
  1173.     if (y0 < wk[jwk - 1]) {
  1174.         goto L120;
  1175.     }
  1176.     if (y0 > wk[jwk]) {
  1177.         goto L120;
  1178.     }
  1179.     it0t3 = it0 * 3;
  1180.     ip1 = ipt[it0t3 - 2];
  1181.     x1 = xd[ip1];
  1182.     y1 = yd[ip1];
  1183.     ip2 = ipt[it0t3 - 1];
  1184.     x2 = xd[ip2];
  1185.     y2 = yd[ip2];
  1186.     if ((x1 - x0) * (y2 - y0) - (y1 - y0) * (x2 - x0) < (float)0.) {
  1187.         goto L120;
  1188.     }
  1189.     ip3 = ipt[it0t3];
  1190.     x3 = xd[ip3];
  1191.     y3 = yd[ip3];
  1192.     if ((x2 - x0) * (y3 - y0) - (y2 - y0) * (x3 - x0) < (float)0.) {
  1193.         goto L120;
  1194.     }
  1195.     if ((x3 - x0) * (y1 - y0) - (y3 - y0) * (x1 - x0) < (float)0.) {
  1196.         goto L120;
  1197.     }
  1198.     goto L170;
  1199. L120:
  1200.     ;}
  1201. /* locates outside the data area. */
  1202. L130:
  1203.     i_1 = nl0;
  1204.     for (il1 = 1; il1 <= i_1; ++il1) {
  1205.     il1t3 = il1 * 3;
  1206.     ip1 = ipl[il1t3 - 2];
  1207.     x1 = xd[ip1];
  1208.     y1 = yd[ip1];
  1209.     ip2 = ipl[il1t3 - 1];
  1210.     x2 = xd[ip2];
  1211.     y2 = yd[ip2];
  1212.     if ((x2 - x1) * (x0 - x1) + (y2 - y1) * (y0 - y1) < (float)0.) {
  1213.         goto L150;
  1214.     }
  1215.     if ((x1 - x2) * (x0 - x2) + (y1 - y2) * (y0 - y2) < (float)0.) {
  1216.         goto L140;
  1217.     }
  1218.     if ((x1 - x0) * (y2 - y0) - (y1 - y0) * (x2 - x0) > (float)0.) {
  1219.         goto L150;
  1220.     }
  1221.     il2 = il1;
  1222.     goto L160;
  1223. L140:
  1224.     il2 = il1 % nl0 + 1;
  1225.     ip3 = ipl[il2 * 3 - 1];
  1226.     x3 = xd[ip3];
  1227.     y3 = yd[ip3];
  1228.     if ((x3 - x2) * (x0 - x2) + (y3 - y2) * (y0 - y2) <= (float)0.) {
  1229.         goto L160;
  1230.     }
  1231. L150:
  1232.     ;}
  1233.     it0 = 1;
  1234.     goto L170;
  1235. L160:
  1236.     it0 = il1 * ntl + il2;
  1237. /* normal exit */
  1238. L170:
  1239.     *iti = it0;
  1240.     itipv = it0;
  1241.     return 0;
  1242. } /* idlctn_ */
  1243.  
  1244. /* Subroutine */ int idpdrv_(ndp, xd, yd, zd, ncp, ipc, pd)
  1245. integer *ndp;
  1246. real *xd, *yd, *zd;
  1247. integer *ncp, *ipc;
  1248. real *pd;
  1249. {
  1250.     /* System generated locals */
  1251.     integer i_1, i_2, i_3;
  1252.  
  1253.     /* Local variables */
  1254.     static integer jipc;
  1255.     static real dnmx, dnmy, dnmz, nmxx, nmxy, nmyx, nmyy;
  1256.     static integer jipc0, ic2mn, ncpm1;
  1257.     static real dnmxx, dnmxy, dnmyx, dnmyy, x0, y0, z0;
  1258.     static integer ic1, ic2, ip0;
  1259.     static real dx1, dy1, dz1, dx2, dy2, dz2, zx0, zy0;
  1260.     static integer jpd, ipi;
  1261.     static real nmx, nmy, nmz;
  1262.     static integer jpd0, ncp0, ndp0;
  1263.     static real dzx1, dzy1, dzx2, dzy2;
  1264.  
  1265. /* this subroutine estimates partial derivatives of the first and */
  1266. /* second order at the data points. */
  1267. /* the input parameters are */
  1268. /*     ndp = number of data points, */
  1269. /*     xd,yd,zd = arrays of dimension ndp containing the x, */
  1270. /*           y, and z coordinates of the data points, */
  1271. /*     ncp = number of additional data points used for esti- */
  1272. /*           mating partial derivatives at each data point, */
  1273. /*     ipc = integer array of dimension ncp*ndp containing */
  1274. /*           the point numbers of ncp data points closest to */
  1275. /*           each of the ndp data points. */
  1276. /* the output parameter is */
  1277. /*     pd  = array of dimension 5*ndp, where the estimated */
  1278. /*           zx, zy, zxx, zxy, and zyy values at the data */
  1279. /*           points are to be stored. */
  1280. /* declaration statements */
  1281. /* preliminary processing */
  1282.     /* Parameter adjustments */
  1283.     --pd;
  1284.     --ipc;
  1285.     --zd;
  1286.     --yd;
  1287.     --xd;
  1288.  
  1289.     /* Function Body */
  1290. /* L10: */
  1291.     ndp0 = *ndp;
  1292.     ncp0 = *ncp;
  1293.     ncpm1 = ncp0 - 1;
  1294. /* estimation of zx and zy */
  1295. /* L20: */
  1296.     i_1 = ndp0;
  1297.     for (ip0 = 1; ip0 <= i_1; ++ip0) {
  1298.     x0 = xd[ip0];
  1299.     y0 = yd[ip0];
  1300.     z0 = zd[ip0];
  1301.     nmx = (float)0.;
  1302.     nmy = (float)0.;
  1303.     nmz = (float)0.;
  1304.     jipc0 = ncp0 * (ip0 - 1);
  1305.     i_2 = ncpm1;
  1306.     for (ic1 = 1; ic1 <= i_2; ++ic1) {
  1307.         jipc = jipc0 + ic1;
  1308.         ipi = ipc[jipc];
  1309.         dx1 = xd[ipi] - x0;
  1310.         dy1 = yd[ipi] - y0;
  1311.         dz1 = zd[ipi] - z0;
  1312.         ic2mn = ic1 + 1;
  1313.         i_3 = ncp0;
  1314.         for (ic2 = ic2mn; ic2 <= i_3; ++ic2) {
  1315.         jipc = jipc0 + ic2;
  1316.         ipi = ipc[jipc];
  1317.         dx2 = xd[ipi] - x0;
  1318.         dy2 = yd[ipi] - y0;
  1319.         dnmz = dx1 * dy2 - dy1 * dx2;
  1320.         if (dnmz == (float)0.) {
  1321.             goto L22;
  1322.         }
  1323.         dz2 = zd[ipi] - z0;
  1324.         dnmx = dy1 * dz2 - dz1 * dy2;
  1325.         dnmy = dz1 * dx2 - dx1 * dz2;
  1326.         if (dnmz >= (float)0.) {
  1327.             goto L21;
  1328.         }
  1329.         dnmx = -(doublereal)dnmx;
  1330.         dnmy = -(doublereal)dnmy;
  1331.         dnmz = -(doublereal)dnmz;
  1332. L21:
  1333.         nmx += dnmx;
  1334.         nmy += dnmy;
  1335.         nmz += dnmz;
  1336. L22:
  1337.         ;}
  1338. /* L23: */
  1339.     }
  1340.     jpd0 = ip0 * 5;
  1341.     pd[jpd0 - 4] = -(doublereal)nmx / nmz;
  1342.     pd[jpd0 - 3] = -(doublereal)nmy / nmz;
  1343. /* L24: */
  1344.     }
  1345. /* estimation of zxx, zxy, and zyy */
  1346. /* L30: */
  1347.     i_1 = ndp0;
  1348.     for (ip0 = 1; ip0 <= i_1; ++ip0) {
  1349.     jpd0 += 5;
  1350.     x0 = xd[ip0];
  1351.     jpd0 = ip0 * 5;
  1352.     y0 = yd[ip0];
  1353.     zx0 = pd[jpd0 - 4];
  1354.     zy0 = pd[jpd0 - 3];
  1355.     nmxx = (float)0.;
  1356.     nmxy = (float)0.;
  1357.     nmyx = (float)0.;
  1358.     nmyy = (float)0.;
  1359.     nmz = (float)0.;
  1360.     jipc0 = ncp0 * (ip0 - 1);
  1361.     i_2 = ncpm1;
  1362.     for (ic1 = 1; ic1 <= i_2; ++ic1) {
  1363.         jipc = jipc0 + ic1;
  1364.         ipi = ipc[jipc];
  1365.         dx1 = xd[ipi] - x0;
  1366.         dy1 = yd[ipi] - y0;
  1367.         jpd = ipi * 5;
  1368.         dzx1 = pd[jpd - 4] - zx0;
  1369.         dzy1 = pd[jpd - 3] - zy0;
  1370.         ic2mn = ic1 + 1;
  1371.         i_3 = ncp0;
  1372.         for (ic2 = ic2mn; ic2 <= i_3; ++ic2) {
  1373.         jipc = jipc0 + ic2;
  1374.         ipi = ipc[jipc];
  1375.         dx2 = xd[ipi] - x0;
  1376.         dy2 = yd[ipi] - y0;
  1377.         dnmz = dx1 * dy2 - dy1 * dx2;
  1378.         if (dnmz == (float)0.) {
  1379.             goto L32;
  1380.         }
  1381.         jpd = ipi * 5;
  1382.         dzx2 = pd[jpd - 4] - zx0;
  1383.         dzy2 = pd[jpd - 3] - zy0;
  1384.         dnmxx = dy1 * dzx2 - dzx1 * dy2;
  1385.         dnmxy = dzx1 * dx2 - dx1 * dzx2;
  1386.         dnmyx = dy1 * dzy2 - dzy1 * dy2;
  1387.         dnmyy = dzy1 * dx2 - dx1 * dzy2;
  1388.         if (dnmz >= (float)0.) {
  1389.             goto L31;
  1390.         }
  1391.         dnmxx = -(doublereal)dnmxx;
  1392.         dnmxy = -(doublereal)dnmxy;
  1393.         dnmyx = -(doublereal)dnmyx;
  1394.         dnmyy = -(doublereal)dnmyy;
  1395.         dnmz = -(doublereal)dnmz;
  1396. L31:
  1397.         nmxx += dnmxx;
  1398.         nmxy += dnmxy;
  1399.         nmyx += dnmyx;
  1400.         nmyy += dnmyy;
  1401.         nmz += dnmz;
  1402. L32:
  1403.         ;}
  1404. /* L33: */
  1405.     }
  1406.     pd[jpd0 - 2] = -(doublereal)nmxx / nmz;
  1407.     pd[jpd0 - 1] = -(doublereal)(nmxy + nmyx) / (nmz * (float)2.);
  1408.     pd[jpd0] = -(doublereal)nmyy / nmz;
  1409. /* L34: */
  1410.     }
  1411.     return 0;
  1412. } /* idpdrv_ */
  1413.  
  1414. /* Subroutine */ int idptip_(xd, yd, zd, nt, ipt, nl, ipl, pdd, iti, xii, yii,
  1415.      zii)
  1416. real *xd, *yd, *zd;
  1417. integer *nt, *ipt, *nl, *ipl;
  1418. real *pdd;
  1419. integer *iti;
  1420. real *xii, *yii, *zii;
  1421. {
  1422.     /* System generated locals */
  1423.     static real equiv_0[1];
  1424.  
  1425.     /* Builtin functions */
  1426.     double sqrt(), atan2(), cos(), sin();
  1427.  
  1428.     /* Local variables */
  1429.     static integer jpdd, jipl, jipt;
  1430.     static real csuv, thus, thsv, thuv, thxu, a, b, c, d;
  1431.     static integer i;
  1432.     static real u, v, x[3], y[3], z[3], g1, h1, h2, h3, g2, p0, p1, p2, p3, 
  1433.         p4;
  1434. #define p5 (equiv_0)
  1435.     static real x0, y0, aa, ab, bb, ad, bc, cc, cd, dd, ac, p00, ap, bp, cp, 
  1436.         pd[15];
  1437. #define p50 (equiv_0)
  1438.     static real dp, p10, p01, p20, p11, p02, p30, lu, lv, p40, p03, p04, p05, 
  1439.         p41, p14, p21, p31, p12, p13, p22, zu[3], zv[3], p32, p23, dx, dy;
  1440.  
  1441.     static integer il1, il2, it0, idp, jpd, kpd;
  1442.     static real dlt;
  1443.     static integer ntl;
  1444.     static real zuu[3], zuv[3], zvv[3], act2, bdt2, adbc;
  1445.  
  1446. /* this subroutine performs punctual interpolation or extrapola- */
  1447. /* tion, i.e., determines the z value at a point. */
  1448. /* the input parameters are */
  1449. /*     xd,yd,zd = arrays of dimension ndp containing the x, */
  1450. /*           y, and z coordinates of the data points, where */
  1451. /*           ndp is the number of the data points, */
  1452. /*     nt  = number of triangles, */
  1453. /*     ipt = integer array of dimension 3*nt containing the */
  1454. /*           point numbers of the vertexes of the triangles, */
  1455. /*     nl  = number of border line segments, */
  1456. /*     ipl = integer array of dimension 3*nl containing the */
  1457. /*           point numbers of the end points of the border */
  1458. /*           line segments and their respective triangle */
  1459. /*           numbers, */
  1460. /*     pdd = array of dimension 5*ndp containing the partial */
  1461. /*           derivatives at the data points, */
  1462. /*     iti = triangle number of the triangle in which lies */
  1463. /*           the point for which interpolation is to be */
  1464. /*           performed, */
  1465. /*     xii,yii = x and y coordinates of the point for which */
  1466. /*           interpolation is to be performed. */
  1467. /* the output parameter is */
  1468. /*     zii = interpolated z value. */
  1469. /* declaration statements */
  1470. /* preliminary processing */
  1471.     /* Parameter adjustments */
  1472.     --pdd;
  1473.     --ipl;
  1474.     --ipt;
  1475.     --zd;
  1476.     --yd;
  1477.     --xd;
  1478.  
  1479.     /* Function Body */
  1480. /* L10: */
  1481.     it0 = *iti;
  1482.     ntl = *nt + *nl;
  1483.     if (it0 <= ntl) {
  1484.     goto L20;
  1485.     }
  1486.     il1 = it0 / ntl;
  1487.     il2 = it0 - il1 * ntl;
  1488.     if (il1 == il2) {
  1489.     goto L40;
  1490.     }
  1491.     goto L60;
  1492. /* calculation of zii by interpolation. */
  1493. /* checks if the necessary coefficients have been calculated. */
  1494. L20:
  1495.     if (it0 == idpi_1.itpv) {
  1496.     goto L30;
  1497.     }
  1498. /* loads coordinate and partial derivative values at the */
  1499. /* vertexes. */
  1500. /* L21: */
  1501.     jipt = (it0 - 1) * 3;
  1502.     jpd = 0;
  1503.     for (i = 1; i <= 3; ++i) {
  1504.     ++jipt;
  1505.     idp = ipt[jipt];
  1506.     x[i - 1] = xd[idp];
  1507.     y[i - 1] = yd[idp];
  1508.     z[i - 1] = zd[idp];
  1509.     jpdd = (idp - 1) * 5;
  1510.     for (kpd = 1; kpd <= 5; ++kpd) {
  1511.         ++jpd;
  1512.         ++jpdd;
  1513.         pd[jpd - 1] = pdd[jpdd];
  1514. /* L22: */
  1515.     }
  1516. /* L23: */
  1517.     }
  1518. /* determines the coefficients for the coordinate system */
  1519. /* transformation from the x-y system to the u-v system */
  1520. /* and vice versa. */
  1521. /* L24: */
  1522.     x0 = x[0];
  1523.     y0 = y[0];
  1524.     a = x[1] - x0;
  1525.     b = x[2] - x0;
  1526.     c = y[1] - y0;
  1527.     d = y[2] - y0;
  1528.     ad = a * d;
  1529.     bc = b * c;
  1530.     dlt = ad - bc;
  1531.     ap = d / dlt;
  1532.     bp = -(doublereal)b / dlt;
  1533.     cp = -(doublereal)c / dlt;
  1534.     dp = a / dlt;
  1535. /* converts the partial derivatives at the vertexes of the */
  1536. /* triangle for the u-v coordinate system. */
  1537. /* L25: */
  1538.     aa = a * a;
  1539.     act2 = a * (float)2. * c;
  1540.     cc = c * c;
  1541.     ab = a * b;
  1542.     adbc = ad + bc;
  1543.     cd = c * d;
  1544.     bb = b * b;
  1545.     bdt2 = b * (float)2. * d;
  1546.     dd = d * d;
  1547.     for (i = 1; i <= 3; ++i) {
  1548.     jpd = i * 5;
  1549.     zu[i - 1] = a * pd[jpd - 5] + c * pd[jpd - 4];
  1550.     zv[i - 1] = b * pd[jpd - 5] + d * pd[jpd - 4];
  1551.     zuu[i - 1] = aa * pd[jpd - 3] + act2 * pd[jpd - 2] + cc * pd[jpd - 1];
  1552.  
  1553.     zuv[i - 1] = ab * pd[jpd - 3] + adbc * pd[jpd - 2] + cd * pd[jpd - 1];
  1554.  
  1555.     zvv[i - 1] = bb * pd[jpd - 3] + bdt2 * pd[jpd - 2] + dd * pd[jpd - 1];
  1556.  
  1557. /* L26: */
  1558.     }
  1559. /* calculates the coefficients of the polynomial. */
  1560. /* L27: */
  1561.     p00 = z[0];
  1562.     p10 = zu[0];
  1563.     p01 = zv[0];
  1564.     p20 = zuu[0] * (float).5;
  1565.     p11 = zuv[0];
  1566.     p02 = zvv[0] * (float).5;
  1567.     h1 = z[1] - p00 - p10 - p20;
  1568.     h2 = zu[1] - p10 - zuu[0];
  1569.     h3 = zuu[1] - zuu[0];
  1570.     p30 = h1 * (float)10. - h2 * (float)4. + h3 * (float).5;
  1571.     p40 = h1 * (float)-15. + h2 * (float)7. - h3;
  1572.     *p50 = h1 * (float)6. - h2 * (float)3. + h3 * (float).5;
  1573.     h1 = z[2] - p00 - p01 - p02;
  1574.     h2 = zv[2] - p01 - zvv[0];
  1575.     h3 = zvv[2] - zvv[0];
  1576.     p03 = h1 * (float)10. - h2 * (float)4. + h3 * (float).5;
  1577.     p04 = h1 * (float)-15. + h2 * (float)7. - h3;
  1578.     p05 = h1 * (float)6. - h2 * (float)3. + h3 * (float).5;
  1579.     lu = sqrt(aa + cc);
  1580.     lv = sqrt(bb + dd);
  1581.     thxu = atan2(c, a);
  1582.     thuv = atan2(d, b) - thxu;
  1583.     csuv = cos(thuv);
  1584.     p41 = lv * (float)5. * csuv / lu * *p50;
  1585.     p14 = lu * (float)5. * csuv / lv * p05;
  1586.     h1 = zv[1] - p01 - p11 - p41;
  1587.     h2 = zuv[1] - p11 - p41 * (float)4.;
  1588.     p21 = h1 * (float)3. - h2;
  1589.     p31 = h1 * (float)-2. + h2;
  1590.     h1 = zu[2] - p10 - p11 - p14;
  1591.     h2 = zuv[2] - p11 - p14 * (float)4.;
  1592.     p12 = h1 * (float)3. - h2;
  1593.     p13 = h1 * (float)-2. + h2;
  1594.     thus = atan2(d - c, b - a) - thxu;
  1595.     thsv = thuv - thus;
  1596.     aa = sin(thsv) / lu;
  1597.     bb = -(doublereal)cos(thsv) / lu;
  1598.     cc = sin(thus) / lv;
  1599.     dd = cos(thus) / lv;
  1600.     ac = aa * cc;
  1601.     ad = aa * dd;
  1602.     bc = bb * cc;
  1603.     g1 = aa * ac * (bc * (float)3. + ad * (float)2.);
  1604.     g2 = cc * ac * (ad * (float)3. + bc * (float)2.);
  1605.     h1 = -(doublereal)aa * aa * aa * (aa * (float)5. * bb * *p50 + (bc * (
  1606.         float)4. + ad) * p41) - cc * cc * cc * (cc * (float)5. * dd * p05 
  1607.         + (ad * (float)4. + bc) * p14);
  1608.     h2 = zvv[1] * (float).5 - p02 - p12;
  1609.     h3 = zuu[2] * (float).5 - p20 - p21;
  1610.     p22 = (g1 * h2 + g2 * h3 - h1) / (g1 + g2);
  1611.     p32 = h2 - p22;
  1612.     p23 = h3 - p22;
  1613.     idpi_1.itpv = it0;
  1614. /* converts xii and yii to u-v system. */
  1615. L30:
  1616.     dx = *xii - x0;
  1617.     dy = *yii - y0;
  1618.     u = ap * dx + bp * dy;
  1619.     v = cp * dx + dp * dy;
  1620. /* evaluates the polynomial. */
  1621. /* L31: */
  1622.     p0 = p00 + v * (p01 + v * (p02 + v * (p03 + v * (p04 + v * p05))));
  1623.     p1 = p10 + v * (p11 + v * (p12 + v * (p13 + v * p14)));
  1624.     p2 = p20 + v * (p21 + v * (p22 + v * p23));
  1625.     p3 = p30 + v * (p31 + v * p32);
  1626.     p4 = p40 + v * p41;
  1627.     *zii = p0 + u * (p1 + u * (p2 + u * (p3 + u * (p4 + u * *p5))));
  1628.     return 0;
  1629. /* calculation of zii by extrapolation in the rectangle. */
  1630. /* checks if the necessary coefficients have been calculated. */
  1631. L40:
  1632.     if (it0 == idpi_1.itpv) {
  1633.     goto L50;
  1634.     }
  1635. /* loads coordinate and partial derivative values at the end */
  1636. /* points of the border line segment. */
  1637. /* L41: */
  1638.     jipl = (il1 - 1) * 3;
  1639.     jpd = 0;
  1640.     for (i = 1; i <= 2; ++i) {
  1641.     ++jipl;
  1642.     idp = ipl[jipl];
  1643.     x[i - 1] = xd[idp];
  1644.     y[i - 1] = yd[idp];
  1645.     z[i - 1] = zd[idp];
  1646.     jpdd = (idp - 1) * 5;
  1647.     for (kpd = 1; kpd <= 5; ++kpd) {
  1648.         ++jpd;
  1649.         ++jpdd;
  1650.         pd[jpd - 1] = pdd[jpdd];
  1651. /* L42: */
  1652.     }
  1653. /* L43: */
  1654.     }
  1655. /* determines the coefficients for the coordinate system */
  1656. /* transformation from the x-y system to the u-v system */
  1657. /* and vice versa. */
  1658. /* L44: */
  1659.     x0 = x[0];
  1660.     y0 = y[0];
  1661.     a = y[1] - y[0];
  1662.     b = x[1] - x[0];
  1663.     c = -(doublereal)b;
  1664.     d = a;
  1665.     ad = a * d;
  1666.     bc = b * c;
  1667.     dlt = ad - bc;
  1668.     ap = d / dlt;
  1669.     bp = -(doublereal)b / dlt;
  1670.     cp = -(doublereal)bp;
  1671.     dp = ap;
  1672. /* converts the partial derivatives at the end points of the */
  1673. /* border line segment for the u-v coordinate system. */
  1674. /* L45: */
  1675.     aa = a * a;
  1676.     act2 = a * (float)2. * c;
  1677.     cc = c * c;
  1678.     ab = a * b;
  1679.     adbc = ad + bc;
  1680.     cd = c * d;
  1681.     bb = b * b;
  1682.     bdt2 = b * (float)2. * d;
  1683.     dd = d * d;
  1684.     for (i = 1; i <= 2; ++i) {
  1685.     jpd = i * 5;
  1686.     zu[i - 1] = a * pd[jpd - 5] + c * pd[jpd - 4];
  1687.     zv[i - 1] = b * pd[jpd - 5] + d * pd[jpd - 4];
  1688.     zuu[i - 1] = aa * pd[jpd - 3] + act2 * pd[jpd - 2] + cc * pd[jpd - 1];
  1689.  
  1690.     zuv[i - 1] = ab * pd[jpd - 3] + adbc * pd[jpd - 2] + cd * pd[jpd - 1];
  1691.  
  1692.     zvv[i - 1] = bb * pd[jpd - 3] + bdt2 * pd[jpd - 2] + dd * pd[jpd - 1];
  1693.  
  1694. /* L46: */
  1695.     }
  1696. /* calculates the coefficients of the polynomial. */
  1697. /* L47: */
  1698.     p00 = z[0];
  1699.     p10 = zu[0];
  1700.     p01 = zv[0];
  1701.     p20 = zuu[0] * (float).5;
  1702.     p11 = zuv[0];
  1703.     p02 = zvv[0] * (float).5;
  1704.     h1 = z[1] - p00 - p01 - p02;
  1705.     h2 = zv[1] - p01 - zvv[0];
  1706.     h3 = zvv[1] - zvv[0];
  1707.     p03 = h1 * (float)10. - h2 * (float)4. + h3 * (float).5;
  1708.     p04 = h1 * (float)-15. + h2 * (float)7. - h3;
  1709.     p05 = h1 * (float)6. - h2 * (float)3. + h3 * (float).5;
  1710.     h1 = zu[1] - p10 - p11;
  1711.     h2 = zuv[1] - p11;
  1712.     p12 = h1 * (float)3. - h2;
  1713.     p13 = h1 * (float)-2. + h2;
  1714.     p21 = (float)0.;
  1715.     p23 = -(doublereal)zuu[1] + zuu[0];
  1716.     p22 = p23 * (float)-1.5;
  1717.     idpi_1.itpv = it0;
  1718. /* converts xii and yii to u-v system. */
  1719. L50:
  1720.     dx = *xii - x0;
  1721.     dy = *yii - y0;
  1722.     u = ap * dx + bp * dy;
  1723.     v = cp * dx + dp * dy;
  1724. /* evaluates the polynomial. */
  1725. /* L51: */
  1726.     p0 = p00 + v * (p01 + v * (p02 + v * (p03 + v * (p04 + v * p05))));
  1727.     p1 = p10 + v * (p11 + v * (p12 + v * p13));
  1728.     p2 = p20 + v * (p21 + v * (p22 + v * p23));
  1729.     *zii = p0 + u * (p1 + u * p2);
  1730.     return 0;
  1731. /* calculation of zii by extrapolation in the triangle. */
  1732. /* checks if the necessary coefficients have been calculated. */
  1733. L60:
  1734.     if (it0 == idpi_1.itpv) {
  1735.     goto L70;
  1736.     }
  1737. /* loads coordinate and partial derivative values at the vertex */
  1738. /* of the triangle. */
  1739. /* L61: */
  1740.     jipl = il2 * 3 - 2;
  1741.     idp = ipl[jipl];
  1742.     x[0] = xd[idp];
  1743.     y[0] = yd[idp];
  1744.     z[0] = zd[idp];
  1745.     jpdd = (idp - 1) * 5;
  1746.     for (kpd = 1; kpd <= 5; ++kpd) {
  1747.     ++jpdd;
  1748.     pd[kpd - 1] = pdd[jpdd];
  1749. /* L62: */
  1750.     }
  1751. /* calculates the coefficients of the polynomial. */
  1752. /* L67: */
  1753.     p00 = z[0];
  1754.     p10 = pd[0];
  1755.     p01 = pd[1];
  1756.     p20 = pd[2] * (float).5;
  1757.     p11 = pd[3];
  1758.     p02 = pd[4] * (float).5;
  1759.     idpi_1.itpv = it0;
  1760. /* converts xii and yii to u-v system. */
  1761. L70:
  1762.     u = *xii - x[0];
  1763.     v = *yii - y[0];
  1764. /* evaluates the polynomial. */
  1765. /* L71: */
  1766.     p0 = p00 + v * (p01 + v * p02);
  1767.     p1 = p10 + v * p11;
  1768.     *zii = p0 + u * (p1 + u * p20);
  1769.     return 0;
  1770. } /* idptip_ */
  1771.  
  1772. #undef p50
  1773. #undef p5
  1774.  
  1775.  
  1776. /* Subroutine */ int idsfft_(md, ncp, ndp, xd, yd, zd, nxi, nyi, xi, yi, zi, 
  1777.     iwk, wk)
  1778. integer *md, *ncp, *ndp;
  1779. real *xd, *yd, *zd;
  1780. integer *nxi, *nyi;
  1781. real *xi, *yi, *zi;
  1782. integer *iwk;
  1783. real *wk;
  1784. {
  1785.  
  1786.     /* System generated locals */
  1787.     integer i_1, i_2;
  1788.  
  1789.     /* Local variables */
  1790.     static integer jigp, jngp, nngp;
  1791.     extern /* Subroutine */ int err2090_();
  1792.     static integer jwipc, jwigp, jwipl, ncppv, ndppv, jwngp, jwiwl, jwipt, 
  1793.         jwiwp, nxipv, nyipv, jig0mn, jig1mn, jig0mx, jig1mx, jwigp0, 
  1794.         jwngp0, nl;
  1795.     extern /* Subroutine */ int idcldp_();
  1796.     static integer nt;
  1797.     extern /* Subroutine */ int idtang_(), idgrid_(), idpdrv_(), idptip_();
  1798.     static integer md0, il1, il2, iti, ixi, izi, iyi, ncp0, ndp0, ngp0, ngp1, 
  1799.         nxi0, nyi0;
  1800.  
  1801. /* this subroutine performs smooth surface fitting when the pro- */
  1802. /* jections of the data points in the x-y plane are irregularly */
  1803. /* distributed in the plane. */
  1804. /* the input parameters are */
  1805. /*     md  = mode of computation (must be 1, 2, or 3), */
  1806. /*         = 1 for new ncp and/or new xd-yd, */
  1807. /*         = 2 for old ncp, old xd-yd, new xi-yi, */
  1808. /*         = 3 for old ncp, old xd-yd, old xi-yi, */
  1809. /*     ncp = number of additional data points used for esti- */
  1810. /*           mating partial derivatives at each data point */
  1811. /*           (must be 2 or greater, but smaller than ndp), */
  1812. /*     ndp = number of data points (must be 4 or greater), */
  1813. /*     xd  = array of dimension ndp containing the x */
  1814. /*           coordinates of the data points, */
  1815. /*     yd  = array of dimension ndp containing the y */
  1816. /*           coordinates of the data points, */
  1817. /*     zd  = array of dimension ndp containing the z */
  1818. /*           coordinates of the data points, */
  1819. /*     nxi = number of output grid points in the x coordinate */
  1820. /*           (must be 1 or greater), */
  1821. /*     nyi = number of output grid points in the y coordinate */
  1822. /*           (must be 1 or greater), */
  1823. /*     xi  = array of dimension nxi containing the x */
  1824. /*           coordinates of the output grid points, */
  1825. /*     yi  = array of dimension nyi containing the y */
  1826. /*           coordinates of the output grid points. */
  1827. /* the output parameter is */
  1828. /*     zi  = doubly-dimensioned array of dimension (nxi,nyi), */
  1829. /*           where the interpolated z values at the output */
  1830. /*           grid points are to be stored. */
  1831. /* the other parameters are */
  1832. /*     iwk = integer array of dimension */
  1833. /*              max0(31,27+ncp)*ndp+nxi*nyi */
  1834. /*           used internally as a work area, */
  1835. /*     wk  = array of dimension 5*ndp used internally as a */
  1836. /*           work area. */
  1837. /* the very first call to this subroutine and the call with a new */
  1838. /* ncp value, a new ndp value, and/or new contents of the xd and */
  1839. /* yd arrays must be made with md=1.  the call with md=2 must be */
  1840. /* preceded by another call with the same ncp and ndp values and */
  1841. /* with the same contents of the xd and yd arrays.  the call with */
  1842. /* md=3 must be preceded by another call with the same ncp, ndp, */
  1843. /* nxi, and nyi values and with the same contents of the xd, yd, */
  1844. /* xi, and yi arrays.  between the call with md=2 or md=3 and its */
  1845. /* preceding call, the iwk and wk arrays must not be disturbed. */
  1846. /* use of a value between 3 and 5 (inclusive) for ncp is recom- */
  1847. /* mended unless there are evidences that dictate otherwise. */
  1848. /* the lun constant in the data initialization statement is the */
  1849. /* logical unit number of the standard output unit and is, */
  1850. /* therefore, system dependent. */
  1851. /* this subroutine calls the idcldp, idgrid, idpdrv, idptip, and */
  1852. /* idtang subroutines. */
  1853. /* declaration statements */
  1854.     /* Parameter adjustments */
  1855.     --wk;
  1856.     --iwk;
  1857.     --zi;
  1858.     --yi;
  1859.     --xi;
  1860.     --zd;
  1861.     --yd;
  1862.     --xd;
  1863.  
  1864.     /* Function Body */
  1865. /* setting of some input parameters to local variables. */
  1866. /* (for md=1,2,3) */
  1867. /* L10: */
  1868.     md0 = *md;
  1869.     ncp0 = *ncp;
  1870.     ndp0 = *ndp;
  1871.     nxi0 = *nxi;
  1872.     nyi0 = *nyi;
  1873. /* error check.  (for md=1,2,3) */
  1874. /* L20: */
  1875.     if (md0 < 1 || md0 > 3) {
  1876.     goto L90;
  1877.     }
  1878.     if (ncp0 < 2 || ncp0 >= ndp0) {
  1879.     goto L90;
  1880.     }
  1881.     if (ndp0 < 4) {
  1882.     goto L90;
  1883.     }
  1884.     if (nxi0 < 1 || nyi0 < 1) {
  1885.     goto L90;
  1886.     }
  1887.     if (md0 >= 2) {
  1888.     goto L21;
  1889.     }
  1890.     iwk[1] = ncp0;
  1891.     iwk[2] = ndp0;
  1892.     goto L22;
  1893. L21:
  1894.     ncppv = iwk[1];
  1895.     ndppv = iwk[2];
  1896.     if (ncp0 != ncppv) {
  1897.     goto L90;
  1898.     }
  1899.     if (ndp0 != ndppv) {
  1900.     goto L90;
  1901.     }
  1902. L22:
  1903.     if (md0 >= 3) {
  1904.     goto L23;
  1905.     }
  1906.     iwk[3] = nxi0;
  1907.     iwk[4] = nyi0;
  1908.     goto L30;
  1909. L23:
  1910.     nxipv = iwk[3];
  1911.     nyipv = iwk[4];
  1912.     if (nxi0 != nxipv) {
  1913.     goto L90;
  1914.     }
  1915.     if (nyi0 != nyipv) {
  1916.     goto L90;
  1917.     }
  1918. /* allocation of storage areas in the iwk array.  (for md=1,2,3) */
  1919. L30:
  1920.     jwipt = 16;
  1921.     jwiwl = ndp0 * 6 + 1;
  1922.     jwngp0 = jwiwl - 1;
  1923.     jwipl = ndp0 * 24 + 1;
  1924.     jwiwp = ndp0 * 30 + 1;
  1925.     jwipc = ndp0 * 27 + 1;
  1926. /* Computing MAX */
  1927.     i_1 = 31, i_2 = ncp0 + 27;
  1928.     jwigp0 = max(i_1,i_2) * ndp0;
  1929. /* triangulates the x-y plane.  (for md=1) */
  1930. /* L40: */
  1931.     if (md0 > 1) {
  1932.     goto L50;
  1933.     }
  1934.     idtang_(&ndp0, &xd[1], &yd[1], &nt, &iwk[jwipt], &nl, &iwk[jwipl], &iwk[
  1935.         jwiwl], &iwk[jwiwp], &wk[1]);
  1936.     iwk[5] = nt;
  1937.     iwk[6] = nl;
  1938.     if (nt == 0) {
  1939.     return 0;
  1940.     }
  1941. /* determines ncp points closest to each data point.  (for md=1) */
  1942. L50:
  1943.     if (md0 > 1) {
  1944.     goto L60;
  1945.     }
  1946.     idcldp_(&ndp0, &xd[1], &yd[1], &ncp0, &iwk[jwipc]);
  1947.     if (iwk[jwipc] == 0) {
  1948.     return 0;
  1949.     }
  1950. /* sorts output grid points in ascending order of the triangle */
  1951. /* number and the border line segment number.  (for md=1,2) */
  1952. L60:
  1953.     if (md0 == 3) {
  1954.     goto L70;
  1955.     }
  1956.     idgrid_(&xd[1], &yd[1], &nt, &iwk[jwipt], &nl, &iwk[jwipl], &nxi0, &nyi0, 
  1957.         &xi[1], &yi[1], &iwk[jwngp0 + 1], &iwk[jwigp0 + 1]);
  1958. /* estimates partial derivatives at all data points. */
  1959. /* (for md=1,2,3) */
  1960. L70:
  1961.     idpdrv_(&ndp0, &xd[1], &yd[1], &zd[1], &ncp0, &iwk[jwipc], &wk[1]);
  1962. /* interpolates the zi values.  (for md=1,2,3) */
  1963. /* L80: */
  1964.     idpi_1.itpv = 0;
  1965.     jig0mx = 0;
  1966.     jig1mn = nxi0 * nyi0 + 1;
  1967.     nngp = nt + (nl << 1);
  1968.     i_1 = nngp;
  1969.     for (jngp = 1; jngp <= i_1; ++jngp) {
  1970.     iti = jngp;
  1971.     if (jngp <= nt) {
  1972.         goto L81;
  1973.     }
  1974.     il1 = (jngp - nt + 1) / 2;
  1975.     il2 = (jngp - nt + 2) / 2;
  1976.     if (il2 > nl) {
  1977.         il2 = 1;
  1978.     }
  1979.     iti = il1 * (nt + nl) + il2;
  1980. L81:
  1981.     jwngp = jwngp0 + jngp;
  1982.     ngp0 = iwk[jwngp];
  1983.     if (ngp0 == 0) {
  1984.         goto L86;
  1985.     }
  1986.     jig0mn = jig0mx + 1;
  1987.     jig0mx += ngp0;
  1988.     i_2 = jig0mx;
  1989.     for (jigp = jig0mn; jigp <= i_2; ++jigp) {
  1990.         jwigp = jwigp0 + jigp;
  1991.         izi = iwk[jwigp];
  1992.         iyi = (izi - 1) / nxi0 + 1;
  1993.         ixi = izi - nxi0 * (iyi - 1);
  1994.         idptip_(&xd[1], &yd[1], &zd[1], &nt, &iwk[jwipt], &nl, &iwk[jwipl]
  1995.             , &wk[1], &iti, &xi[ixi], &yi[iyi], &zi[izi]);
  1996. /* L82: */
  1997.     }
  1998. L86:
  1999.     jwngp = jwngp0 + (nngp << 1) + 1 - jngp;
  2000.     ngp1 = iwk[jwngp];
  2001.     if (ngp1 == 0) {
  2002.         goto L89;
  2003.     }
  2004.     jig1mx = jig1mn - 1;
  2005.     jig1mn -= ngp1;
  2006.     i_2 = jig1mx;
  2007.     for (jigp = jig1mn; jigp <= i_2; ++jigp) {
  2008.         jwigp = jwigp0 + jigp;
  2009.         izi = iwk[jwigp];
  2010.         iyi = (izi - 1) / nxi0 + 1;
  2011.         ixi = izi - nxi0 * (iyi - 1);
  2012.         idptip_(&xd[1], &yd[1], &zd[1], &nt, &iwk[jwipt], &nl, &iwk[jwipl]
  2013.             , &wk[1], &iti, &xi[ixi], &yi[iyi], &zi[izi]);
  2014. /* L87: */
  2015.     }
  2016. L89:
  2017.     ;}
  2018.     return 0;
  2019. /* error exit */
  2020. L90:
  2021.     err2090_();
  2022.     return 0;
  2023. /* format statement for error message */
  2024. /* L2090: */
  2025. } /* idsfft_ */
  2026.  
  2027. /* Subroutine */ int idtang_(ndp, xd, yd, nt, ipt, nl, ipl, iwl, iwp, wk)
  2028. integer *ndp;
  2029. real *xd, *yd;
  2030. integer *nt, *ipt, *nl, *ipl, *iwl, *iwp;
  2031. real *wk;
  2032. {
  2033.     /* Initialized data */
  2034.  
  2035.     static real ratio = (float)1e-6;
  2036.     static integer nrep = 100;
  2037.  
  2038.     /* System generated locals */
  2039.     integer i_1, i_2, i_3, i_4;
  2040.     real r_1, r_2;
  2041.  
  2042.     /* Local variables */
  2043.     static integer nlfc, ip1p1;
  2044.     static real dsq12, armn;
  2045.     static integer irep;
  2046.     static real dsqi;
  2047.     static integer jp2t3, jp3t3, jpmn;
  2048.     static real dxmn, dymn, xdmp, ydmp, armx;
  2049.     static integer ipti, it1t3, it2t3, jpmx;
  2050.     static real dxmx, dymx;
  2051.     static integer ndpm1, ilft2, iplj1, iplj2, ipmn1, ipmn2, ipti1, ipti2, 
  2052.         nlft2, nlnt3, nsht3, itt3r;
  2053.     extern /* Subroutine */ int err2090_(), err2091_(), err2092_(), err2093_()
  2054.         ;
  2055.     static real dsqmn;
  2056.     static integer ntt3p3;
  2057.     static real dsqmx, x1, y1;
  2058.     static integer jwl1mn;
  2059.     static real ar;
  2060.     static integer ip, jp;
  2061.     static real dx, dy;
  2062.     static integer it;
  2063.     extern integer idxchg_();
  2064.     static integer ip1, ip2, jp1, jp2, ip3, nl0, nt0, ilf, jpc;
  2065.     static real dx21, dy21;
  2066.     static integer nlf, itf[2], nln, nsh, ntf, jwl, its, ndp0, ipl1, ipl2, 
  2067.         jlt3, ipt1, ipt2, ipt3, nlt3, jwl1, itt3, ntt3;
  2068.  
  2069. /* this subroutine performs triangulation.  it divides the x-y */
  2070. /* plane into a number of triangles according to given data */
  2071. /* points in the plane, determines line segments that form the */
  2072. /* border of data area, and determines the triangle numbers */
  2073. /* corresponding to the border line segments. */
  2074. /* at completion, point numbers of the vertexes of each triangle */
  2075. /* are listed counter-clockwise.  point numbers of the end points */
  2076. /* of each border line segment are listed counter-clockwise, */
  2077. /* listing order of the line segments being counter-clockwise. */
  2078. /* the lun constant in the data initialization statement is the */
  2079. /* logical unit number of the standard output unit and is, */
  2080. /* therefore, system dependent. */
  2081. /* this subroutine calls the idxchg function. */
  2082. /* the input parameters are */
  2083. /*     ndp = number of data points, */
  2084. /*     xd  = array of dimension ndp containing the */
  2085. /*           x coordinates of the data points, */
  2086. /*     yd  = array of dimension ndp containing the */
  2087. /*           y coordinates of the data points. */
  2088. /* the output parameters are */
  2089. /*     nt  = number of triangles, */
  2090. /*     ipt = integer array of dimension 6*ndp-15, where the */
  2091. /*           point numbers of the vertexes of the (it)th */
  2092. /*           triangle are to be stored as the (3*it-2)nd, */
  2093. /*           (3*it-1)st, and (3*it)th elements, */
  2094. /*           it=1,2,...,nt, */
  2095. /*     nl  = number of border line segments, */
  2096. /*     ipl = integer array of dimension 6*ndp, where the */
  2097. /*           point numbers of the end points of the (il)th */
  2098. /*           border line segment and its respective triangle */
  2099. /*           number are to be stored as the (3*il-2)nd, */
  2100. /*           (3*il-1)st, and (3*il)th elements, */
  2101. /*           il=1,2,..., nl. */
  2102. /* the other parameters are */
  2103. /*     iwl = integer array of dimension 18*ndp used */
  2104. /*           internally as a work area, */
  2105. /*     iwp = integer array of dimension ndp used */
  2106. /*           internally as a work area, */
  2107. /*     wk  = array of dimension ndp used internally as a */
  2108. /*           work area. */
  2109. /* declaration statements */
  2110.     /* Parameter adjustments */
  2111.     --wk;
  2112.     --iwp;
  2113.     --iwl;
  2114.     --ipl;
  2115.     --ipt;
  2116.     --yd;
  2117.     --xd;
  2118.  
  2119.     /* Function Body */
  2120. /* statement functions */
  2121. /* preliminary processing */
  2122. /* L10: */
  2123.     ndp0 = *ndp;
  2124.     ndpm1 = ndp0 - 1;
  2125.     if (ndp0 < 4) {
  2126.     goto L90;
  2127.     }
  2128. /* determines the closest pair of data points and their midpoint. */
  2129. /* L20: */
  2130. /* Computing 2nd power */
  2131.     r_1 = xd[2] - xd[1];
  2132. /* Computing 2nd power */
  2133.     r_2 = yd[2] - yd[1];
  2134.     dsqmn = r_1 * r_1 + r_2 * r_2;
  2135.     ipmn1 = 1;
  2136.     ipmn2 = 2;
  2137.     i_1 = ndpm1;
  2138.     for (ip1 = 1; ip1 <= i_1; ++ip1) {
  2139.     x1 = xd[ip1];
  2140.     y1 = yd[ip1];
  2141.     ip1p1 = ip1 + 1;
  2142.     i_2 = ndp0;
  2143.     for (ip2 = ip1p1; ip2 <= i_2; ++ip2) {
  2144. /* Computing 2nd power */
  2145.         r_1 = xd[ip2] - x1;
  2146. /* Computing 2nd power */
  2147.         r_2 = yd[ip2] - y1;
  2148.         dsqi = r_1 * r_1 + r_2 * r_2;
  2149.         if (dsqi == (float)0.) {
  2150.         goto L91;
  2151.         }
  2152.         if (dsqi >= dsqmn) {
  2153.         goto L21;
  2154.         }
  2155.         dsqmn = dsqi;
  2156.         ipmn1 = ip1;
  2157.         ipmn2 = ip2;
  2158. L21:
  2159.     ;}
  2160. /* L22: */
  2161.     }
  2162.     dsq12 = dsqmn;
  2163.     xdmp = (xd[ipmn1] + xd[ipmn2]) / (float)2.;
  2164.     ydmp = (yd[ipmn1] + yd[ipmn2]) / (float)2.;
  2165. /* sorts the other (ndp-2) data points in ascending order of */
  2166. /* distance from the midpoint and stores the sorted data point */
  2167. /* numbers in the iwp array. */
  2168. /* L30: */
  2169.     jp1 = 2;
  2170.     i_1 = ndp0;
  2171.     for (ip1 = 1; ip1 <= i_1; ++ip1) {
  2172.     if (ip1 == ipmn1 || ip1 == ipmn2) {
  2173.         goto L31;
  2174.     }
  2175.     ++jp1;
  2176.     iwp[jp1] = ip1;
  2177. /* Computing 2nd power */
  2178.     r_1 = xd[ip1] - xdmp;
  2179. /* Computing 2nd power */
  2180.     r_2 = yd[ip1] - ydmp;
  2181.     wk[jp1] = r_1 * r_1 + r_2 * r_2;
  2182. L31:
  2183.     ;}
  2184.     i_1 = ndpm1;
  2185.     for (jp1 = 3; jp1 <= i_1; ++jp1) {
  2186.     dsqmn = wk[jp1];
  2187.     jpmn = jp1;
  2188.     i_2 = ndp0;
  2189.     for (jp2 = jp1; jp2 <= i_2; ++jp2) {
  2190.         if (wk[jp2] >= dsqmn) {
  2191.         goto L32;
  2192.         }
  2193.         dsqmn = wk[jp2];
  2194.         jpmn = jp2;
  2195. L32:
  2196.     ;}
  2197.     its = iwp[jp1];
  2198.     iwp[jp1] = iwp[jpmn];
  2199.     iwp[jpmn] = its;
  2200.     wk[jpmn] = wk[jp1];
  2201. /* L33: */
  2202.     }
  2203. /* if necessary, modifies the ordering in such a way that the */
  2204. /* first three data points are not collinear. */
  2205. /* L35: */
  2206.     ar = dsq12 * ratio;
  2207.     x1 = xd[ipmn1];
  2208.     y1 = yd[ipmn1];
  2209.     dx21 = xd[ipmn2] - x1;
  2210.     dy21 = yd[ipmn2] - y1;
  2211.     i_1 = ndp0;
  2212.     for (jp = 3; jp <= i_1; ++jp) {
  2213.     ip = iwp[jp];
  2214.     if ((r_1 = (yd[ip] - y1) * dx21 - (xd[ip] - x1) * dy21, dabs(r_1)) > 
  2215.         ar) {
  2216.         goto L37;
  2217.     }
  2218. /* L36: */
  2219.     }
  2220.     goto L92;
  2221. L37:
  2222.     if (jp == 3) {
  2223.     goto L40;
  2224.     }
  2225.     jpmx = jp;
  2226.     jp = jpmx + 1;
  2227.     i_1 = jpmx;
  2228.     for (jpc = 4; jpc <= i_1; ++jpc) {
  2229.     --jp;
  2230.     iwp[jp] = iwp[jp - 1];
  2231. /* L38: */
  2232.     }
  2233.     iwp[3] = ip;
  2234. /* forms the first triangle.  stores point numbers of the ver- */
  2235. /* texes of the triangle in the ipt array, and stores point num- */
  2236. /* bers of the border line segments and the triangle number in */
  2237. /* the ipl array. */
  2238. L40:
  2239.     ip1 = ipmn1;
  2240.     ip2 = ipmn2;
  2241.     ip3 = iwp[3];
  2242.     if ((yd[ip3] - yd[ip1]) * (xd[ip2] - xd[ip1]) - (xd[ip3] - xd[ip1]) * (yd[
  2243.         ip2] - yd[ip1]) >= (float)0.) {
  2244.     goto L41;
  2245.     }
  2246.     ip1 = ipmn2;
  2247.     ip2 = ipmn1;
  2248. L41:
  2249.     nt0 = 1;
  2250.     ntt3 = 3;
  2251.     ipt[1] = ip1;
  2252.     ipt[2] = ip2;
  2253.     ipt[3] = ip3;
  2254.     nl0 = 3;
  2255.     nlt3 = 9;
  2256.     ipl[1] = ip1;
  2257.     ipl[2] = ip2;
  2258.     ipl[3] = 1;
  2259.     ipl[4] = ip2;
  2260.     ipl[5] = ip3;
  2261.     ipl[6] = 1;
  2262.     ipl[7] = ip3;
  2263.     ipl[8] = ip1;
  2264.     ipl[9] = 1;
  2265. /* adds the remaining (ndp-3) data points, one by one. */
  2266. /* L50: */
  2267.     i_1 = ndp0;
  2268.     for (jp1 = 4; jp1 <= i_1; ++jp1) {
  2269.     ip1 = iwp[jp1];
  2270.     x1 = xd[ip1];
  2271.     y1 = yd[ip1];
  2272. /* - determines the visible border line segments. */
  2273.     ip2 = ipl[1];
  2274.     jpmn = 1;
  2275.     dxmn = xd[ip2] - x1;
  2276.     dymn = yd[ip2] - y1;
  2277. /* Computing 2nd power */
  2278.     r_1 = dxmn;
  2279. /* Computing 2nd power */
  2280.     r_2 = dymn;
  2281.     dsqmn = r_1 * r_1 + r_2 * r_2;
  2282.     armn = dsqmn * ratio;
  2283.     jpmx = 1;
  2284.     dxmx = dxmn;
  2285.     dymx = dymn;
  2286.     dsqmx = dsqmn;
  2287.     armx = armn;
  2288.     i_2 = nl0;
  2289.     for (jp2 = 2; jp2 <= i_2; ++jp2) {
  2290.         ip2 = ipl[jp2 * 3 - 2];
  2291.         dx = xd[ip2] - x1;
  2292.         dy = yd[ip2] - y1;
  2293.         ar = dy * dxmn - dx * dymn;
  2294.         if (ar > armn) {
  2295.         goto L51;
  2296.         }
  2297. /* Computing 2nd power */
  2298.         r_1 = dx;
  2299. /* Computing 2nd power */
  2300.         r_2 = dy;
  2301.         dsqi = r_1 * r_1 + r_2 * r_2;
  2302.         if (ar >= -(doublereal)armn && dsqi >= dsqmn) {
  2303.         goto L51;
  2304.         }
  2305.         jpmn = jp2;
  2306.         dxmn = dx;
  2307.         dymn = dy;
  2308.         dsqmn = dsqi;
  2309.         armn = dsqmn * ratio;
  2310. L51:
  2311.         ar = dy * dxmx - dx * dymx;
  2312.         if (ar < -(doublereal)armx) {
  2313.         goto L52;
  2314.         }
  2315. /* Computing 2nd power */
  2316.         r_1 = dx;
  2317. /* Computing 2nd power */
  2318.         r_2 = dy;
  2319.         dsqi = r_1 * r_1 + r_2 * r_2;
  2320.         if (ar <= armx && dsqi >= dsqmx) {
  2321.         goto L52;
  2322.         }
  2323.         jpmx = jp2;
  2324.         dxmx = dx;
  2325.         dymx = dy;
  2326.         dsqmx = dsqi;
  2327.         armx = dsqmx * ratio;
  2328. L52:
  2329.     ;}
  2330.     if (jpmx < jpmn) {
  2331.         jpmx += nl0;
  2332.     }
  2333.     nsh = jpmn - 1;
  2334.     if (nsh <= 0) {
  2335.         goto L60;
  2336.     }
  2337. /* - shifts (rotates) the ipl array to have the invisible border */
  2338. /* - line segments contained in the first part of the ipl array. */
  2339.     nsht3 = nsh * 3;
  2340.     i_2 = nsht3;
  2341.     for (jp2t3 = 3; jp2t3 <= i_2; jp2t3 += 3) {
  2342.         jp3t3 = jp2t3 + nlt3;
  2343.         ipl[jp3t3 - 2] = ipl[jp2t3 - 2];
  2344.         ipl[jp3t3 - 1] = ipl[jp2t3 - 1];
  2345.         ipl[jp3t3] = ipl[jp2t3];
  2346. /* L53: */
  2347.     }
  2348.     i_2 = nlt3;
  2349.     for (jp2t3 = 3; jp2t3 <= i_2; jp2t3 += 3) {
  2350.         jp3t3 = jp2t3 + nsht3;
  2351.         ipl[jp2t3 - 2] = ipl[jp3t3 - 2];
  2352.         ipl[jp2t3 - 1] = ipl[jp3t3 - 1];
  2353.         ipl[jp2t3] = ipl[jp3t3];
  2354. /* L54: */
  2355.     }
  2356.     jpmx -= nsh;
  2357. /* - adds triangles to the ipt array, updates border line */
  2358. /* - segments in the ipl array, and sets flags for the border */
  2359. /* - line segments to be reexamined in the iwl array. */
  2360. L60:
  2361.     jwl = 0;
  2362.     i_2 = nl0;
  2363.     for (jp2 = jpmx; jp2 <= i_2; ++jp2) {
  2364.         jp2t3 = jp2 * 3;
  2365.         ipl1 = ipl[jp2t3 - 2];
  2366.         ipl2 = ipl[jp2t3 - 1];
  2367.         it = ipl[jp2t3];
  2368. /* - - adds a triangle to the ipt array. */
  2369.         ++nt0;
  2370.         ntt3 += 3;
  2371.         ipt[ntt3 - 2] = ipl2;
  2372.         ipt[ntt3 - 1] = ipl1;
  2373.         ipt[ntt3] = ip1;
  2374. /* - - updates border line segments in the ipl array. */
  2375.         if (jp2 != jpmx) {
  2376.         goto L61;
  2377.         }
  2378.         ipl[jp2t3 - 1] = ip1;
  2379.         ipl[jp2t3] = nt0;
  2380. L61:
  2381.         if (jp2 != nl0) {
  2382.         goto L62;
  2383.         }
  2384.         nln = jpmx + 1;
  2385.         nlnt3 = nln * 3;
  2386.         ipl[nlnt3 - 2] = ip1;
  2387.         ipl[nlnt3 - 1] = ipl[1];
  2388.         ipl[nlnt3] = nt0;
  2389. /* - - determines the vertex that does not lie on the border */
  2390. /* - - line segments. */
  2391. L62:
  2392.         itt3 = it * 3;
  2393.         ipti = ipt[itt3 - 2];
  2394.         if (ipti != ipl1 && ipti != ipl2) {
  2395.         goto L63;
  2396.         }
  2397.         ipti = ipt[itt3 - 1];
  2398.         if (ipti != ipl1 && ipti != ipl2) {
  2399.         goto L63;
  2400.         }
  2401.         ipti = ipt[itt3];
  2402. /* - - checks if the exchange is necessary. */
  2403. L63:
  2404.         if (idxchg_(&xd[1], &yd[1], &ip1, &ipti, &ipl1, &ipl2) == 0) {
  2405.         goto L64;
  2406.         }
  2407. /* - - modifies the ipt array when necessary. */
  2408.         ipt[itt3 - 2] = ipti;
  2409.         ipt[itt3 - 1] = ipl1;
  2410.         ipt[itt3] = ip1;
  2411.         ipt[ntt3 - 1] = ipti;
  2412.         if (jp2 == jpmx) {
  2413.         ipl[jp2t3] = it;
  2414.         }
  2415.         if (jp2 == nl0 && ipl[3] == it) {
  2416.         ipl[3] = nt0;
  2417.         }
  2418. /* - - sets flags in the iwl array. */
  2419.         jwl += 4;
  2420.         iwl[jwl - 3] = ipl1;
  2421.         iwl[jwl - 2] = ipti;
  2422.         iwl[jwl - 1] = ipti;
  2423.         iwl[jwl] = ipl2;
  2424. L64:
  2425.     ;}
  2426.     nl0 = nln;
  2427.     nlt3 = nlnt3;
  2428.     nlf = jwl / 2;
  2429.     if (nlf == 0) {
  2430.         goto L79;
  2431.     }
  2432. /* - improves triangulation. */
  2433. /* L70: */
  2434.     ntt3p3 = ntt3 + 3;
  2435.     i_2 = nrep;
  2436.     for (irep = 1; irep <= i_2; ++irep) {
  2437.         i_3 = nlf;
  2438.         for (ilf = 1; ilf <= i_3; ++ilf) {
  2439.         ilft2 = ilf << 1;
  2440.         ipl1 = iwl[ilft2 - 1];
  2441.         ipl2 = iwl[ilft2];
  2442. /* - - locates in the ipt array two triangles on both sides 
  2443. of */
  2444. /* - - the flagged line segment. */
  2445.         ntf = 0;
  2446.         i_4 = ntt3;
  2447.         for (itt3r = 3; itt3r <= i_4; itt3r += 3) {
  2448.             itt3 = ntt3p3 - itt3r;
  2449.             ipt1 = ipt[itt3 - 2];
  2450.             ipt2 = ipt[itt3 - 1];
  2451.             ipt3 = ipt[itt3];
  2452.             if (ipl1 != ipt1 && ipl1 != ipt2 && ipl1 != ipt3) {
  2453.             goto L71;
  2454.             }
  2455.             if (ipl2 != ipt1 && ipl2 != ipt2 && ipl2 != ipt3) {
  2456.             goto L71;
  2457.             }
  2458.             ++ntf;
  2459.             itf[ntf - 1] = itt3 / 3;
  2460.             if (ntf == 2) {
  2461.             goto L72;
  2462.             }
  2463. L71:
  2464.         ;}
  2465.         if (ntf < 2) {
  2466.             goto L76;
  2467.         }
  2468. /* - - determines the vertexes of the triangles that do not 
  2469. lie */
  2470. /* - - on the line segment. */
  2471. L72:
  2472.         it1t3 = itf[0] * 3;
  2473.         ipti1 = ipt[it1t3 - 2];
  2474.         if (ipti1 != ipl1 && ipti1 != ipl2) {
  2475.             goto L73;
  2476.         }
  2477.         ipti1 = ipt[it1t3 - 1];
  2478.         if (ipti1 != ipl1 && ipti1 != ipl2) {
  2479.             goto L73;
  2480.         }
  2481.         ipti1 = ipt[it1t3];
  2482. L73:
  2483.         it2t3 = itf[1] * 3;
  2484.         ipti2 = ipt[it2t3 - 2];
  2485.         if (ipti2 != ipl1 && ipti2 != ipl2) {
  2486.             goto L74;
  2487.         }
  2488.         ipti2 = ipt[it2t3 - 1];
  2489.         if (ipti2 != ipl1 && ipti2 != ipl2) {
  2490.             goto L74;
  2491.         }
  2492.         ipti2 = ipt[it2t3];
  2493. /* - - checks if the exchange is necessary. */
  2494. L74:
  2495.         if (idxchg_(&xd[1], &yd[1], &ipti1, &ipti2, &ipl1, &ipl2) == 
  2496.             0) {
  2497.             goto L76;
  2498.         }
  2499. /* - - modifies the ipt array when necessary. */
  2500.         ipt[it1t3 - 2] = ipti1;
  2501.         ipt[it1t3 - 1] = ipti2;
  2502.         ipt[it1t3] = ipl1;
  2503.         ipt[it2t3 - 2] = ipti2;
  2504.         ipt[it2t3 - 1] = ipti1;
  2505.         ipt[it2t3] = ipl2;
  2506. /* - - sets new flags. */
  2507.         jwl += 8;
  2508.         iwl[jwl - 7] = ipl1;
  2509.         iwl[jwl - 6] = ipti1;
  2510.         iwl[jwl - 5] = ipti1;
  2511.         iwl[jwl - 4] = ipl2;
  2512.         iwl[jwl - 3] = ipl2;
  2513.         iwl[jwl - 2] = ipti2;
  2514.         iwl[jwl - 1] = ipti2;
  2515.         iwl[jwl] = ipl1;
  2516.         i_4 = nlt3;
  2517.         for (jlt3 = 3; jlt3 <= i_4; jlt3 += 3) {
  2518.             iplj1 = ipl[jlt3 - 2];
  2519.             iplj2 = ipl[jlt3 - 1];
  2520.             if (iplj1 == ipl1 && iplj2 == ipti2 || iplj2 == ipl1 && 
  2521.                 iplj1 == ipti2) {
  2522.             ipl[jlt3] = itf[0];
  2523.             }
  2524.             if (iplj1 == ipl2 && iplj2 == ipti1 || iplj2 == ipl2 && 
  2525.                 iplj1 == ipti1) {
  2526.             ipl[jlt3] = itf[1];
  2527.             }
  2528. /* L75: */
  2529.         }
  2530. L76:
  2531.         ;}
  2532.         nlfc = nlf;
  2533.         nlf = jwl / 2;
  2534.         if (nlf == nlfc) {
  2535.         goto L79;
  2536.         }
  2537. /* - - resets the iwl array for the next round. */
  2538.         jwl = 0;
  2539.         jwl1mn = nlfc + 1 << 1;
  2540.         nlft2 = nlf << 1;
  2541.         i_3 = nlft2;
  2542.         for (jwl1 = jwl1mn; jwl1 <= i_3; jwl1 += 2) {
  2543.         jwl += 2;
  2544.         iwl[jwl - 1] = iwl[jwl1 - 1];
  2545.         iwl[jwl] = iwl[jwl1];
  2546. /* L77: */
  2547.         }
  2548.         nlf = jwl / 2;
  2549. /* L78: */
  2550.     }
  2551. L79:
  2552.     ;}
  2553. /* rearranges the ipt array so that the vertexes of each triangle */
  2554. /* are listed counter-clockwise. */
  2555. /* L80: */
  2556.     i_1 = ntt3;
  2557.     for (itt3 = 3; itt3 <= i_1; itt3 += 3) {
  2558.     ip1 = ipt[itt3 - 2];
  2559.     ip2 = ipt[itt3 - 1];
  2560.     ip3 = ipt[itt3];
  2561.     if ((yd[ip3] - yd[ip1]) * (xd[ip2] - xd[ip1]) - (xd[ip3] - xd[ip1]) * 
  2562.         (yd[ip2] - yd[ip1]) >= (float)0.) {
  2563.         goto L81;
  2564.     }
  2565.     ipt[itt3 - 2] = ip2;
  2566.     ipt[itt3 - 1] = ip1;
  2567. L81:
  2568.     ;}
  2569.     *nt = nt0;
  2570.     *nl = nl0;
  2571.     return 0;
  2572. /* error exit */
  2573. L90:
  2574.     err2090_();
  2575.     goto L93;
  2576. L91:
  2577.     err2091b_();
  2578.     goto L93;
  2579. L92:
  2580.     err2092_();
  2581. L93:
  2582.     err2093_();
  2583.     *nt = 0;
  2584.     return 0;
  2585. /* format statements */
  2586. /* L2090: */
  2587. /* L2091: */
  2588. /* L2092: */
  2589. /* L2093: */
  2590. } /* idtang_ */
  2591.  
  2592. integer idxchg_(x, y, i1, i2, i3, i4)
  2593. real *x, *y;
  2594. integer *i1, *i2, *i3, *i4;
  2595. {
  2596.     /* System generated locals */
  2597.     integer ret_val;
  2598.     real r_1, r_2;
  2599.     static real equiv_0[1], equiv_1[1], equiv_2[1], equiv_3[1], equiv_4[1], 
  2600.         equiv_5[1];
  2601.  
  2602.     /* Local variables */
  2603.     static real u1, u2, u3, x1, y1, x2, y2, x3, y3, x4, y4, u4;
  2604.     static integer idx;
  2605. #define a1sq (equiv_2)
  2606. #define b1sq (equiv_3)
  2607. #define c1sq (equiv_0)
  2608. #define c2sq (equiv_0)
  2609. #define a3sq (equiv_1)
  2610. #define b2sq (equiv_1)
  2611. #define b3sq (equiv_2)
  2612. #define a4sq (equiv_3)
  2613. #define b4sq (equiv_4)
  2614. #define a2sq (equiv_4)
  2615. #define c4sq (equiv_5)
  2616. #define c3sq (equiv_5)
  2617.     static real s1sq, s2sq, s3sq, s4sq;
  2618.  
  2619. /* this function determines whether or not the exchange of two */
  2620. /* triangles is necessary on the basis of max-min-angle criterion */
  2621. /* by c. l. lawson. */
  2622. /* the input parameters are */
  2623. /*     x,y = arrays containing the coordinates of the data */
  2624. /*           points, */
  2625. /*     i1,i2,i3,i4 = point numbers of four points p1, p2, */
  2626. /*           p3, and p4 that form a quadrilateral with p3 */
  2627. /*           and p4 connected diagonally. */
  2628. /* this function returns an integer value 1 (one) when an ex- */
  2629. /* change is necessary, and 0 (zero) otherwise. */
  2630. /* declaration statements */
  2631. /* preliminary processing */
  2632.     /* Parameter adjustments */
  2633.     --y;
  2634.     --x;
  2635.  
  2636.     /* Function Body */
  2637. /* L10: */
  2638.     x1 = x[*i1];
  2639.     y1 = y[*i1];
  2640.     x2 = x[*i2];
  2641.     y2 = y[*i2];
  2642.     x3 = x[*i3];
  2643.     y3 = y[*i3];
  2644.     x4 = x[*i4];
  2645.     y4 = y[*i4];
  2646. /* calculation */
  2647. /* L20: */
  2648.     idx = 0;
  2649.     u3 = (y2 - y3) * (x1 - x3) - (x2 - x3) * (y1 - y3);
  2650.     u4 = (y1 - y4) * (x2 - x4) - (x1 - x4) * (y2 - y4);
  2651.     if (u3 * u4 <= (float)0.) {
  2652.     goto L30;
  2653.     }
  2654.     u1 = (y3 - y1) * (x4 - x1) - (x3 - x1) * (y4 - y1);
  2655.     u2 = (y4 - y2) * (x3 - x2) - (x4 - x2) * (y3 - y2);
  2656. /* Computing 2nd power */
  2657.     r_1 = x1 - x3;
  2658. /* Computing 2nd power */
  2659.     r_2 = y1 - y3;
  2660.     *a1sq = r_1 * r_1 + r_2 * r_2;
  2661. /* Computing 2nd power */
  2662.     r_1 = x4 - x1;
  2663. /* Computing 2nd power */
  2664.     r_2 = y4 - y1;
  2665.     *b1sq = r_1 * r_1 + r_2 * r_2;
  2666. /* Computing 2nd power */
  2667.     r_1 = x3 - x4;
  2668. /* Computing 2nd power */
  2669.     r_2 = y3 - y4;
  2670.     *c1sq = r_1 * r_1 + r_2 * r_2;
  2671. /* Computing 2nd power */
  2672.     r_1 = x2 - x4;
  2673. /* Computing 2nd power */
  2674.     r_2 = y2 - y4;
  2675.     *a2sq = r_1 * r_1 + r_2 * r_2;
  2676. /* Computing 2nd power */
  2677.     r_1 = x3 - x2;
  2678. /* Computing 2nd power */
  2679.     r_2 = y3 - y2;
  2680.     *b2sq = r_1 * r_1 + r_2 * r_2;
  2681. /* Computing 2nd power */
  2682.     r_1 = x2 - x1;
  2683. /* Computing 2nd power */
  2684.     r_2 = y2 - y1;
  2685.     *c3sq = r_1 * r_1 + r_2 * r_2;
  2686.     s1sq = u1 * u1 / (*c1sq * dmax(*a1sq,*b1sq));
  2687.     s2sq = u2 * u2 / (*c2sq * dmax(*a2sq,*b2sq));
  2688.     s3sq = u3 * u3 / (*c3sq * dmax(*a3sq,*b3sq));
  2689.     s4sq = u4 * u4 / (*c4sq * dmax(*a4sq,*b4sq));
  2690.     if (dmin(s1sq,s2sq) < dmin(s3sq,s4sq)) {
  2691.     idx = 1;
  2692.     }
  2693. L30:
  2694.     ret_val = idx;
  2695.     return ret_val;
  2696. } /* idxchg_ */
  2697.  
  2698. #undef c3sq
  2699. #undef c4sq
  2700. #undef a2sq
  2701. #undef b4sq
  2702. #undef a4sq
  2703. #undef b3sq
  2704. #undef b2sq
  2705. #undef a3sq
  2706. #undef c2sq
  2707. #undef c1sq
  2708. #undef b1sq
  2709. #undef a1sq
  2710.  
  2711.  
  2712.