home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / libs / gle / util / contour / fitcf.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-29  |  10.5 KB  |  546 lines
  1. /* fitcf.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. #define gprint printf
  8. /* --   FITCF */
  9. /* Subroutine */ int glefitcf_(mode, x, y, l, m, u, v, n)
  10. integer *mode;
  11. real *x, *y;
  12. integer *l, *m;
  13. real *u, *v;
  14. integer *n;
  15. {
  16.     /* System generated locals */
  17.     integer i_1, i_2;
  18.     real r_1;
  19.     static real equiv_1[1], equiv_2[1], equiv_3[1], equiv_4[1], equiv_5[1],
  20.         equiv_6[1], equiv_7[1], equiv_8[1], equiv_9[1], equiv_10[1],
  21.         equiv_11[1], equiv_12[1];
  22.  
  23.     /* Builtin functions */
  24.     double sqrt();
  25.  
  26.     /* Local variables */
  27.     static integer mode0, i, j, k;
  28. #define r (equiv_10)
  29. #define z (equiv_10)
  30. #define a1 (equiv_7)
  31. #define a2 (equiv_8)
  32.     static real a3, a4;
  33. #define b1 (equiv_9)
  34. #define b2 (equiv_1)
  35. #define b3 (equiv_2)
  36. #define b4 (equiv_3)
  37.     static integer l0, m0, n0;
  38. #define m1 (equiv_9)
  39.     static integer k5;
  40. #define m2 (equiv_1)
  41. #define m3 (equiv_2)
  42. #define m4 (equiv_3)
  43. #define p0 (equiv_5)
  44.     static real p1;
  45. #define p2 (equiv_7)
  46. #define p3 (equiv_9)
  47. #define q0 (equiv_6)
  48. #define q1 (equiv_4)
  49. #define q2 (equiv_11)
  50. #define q3 (equiv_12)
  51. #define t2 (equiv_4)
  52.     static real t3;
  53. #define w2 (equiv_11)
  54. #define w3 (equiv_12)
  55. #define x2 (equiv_5)
  56.     static real x3;
  57.     static integer modem1;
  58.     static real x4, x5;
  59. #define y2 (equiv_6)
  60.     static real y3, y4, y5;
  61.     extern doublereal gutre2_();
  62. #define dz (equiv_8)
  63.     static real rm;
  64. #define sw (equiv_10)
  65.     extern /* Subroutine */ int gd_message__();
  66.     static integer ierval[1], lm1, mm1;
  67.     static real cos2, cos3, sin2, sin3;
  68.  
  69. /*     (Smooth Curve Fitting) */
  70. /*     This subroutine fits a smooth curve to a given set of input */
  71. /*     data points in  an X-Y  plane.  It  interpolates points  in */
  72. /*     each interval between a pair of data points and generates a */
  73. /*     set of output  points consisting of  the input data  points */
  74. /*     and the  interpolated  points.   It can  process  either  a */
  75. /*     single-valued function or a multiple-valued function. */
  76.  
  77. /*     The input arguments are: */
  78.  
  79. /*     MODE = mode of the curve (must be 1 or 2) */
  80. /*          = 1 for a single-valued function */
  81. /*          = 2 for multiple-valued function */
  82. /*     X  = Array of  dimension L storing  the abscissas  of input */
  83. /*          data points (in ascending or descending order for mode */
  84. /*          = 1) */
  85. /*     Y  = Array of  dimension L storing  the ordinates  of input */
  86. /*          data points */
  87. /*     L  = Number of input data points (must be 2 or greater) */
  88. /*     M  = Number of subintervals between each pair of input data */
  89. /*          points (must be 2 or greater). */
  90. /*     N  = Number of output points */
  91. /*        = (L-1)*M+1 */
  92.  
  93. /*     The output arguments are: */
  94.  
  95. /*     U  = Array of dimension N where the abscissas of output */
  96. /*          points are to be displayed */
  97. /*     V  = Array of dimension N where the ordinates of output */
  98. /*          points are to be displayed */
  99.  
  100. /*     Author:  Hiroshi Akima,  "Interpolation  and  Smooth  Curve */
  101. /*              Fitting Based on Local Procedures", COMM.  ACM 15, */
  102. /*              914-918 (1972), and "A New Method of Interpolation */
  103. /*              and  Smooth   Curve   Fitting   Based   on   Local */
  104. /*              Procedures", J. ACM 17, 589-602 (1970). */
  105.  
  106. /*     Corrections: M.R. Andersen, "Remark on Algorithm 433", ACM */
  107. /*                  Trans. on Math. Software, 2, 208 (1976). */
  108. /*     (30-JAN-82) */
  109.  
  110. /*     Preliminary processing */
  111.  
  112.     /* Parameter adjustments */
  113.     --v;
  114.     --u;
  115.     --y;
  116.     --x;
  117.  
  118.     /* Function Body */
  119.     mode0 = *mode;
  120.     modem1 = mode0 - 1;
  121.     l0 = *l;
  122.     lm1 = l0 - 1;
  123.     m0 = *m;
  124.     mm1 = m0 - 1;
  125.     n0 = *n;
  126.     if (mode0 <= 0) {
  127.     goto L320;
  128.     }
  129.     if (mode0 >= 3) {
  130.     goto L320;
  131.     }
  132.     if (lm1 <= 0) {
  133.     goto L330;
  134.     }
  135.     if (mm1 <= 0) {
  136.     goto L340;
  137.     }
  138.     if (n0 != lm1 * m0 + 1) {
  139.     gprint("Improper n value %ld, wanted %ld \n",n0,lm1 * m0 + 1);
  140.     gprint("n0 %ld, lm1 %ld,  m0 %ld  l0 %ld \n",n0,lm1,m0,l0);
  141.     goto L350;
  142.     }
  143.  
  144.     switch (mode0) {
  145.     case 1:  goto L10;
  146.     case 2:  goto L60;
  147.     }
  148. L10:
  149.     i = 2;
  150.     if ((r_1 = x[1] - x[2]) < (float)0.) {
  151.     goto L20;
  152.     } else if (r_1 == 0) {
  153.     gprint("Identical x values 1, %g %g \n",(double) x[1], (double) x[2]);
  154.     goto L360;
  155.     } else {
  156.     goto L40;
  157.     }
  158. L20:
  159.     if (l0 <= 2) {
  160.     goto L80;
  161.     }
  162.     i_1 = l0;
  163.     for (i = 3; i <= i_1; ++i) {
  164.     if ((r_1 = x[i - 1] - x[i]) < (float)0.) {
  165.         goto L30;
  166.     } else if (r_1 == 0) {
  167.     gprint("Identical x values i, %g %g \n",(int) i,(double) x[i-1], (double) x[i]);
  168.         goto L360;
  169.     } else {
  170.         goto L370;
  171.     }
  172. L30:
  173.     ;}
  174.     goto L80;
  175. L40:
  176.     if (l0 <= 2) {
  177.     goto L80;
  178.     }
  179.     i_1 = l0;
  180.     for (i = 3; i <= i_1; ++i) {
  181.     if ((r_1 = x[i - 1] - x[i]) < (float)0.) {
  182.         goto L370;
  183.     } else if (r_1 == 0) {
  184.     gprint("Identical x aavalues i, %g %g \n",(int) i,(double) x[i-1], (double) x[i]);
  185.         goto L360;
  186.     } else {
  187.         goto L50;
  188.     }
  189. L50:
  190.     ;}
  191.     goto L80;
  192. L60:
  193.     i_1 = l0;
  194.     for (i = 2; i <= i_1; ++i) {
  195.     if (x[i - 1] != x[i]) {
  196.         goto L70;
  197.     }
  198.     if (y[i - 1] == y[i]) {
  199.         goto L380;
  200.     }
  201. L70:
  202.     ;}
  203.  
  204. L80:
  205.     k = n0 + m0;
  206.     i = l0 + 1;
  207.     i_1 = l0;
  208.     for (j = 1; j <= i_1; ++j) {
  209.     k -= m0;
  210.     --i;
  211.     u[k] = x[i];
  212. /* L90: */
  213.     v[k] = y[i];
  214.     }
  215.     rm = (real) m0;
  216.     rm = (float)1. / rm;
  217.  
  218. /*     Main DO-loop */
  219.  
  220.     k5 = m0 + 1;
  221.     i_1 = l0;
  222.     for (i = 1; i <= i_1; ++i) {
  223.  
  224. /*     Routines to  pick  up  necessary  X and  Y  values  and  to */
  225. /*     estimate them if necessary */
  226.  
  227.     if (i > 1) {
  228.         goto L130;
  229.     }
  230.     x3 = u[1];
  231.     y3 = v[1];
  232.     x4 = u[m0 + 1];
  233.     y4 = v[m0 + 1];
  234.     a3 = x4 - x3;
  235.     *b3 = y4 - y3;
  236.     if (modem1 == 0) {
  237.         *m3 = *b3 / a3;
  238.     }
  239.     if (l0 != 2) {
  240.         goto L140;
  241.     }
  242.     a4 = a3;
  243.     *b4 = *b3;
  244. L100:
  245.     switch (mode0) {
  246.         case 1:  goto L120;
  247.         case 2:  goto L110;
  248.     }
  249. L110:
  250.     *a2 = a3 + a3 - a4;
  251.     *a1 = *a2 + *a2 - a3;
  252. L120:
  253.     *b2 = *b3 + *b3 - *b4;
  254.     *b1 = *b2 + *b2 - *b3;
  255.     switch (mode0) {
  256.         case 1:  goto L180;
  257.         case 2:  goto L210;
  258.     }
  259. L130:
  260.     *x2 = x3;
  261.     *y2 = y3;
  262.     x3 = x4;
  263.     y3 = y4;
  264.     x4 = x5;
  265.     y4 = y5;
  266.     *a1 = *a2;
  267.     *b1 = *b2;
  268.     *a2 = a3;
  269.     *b2 = *b3;
  270.     a3 = a4;
  271.     *b3 = *b4;
  272.     if (i >= lm1) {
  273.         goto L150;
  274.     }
  275. L140:
  276.     k5 += m0;
  277.     x5 = u[k5];
  278.     y5 = v[k5];
  279.     a4 = x5 - x4;
  280.     *b4 = y5 - y4;
  281.     if (modem1 == 0) {
  282.         *m4 = *b4 / a4;
  283.     }
  284.     goto L160;
  285. L150:
  286.     if (modem1 != 0) {
  287.         a4 = a3 + a3 - *a2;
  288.     }
  289.     *b4 = *b3 + *b3 - *b2;
  290. L160:
  291.     if (i == 1) {
  292.         goto L100;
  293.     }
  294.     switch (mode0) {
  295.         case 1:  goto L170;
  296.         case 2:  goto L200;
  297.     }
  298.  
  299. /*     Numerical differentiation */
  300.  
  301. L170:
  302.     *t2 = t3;
  303. L180:
  304.     *w2 = (r_1 = *m4 - *m3, dabs(r_1));
  305.     *w3 = (r_1 = *m2 - *m1, dabs(r_1));
  306.     *sw = *w2 + *w3;
  307.     if (*sw != (float)0.) {
  308.         goto L190;
  309.     }
  310.     *w2 = (float).5;
  311.     *w3 = (float).5;
  312.     *sw = (float)1.;
  313. L190:
  314.     t3 = (*w2 * *m2 + *w3 * *m3) / *sw;
  315.     if (i - 1 <= 0) {
  316.         goto L310;
  317.     } else {
  318.         goto L240;
  319.     }
  320.  
  321. L200:
  322.     cos2 = cos3;
  323.     sin2 = sin3;
  324. L210:
  325.     *w2 = (r_1 = a3 * *b4 - a4 * *b3, dabs(r_1));
  326.     *w3 = (r_1 = *a1 * *b2 - *a2 * *b1, dabs(r_1));
  327.     if (*w2 + *w3 != (float)0.) {
  328.         goto L220;
  329.     }
  330.     *w2 = gutre2_(&a3, b3);
  331.     *w3 = gutre2_(a2, b2);
  332. L220:
  333.     cos3 = *w2 * *a2 + *w3 * a3;
  334.     sin3 = *w2 * *b2 + *w3 * *b3;
  335.     *r = cos3 * cos3 + sin3 * sin3;
  336.     if (*r == (float)0.) {
  337.         goto L230;
  338.     }
  339.     *r = sqrt(*r);
  340.     cos3 /= *r;
  341.     sin3 /= *r;
  342. L230:
  343.     if (i - 1 <= 0) {
  344.         goto L310;
  345.     } else {
  346.         goto L250;
  347.     }
  348.  
  349. /*     Determination of the coefficients */
  350.  
  351. L240:
  352.     *q2 = ((*m2 - *t2) * (float)2. + *m2 - t3) / *a2;
  353.     *q3 = (-(doublereal)(*m2) - *m2 + *t2 + t3) / (*a2 * *a2);
  354.     goto L260;
  355.  
  356. L250:
  357.     *r = gutre2_(a2, b2);
  358.     p1 = *r * cos2;
  359.     *p2 = *a2 * (float)3. - *r * (cos2 + cos2 + cos3);
  360.     *p3 = *a2 - p1 - *p2;
  361.     *q1 = *r * sin2;
  362.     *q2 = *b2 * (float)3. - *r * (sin2 + sin2 + sin3);
  363.     *q3 = *b2 - *q1 - *q2;
  364.     goto L280;
  365.  
  366. /*     Computation of the polynomials */
  367.  
  368. L260:
  369.     *dz = *a2 * rm;
  370.     *z = (float)0.;
  371.     i_2 = mm1;
  372.     for (j = 1; j <= i_2; ++j) {
  373.         ++k;
  374.         *z += *dz;
  375.         u[k] = *p0 + *z;
  376. /* L270: */
  377.         v[k] = *q0 + *z * (*q1 + *z * (*q2 + *z * *q3));
  378.     }
  379.     goto L300;
  380.  
  381. L280:
  382.     *z = (float)0.;
  383.     i_2 = mm1;
  384.     for (j = 1; j <= i_2; ++j) {
  385.         ++k;
  386.         *z += rm;
  387.         u[k] = *p0 + *z * (p1 + *z * (*p2 + *z * *p3));
  388. /* L290: */
  389.         v[k] = *q0 + *z * (*q1 + *z * (*q2 + *z * *q3));
  390.     }
  391.  
  392. L300:
  393.     ++k;
  394. L310:
  395.     ;}
  396.     goto L410;
  397.  
  398. /*     Error exit */
  399.  
  400. L320:
  401.     gd_message__("Cannot SMOOTH: Mode out of proper range 1..2", 29L);
  402.     goto L400;
  403. L330:
  404.     gd_message__("Cannot SMOOTH: L = 1 or less", 13L);
  405.     goto L400;
  406. L340:
  407.     gd_message__("Cannot SMOOTH: M = 1 or less", 13L);
  408.     goto L400;
  409. L350:
  410.     gd_message__("Cannot SMOOTH: Improper N value", 16L);
  411.     goto L400;
  412. L360:
  413.     gd_message__("Cannot SMOOTH: Identical X values", 18L);
  414.     goto L390;
  415. L370:
  416.     gd_message__("Cannot SMOOTH: X values out of sequence", 24L);
  417.     goto L390;
  418. L380:
  419.     gd_message__("Cannot SMOOTH: Identical X and Y values", 24L);
  420. L390:
  421.     ierval[0] = i;
  422. L400:
  423.     ierval[0] = mode0;
  424.     ierval[0] = l0;
  425.     ierval[0] = m0;
  426.     ierval[0] = n0;
  427.  
  428. L410:
  429.     return 0;
  430. } /* glefitcf_ */
  431.  
  432. #undef sw
  433. #undef dz
  434. #undef y2
  435. #undef x2
  436. #undef w3
  437. #undef w2
  438. #undef t2
  439. #undef q3
  440. #undef q2
  441. #undef q1
  442. #undef q0
  443. #undef p3
  444. #undef p2
  445. #undef p0
  446. #undef m4
  447. #undef m3
  448. #undef m2
  449. #undef m1
  450. #undef b4
  451. #undef b3
  452. #undef b2
  453. #undef b1
  454. #undef a2
  455. #undef a1
  456. #undef z
  457. #undef r
  458.  
  459.  
  460. /* --   GUTRE2 */
  461. doublereal gutre2_(a, b)
  462. real *a, *b;
  463. {
  464.     /* Initialized data */
  465.  
  466.     static real zero = (float)0.;
  467.     static real two = (float)2.;
  468.     static real four = (float)4.;
  469.  
  470.     /* System generated locals */
  471.     real ret_val, r_1;
  472.  
  473.     /* Local variables */
  474.     static real aabs, babs, r, s;
  475.  
  476. /*     (Euclidean Norm of 2-Vector) */
  477. /*     Return SQRT(A**2+B**2)  without doing  a square  root,  and */
  478. /*     without destructive underflow or overflow. */
  479. /*     (Cleve Moler, University of New Mexico) */
  480. /*     (09-APR-82) */
  481.  
  482.  
  483.     aabs = dabs(*a);
  484.     babs = dabs(*b);
  485. /* --- IF (BABS .GT. AABS) THEN -- SWAP A AND B */
  486.     if (babs <= aabs) {
  487.     goto L10;
  488.     }
  489.     r = babs;
  490.     babs = aabs;
  491.     aabs = r;
  492. /* --- END IF */
  493. /* --- IF (BABS .EQ. ZERO) THEN -- Special case sudden exit */
  494. L10:
  495.     if (babs != zero) {
  496.     goto L20;
  497.     }
  498.     ret_val = aabs;
  499.     goto L40;
  500. /* --- END IF */
  501. /* --- DO FOREVER */
  502. L20:
  503. /* Computing 2nd power */
  504.     r_1 = babs / aabs;
  505.     r = r_1 * r_1;
  506. /* ---      IF ((TWO+R) .EQ. TWO) THEN -- Converged */
  507.     if (two + r != two) {
  508.     goto L30;
  509.     }
  510.     ret_val = aabs;
  511.     goto L40;
  512. /* ---      END IF */
  513. L30:
  514.     s = r / (four + r);
  515.     aabs += two * s * aabs;
  516.     babs = s * babs;
  517. /* --- END FOREVER */
  518.     goto L20;
  519.  
  520. L40:
  521.     return ret_val;
  522. } /* gutre2_ */
  523. gd_message__(char *s,long l)
  524. {
  525.     gprint("%s  %ld \n",s,l);
  526. }
  527.  
  528.  
  529.  
  530.  
  531.  
  532.  
  533.  
  534.  
  535.  
  536.  
  537.  
  538.  
  539.  
  540.  
  541.  
  542.  
  543.  
  544.  
  545.                                  
  546.