home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Photo CD Demo 1
/
Demo.bin
/
gle
/
gle
/
fitcf.c
< prev
next >
Wrap
C/C++ Source or Header
|
1992-11-29
|
11KB
|
546 lines
/* fitcf.f -- translated by f2c (version of 26 February 1990 17:38:00).
You must link the resulting object file with the libraries:
-lF77 -lI77 -lm -lc (in that order)
*/
#include "fitcf.h"
/* -- FITCF */
/* Subroutine */ int glefitcf_(mode, x, y, l, m, u, v, n)
integer *mode;
real *x, *y;
integer *l, *m;
real *u, *v;
integer *n;
{
/* System generated locals */
integer i_1, i_2;
real r_1;
static real equiv_1[1], equiv_2[1], equiv_3[1], equiv_4[1], equiv_5[1],
equiv_6[1], equiv_7[1], equiv_8[1], equiv_9[1], equiv_10[1],
equiv_11[1], equiv_12[1];
/* Builtin functions */
double sqrt();
/* Local variables */
static integer mode0, i, j, k;
#define r (equiv_10)
#define z (equiv_10)
#define a1 (equiv_7)
#define a2 (equiv_8)
static real a3, a4;
#define b1 (equiv_9)
#define b2 (equiv_1)
#define b3 (equiv_2)
#define b4 (equiv_3)
static integer l0, m0, n0;
#define m1 (equiv_9)
static integer k5;
#define m2 (equiv_1)
#define m3 (equiv_2)
#define m4 (equiv_3)
#define p0 (equiv_5)
static real p1;
#define p2 (equiv_7)
#define p3 (equiv_9)
#define q0 (equiv_6)
#define q1 (equiv_4)
#define q2 (equiv_11)
#define q3 (equiv_12)
#define t2 (equiv_4)
static real t3;
#define w2 (equiv_11)
#define w3 (equiv_12)
#define x2 (equiv_5)
static real x3;
static integer modem1;
static real x4, x5;
#define y2 (equiv_6)
static real y3, y4, y5;
extern doublereal gutre2_();
#define dz (equiv_8)
static real rm;
#define sw (equiv_10)
extern /* Subroutine */ int gd_message__();
static integer ierval[1], lm1, mm1;
static real cos2, cos3, sin2, sin3;
/* (Smooth Curve Fitting) */
/* This subroutine fits a smooth curve to a given set of input */
/* data points in an X-Y plane. It interpolates points in */
/* each interval between a pair of data points and generates a */
/* set of output points consisting of the input data points */
/* and the interpolated points. It can process either a */
/* single-valued function or a multiple-valued function. */
/* The input arguments are: */
/* MODE = mode of the curve (must be 1 or 2) */
/* = 1 for a single-valued function */
/* = 2 for multiple-valued function */
/* X = Array of dimension L storing the abscissas of input */
/* data points (in ascending or descending order for mode */
/* = 1) */
/* Y = Array of dimension L storing the ordinates of input */
/* data points */
/* L = Number of input data points (must be 2 or greater) */
/* M = Number of subintervals between each pair of input data */
/* points (must be 2 or greater). */
/* N = Number of output points */
/* = (L-1)*M+1 */
/* The output arguments are: */
/* U = Array of dimension N where the abscissas of output */
/* points are to be displayed */
/* V = Array of dimension N where the ordinates of output */
/* points are to be displayed */
/* Author: Hiroshi Akima, "Interpolation and Smooth Curve */
/* Fitting Based on Local Procedures", COMM. ACM 15, */
/* 914-918 (1972), and "A New Method of Interpolation */
/* and Smooth Curve Fitting Based on Local */
/* Procedures", J. ACM 17, 589-602 (1970). */
/* Corrections: M.R. Andersen, "Remark on Algorithm 433", ACM */
/* Trans. on Math. Software, 2, 208 (1976). */
/* (30-JAN-82) */
/* Preliminary processing */
/* Parameter adjustments */
--v;
--u;
--y;
--x;
/* Function Body */
mode0 = *mode;
modem1 = mode0 - 1;
l0 = *l;
lm1 = l0 - 1;
m0 = *m;
mm1 = m0 - 1;
n0 = *n;
if (mode0 <= 0) {
goto L320;
}
if (mode0 >= 3) {
goto L320;
}
if (lm1 <= 0) {
goto L330;
}
if (mm1 <= 0) {
goto L340;
}
if (n0 != lm1 * m0 + 1) {
gprint("Improper n value %ld, wanted %ld \n",n0,lm1 * m0 + 1);
gprint("n0 %ld, lm1 %ld, m0 %ld l0 %ld \n",n0,lm1,m0,l0);
goto L350;
}
switch (mode0) {
case 1: goto L10;
case 2: goto L60;
}
L10:
i = 2;
if ((r_1 = x[1] - x[2]) < (float)0.) {
goto L20;
} else if (r_1 == 0) {
gprint("Identical x values 1, %g %g \n",(double) x[1], (double) x[2]);
goto L360;
} else {
goto L40;
}
L20:
if (l0 <= 2) {
goto L80;
}
i_1 = l0;
for (i = 3; i <= i_1; ++i) {
if ((r_1 = x[i - 1] - x[i]) < (float)0.) {
goto L30;
} else if (r_1 == 0) {
gprint("Identical x values i, %g %g \n",(int) i,(double) x[i-1], (double) x[i]);
goto L360;
} else {
goto L370;
}
L30:
;}
goto L80;
L40:
if (l0 <= 2) {
goto L80;
}
i_1 = l0;
for (i = 3; i <= i_1; ++i) {
if ((r_1 = x[i - 1] - x[i]) < (float)0.) {
goto L370;
} else if (r_1 == 0) {
gprint("Identical x aavalues i, %g %g \n",(int) i,(double) x[i-1], (double) x[i]);
goto L360;
} else {
goto L50;
}
L50:
;}
goto L80;
L60:
i_1 = l0;
for (i = 2; i <= i_1; ++i) {
if (x[i - 1] != x[i]) {
goto L70;
}
if (y[i - 1] == y[i]) {
goto L380;
}
L70:
;}
L80:
k = n0 + m0;
i = l0 + 1;
i_1 = l0;
for (j = 1; j <= i_1; ++j) {
k -= m0;
--i;
u[k] = x[i];
/* L90: */
v[k] = y[i];
}
rm = (real) m0;
rm = (float)1. / rm;
/* Main DO-loop */
k5 = m0 + 1;
i_1 = l0;
for (i = 1; i <= i_1; ++i) {
/* Routines to pick up necessary X and Y values and to */
/* estimate them if necessary */
if (i > 1) {
goto L130;
}
x3 = u[1];
y3 = v[1];
x4 = u[m0 + 1];
y4 = v[m0 + 1];
a3 = x4 - x3;
*b3 = y4 - y3;
if (modem1 == 0) {
*m3 = *b3 / a3;
}
if (l0 != 2) {
goto L140;
}
a4 = a3;
*b4 = *b3;
L100:
switch (mode0) {
case 1: goto L120;
case 2: goto L110;
}
L110:
*a2 = a3 + a3 - a4;
*a1 = *a2 + *a2 - a3;
L120:
*b2 = *b3 + *b3 - *b4;
*b1 = *b2 + *b2 - *b3;
switch (mode0) {
case 1: goto L180;
case 2: goto L210;
}
L130:
*x2 = x3;
*y2 = y3;
x3 = x4;
y3 = y4;
x4 = x5;
y4 = y5;
*a1 = *a2;
*b1 = *b2;
*a2 = a3;
*b2 = *b3;
a3 = a4;
*b3 = *b4;
if (i >= lm1) {
goto L150;
}
L140:
k5 += m0;
x5 = u[k5];
y5 = v[k5];
a4 = x5 - x4;
*b4 = y5 - y4;
if (modem1 == 0) {
*m4 = *b4 / a4;
}
goto L160;
L150:
if (modem1 != 0) {
a4 = a3 + a3 - *a2;
}
*b4 = *b3 + *b3 - *b2;
L160:
if (i == 1) {
goto L100;
}
switch (mode0) {
case 1: goto L170;
case 2: goto L200;
}
/* Numerical differentiation */
L170:
*t2 = t3;
L180:
*w2 = (r_1 = *m4 - *m3, dabs(r_1));
*w3 = (r_1 = *m2 - *m1, dabs(r_1));
*sw = *w2 + *w3;
if (*sw != (float)0.) {
goto L190;
}
*w2 = (float).5;
*w3 = (float).5;
*sw = (float)1.;
L190:
t3 = (*w2 * *m2 + *w3 * *m3) / *sw;
if (i - 1 <= 0) {
goto L310;
} else {
goto L240;
}
L200:
cos2 = cos3;
sin2 = sin3;
L210:
*w2 = (r_1 = a3 * *b4 - a4 * *b3, dabs(r_1));
*w3 = (r_1 = *a1 * *b2 - *a2 * *b1, dabs(r_1));
if (*w2 + *w3 != (float)0.) {
goto L220;
}
*w2 = gutre2_(&a3, b3);
*w3 = gutre2_(a2, b2);
L220:
cos3 = *w2 * *a2 + *w3 * a3;
sin3 = *w2 * *b2 + *w3 * *b3;
*r = cos3 * cos3 + sin3 * sin3;
if (*r == (float)0.) {
goto L230;
}
*r = sqrt(*r);
cos3 /= *r;
sin3 /= *r;
L230:
if (i - 1 <= 0) {
goto L310;
} else {
goto L250;
}
/* Determination of the coefficients */
L240:
*q2 = ((*m2 - *t2) * (float)2. + *m2 - t3) / *a2;
*q3 = (-(doublereal)(*m2) - *m2 + *t2 + t3) / (*a2 * *a2);
goto L260;
L250:
*r = gutre2_(a2, b2);
p1 = *r * cos2;
*p2 = *a2 * (float)3. - *r * (cos2 + cos2 + cos3);
*p3 = *a2 - p1 - *p2;
*q1 = *r * sin2;
*q2 = *b2 * (float)3. - *r * (sin2 + sin2 + sin3);
*q3 = *b2 - *q1 - *q2;
goto L280;
/* Computation of the polynomials */
L260:
*dz = *a2 * rm;
*z = (float)0.;
i_2 = mm1;
for (j = 1; j <= i_2; ++j) {
++k;
*z += *dz;
u[k] = *p0 + *z;
/* L270: */
v[k] = *q0 + *z * (*q1 + *z * (*q2 + *z * *q3));
}
goto L300;
L280:
*z = (float)0.;
i_2 = mm1;
for (j = 1; j <= i_2; ++j) {
++k;
*z += rm;
u[k] = *p0 + *z * (p1 + *z * (*p2 + *z * *p3));
/* L290: */
v[k] = *q0 + *z * (*q1 + *z * (*q2 + *z * *q3));
}
L300:
++k;
L310:
;}
goto L410;
/* Error exit */
L320:
gd_message__("Cannot SMOOTH: Mode out of proper range 1..2", 29L);
goto L400;
L330:
gd_message__("Cannot SMOOTH: L = 1 or less", 13L);
goto L400;
L340:
gd_message__("Cannot SMOOTH: M = 1 or less", 13L);
goto L400;
L350:
gd_message__("Cannot SMOOTH: Improper N value", 16L);
goto L400;
L360:
gd_message__("Cannot SMOOTH: Identical X values", 18L);
goto L390;
L370:
gd_message__("Cannot SMOOTH: X values out of sequence", 24L);
goto L390;
L380:
gd_message__("Cannot SMOOTH: Identical X and Y values", 24L);
L390:
ierval[0] = i;
L400:
ierval[0] = mode0;
ierval[0] = l0;
ierval[0] = m0;
ierval[0] = n0;
L410:
return 0;
} /* glefitcf_ */
#undef sw
#undef dz
#undef y2
#undef x2
#undef w3
#undef w2
#undef t2
#undef q3
#undef q2
#undef q1
#undef q0
#undef p3
#undef p2
#undef p0
#undef m4
#undef m3
#undef m2
#undef m1
#undef b4
#undef b3
#undef b2
#undef b1
#undef a2
#undef a1
#undef z
#undef r
/* -- GUTRE2 */
doublereal gutre2_(a, b)
real *a, *b;
{
/* Initialized data */
static real zero = (float)0.;
static real two = (float)2.;
static real four = (float)4.;
/* System generated locals */
real ret_val, r_1;
/* Local variables */
static real aabs, babs, r, s;
/* (Euclidean Norm of 2-Vector) */
/* Return SQRT(A**2+B**2) without doing a square root, and */
/* without destructive underflow or overflow. */
/* (Cleve Moler, University of New Mexico) */
/* (09-APR-82) */
aabs = dabs(*a);
babs = dabs(*b);
/* --- IF (BABS .GT. AABS) THEN -- SWAP A AND B */
if (babs <= aabs) {
goto L10;
}
r = babs;
babs = aabs;
aabs = r;
/* --- END IF */
/* --- IF (BABS .EQ. ZERO) THEN -- Special case sudden exit */
L10:
if (babs != zero) {
goto L20;
}
ret_val = aabs;
goto L40;
/* --- END IF */
/* --- DO FOREVER */
L20:
/* Computing 2nd power */
r_1 = babs / aabs;
r = r_1 * r_1;
/* --- IF ((TWO+R) .EQ. TWO) THEN -- Converged */
if (two + r != two) {
goto L30;
}
ret_val = aabs;
goto L40;
/* --- END IF */
L30:
s = r / (four + r);
aabs += two * s * aabs;
babs = s * babs;
/* --- END FOREVER */
goto L20;
L40:
return ret_val;
} /* gutre2_ */
gd_message__(char *s,long l)
{
gprint("%s %ld \n",s,l);
}