Revision: 10037
http://plplot.svn.sourceforge.net/plplot/?rev=10037&view=rev
Author: airwin
Date: 2009-06-06 08:17:29 +0000 (Sat, 06 Jun 2009)
Log Message:
-----------
Fix up a number of issues revealed by deltaT-gen.c and deltaT_test.c. With
these changes, the C and Fortran versions of all deltaT-related code gives
identical results proving the conversion from Fortran to C has been a success.
Modified Paths:
--------------
trunk/lib/qsastime/dspline.c
trunk/lib/qsastime/dsplint.c
Modified: trunk/lib/qsastime/dspline.c
===================================================================
--- trunk/lib/qsastime/dspline.c 2009-06-06 08:12:03 UTC (rev 10036)
+++ trunk/lib/qsastime/dspline.c 2009-06-06 08:17:29 UTC (rev 10037)
@@ -22,20 +22,12 @@
Fortran to C with the aid of f2c and relicensed for PLplot under the LGPL
with the permission of the FreeEOS copyright holder (Alan W. Irwin).
*/
-int dspline_(doublereal *x, doublereal *y, integer *n,
- integer *if1, doublereal *cond1, integer *ifn, doublereal *condn,
- doublereal *y2)
+int dspline(double *x, double *y, int n,
+ int if1, double cond1, int ifn, double condn, double *y2)
{
- /* System generated locals */
- integer i__1;
+ int i__1, i__, k;
+ double p, u[2000], qn, un, sig;
- /* Builtin functions */
- /* Subroutine */ int s_stop(char *, ftnlen);
-
- /* Local variables */
- static integer i__, k;
- static doublereal p, u[2000], qn, un, sig;
-
/* input parameters: */
/* x(n) are the spline knot points */
/* y(n) are the function values at the knot points */
@@ -56,18 +48,18 @@
--x;
/* Function Body */
- if (*n > 2000) {
- s_stop("dspline: dimensions too large", (ftnlen)29);
+ if (n > 2000) {
+ return 1;
}
/* y2(i) = u(i) + d(i)*y2(i+1), where */
/* d(i) is temporarily stored in y2(i) (see below). */
- if (*if1 == 2) {
+ if (if1 == 2) {
/* cond1 is second derivative at first point. */
/* these two values assure that for above equation with d(i) temporarily */
/* stored in y2(i) */
y2[1] = 0.;
- u[0] = *cond1;
- } else if (*if1 == 1) {
+ u[0] = cond1;
+ } else if (if1 == 1) {
/* cond1 is first derivative at first point. */
/* special case (Press et al 3.3.5 with A = 1, and B=0) */
/* of equations below where */
@@ -79,9 +71,9 @@
/* d(i) = -c(i)/b(i) */
/* N.B. d(i) is temporarily stored in y2. */
y2[1] = -.5;
- u[0] = 3. / (x[2] - x[1]) * ((y[2] - y[1]) / (x[2] - x[1]) - *cond1);
+ u[0] = 3. / (x[2] - x[1]) * ((y[2] - y[1]) / (x[2] - x[1]) - cond1);
} else {
- s_stop("dspline: incorrect if1 value specified", (ftnlen)38);
+ return 2;
}
/* if original tri-diagonal system is characterized as */
/* a_j y2_j-1 + b_j y2_j + c_j y2_j+1 = r_j */
@@ -99,7 +91,7 @@
/* u(i) = [r(i) - a(i) u(i-1)]/[b(i) + a(i) d(i-1)] */
/* d(i) = -c(i)/[b(i) + a(i) d(i-1)] */
/* N.B. d(i) is temporarily stored in y2. */
- i__1 = *n - 1;
+ i__1 = n - 1;
for (i__ = 2; i__ <= i__1; ++i__) {
/* sig is scaled a(i) */
sig = (x[i__] - x[i__ - 1]) / (x[i__ + 1] - x[i__ - 1]);
@@ -113,12 +105,12 @@
- y[i__ - 1]) / (x[i__] - x[i__ - 1])) * 6. / (x[i__ + 1] -
x[i__ - 1]) - sig * u[i__ - 2]) / p;
}
- if (*ifn == 2) {
+ if (ifn == 2) {
/* condn is second derivative at nth point. */
/* These two values assure that in the equation below. */
qn = 0.;
- un = *condn;
- } else if (*ifn == 1) {
+ un = condn;
+ } else if (ifn == 1) {
/* specify condn is first derivative at nth point. */
/* special case (Press et al 3.3.5 with A = 0, and B=1) */
/* implies a_n y2(n-1) + b_n y2(n) = r_n, where */
@@ -131,19 +123,19 @@
/* qn is scaled a_n */
qn = .5;
/* un is scaled r_n (N.B. un is not u(n))! Sorry for the mixed notation. */
- un = 3. / (x[*n] - x[*n - 1]) * (*condn - (y[*n] - y[*n - 1]) / (x[*n]
- - x[*n - 1]));
+ un = 3. / (x[n] - x[n - 1]) * (condn - (y[n] - y[n - 1]) / (x[n]
+ - x[n - 1]));
} else {
- s_stop("dspline: incorrect ifn value specified", (ftnlen)38);
+ return 3;
}
/* N.B. d(i) is temporarily stored in y2, and everything is */
/* scaled by b_n. */
/* qn is scaled a_n, 1.d0 is scaled b_n, and un is scaled r_n. */
- y2[*n] = (un - qn * u[*n - 2]) / (qn * y2[*n - 1] + 1.);
+ y2[n] = (un - qn * u[n - 2]) / (qn * y2[n - 1] + 1.);
/* back substitution. */
- for (k = *n - 1; k >= 1; --k) {
+ for (k = n - 1; k >= 1; --k) {
y2[k] = y2[k] * y2[k + 1] + u[k - 1];
}
return 0;
-} /* dspline_ */
+}
Modified: trunk/lib/qsastime/dsplint.c
===================================================================
--- trunk/lib/qsastime/dsplint.c 2009-06-06 08:12:03 UTC (rev 10036)
+++ trunk/lib/qsastime/dsplint.c 2009-06-06 08:17:29 UTC (rev 10037)
@@ -22,25 +22,19 @@
Fortran to C with the aid of f2c and relicensed for PLplot under the LGPL
with the permission of the FreeEOS copyright holder (Alan W. Irwin).
*/
+# define MAX(a,b) (((a) > (b)) ? (a) : (b))
+# define MIN(a,b) (((a) < (b)) ? (a) : (b))
-int dsplint_(doublereal *xa, doublereal *ya, doublereal *y2a,
- integer *n, doublereal *x, doublereal *y, doublereal *dy, doublereal
- *d2y)
+int dsplint(double *xa, double *ya, double *y2a,
+ int n, double x, double *y, double *dy, double *d2y)
{
/* Initialized data */
- static integer nsave = 0;
+ static int nsave = 0, khi, klo;
- /* System generated locals */
- integer i__1, i__2;
+ int i__1, i__2, k;
+ double a, b, h__;
- /* Builtin functions */
- /* Subroutine */ int s_stop(char *, ftnlen);
-
- /* Local variables */
- static doublereal a, b, h__;
- static integer k, khi, klo;
-
/* evaluate spline = y and its derivatives dy and d2y at x given */
/* xa, ya, y2a from dspline. */
/* Parameter adjustments */
@@ -49,49 +43,49 @@
--xa;
/* Function Body */
- if (*n != nsave) {
+ if (n != nsave) {
/* if call with different n value, then redo range */
- nsave = *n;
+ nsave = n;
klo = 1;
- khi = *n;
- if (xa[klo] > *x) {
- s_stop("dsplint: x too low", (ftnlen)18);
+ khi = n;
+ if (xa[klo] > x) {
+ return 1;
}
- if (xa[khi] < *x) {
- s_stop("dsplint: x too high", (ftnlen)19);
+ if (xa[khi] < x) {
+ return 2;
}
} else {
/* optimize range assuming continuous (ascending or */
/* descending x calls. */
- if (xa[klo] > *x) {
+ if (xa[klo] > x) {
/* x is descending so try next range. */
- khi = max(2,klo);
+ khi = MAX(2,klo);
klo = khi - 1;
/* if x smaller than next range try lower limit. */
- if (xa[klo] > *x) {
+ if (xa[klo] > x) {
klo = 1;
}
- if (xa[klo] > *x) {
- s_stop("dsplint: x too low", (ftnlen)18);
+ if (xa[klo] > x) {
+ return 1;
}
- } else if (xa[khi] <= *x) {
+ } else if (xa[khi] <= x) {
/* x is ascending so try next range. */
/* Computing MIN */
- i__1 = khi, i__2 = *n - 1;
- klo = min(i__1,i__2);
+ i__1 = khi, i__2 = n - 1;
+ klo = MIN(i__1,i__2);
khi = klo + 1;
/* if x larger than next range try upper limit. */
- if (xa[khi] <= *x) {
- khi = *n;
+ if (xa[khi] <= x) {
+ khi = n;
}
- if (xa[khi] < *x) {
- s_stop("dsplint: x too high", (ftnlen)19);
+ if (xa[khi] < x) {
+ return 2;
}
}
}
while(khi - klo > 1) {
k = (khi + klo) / 2;
- if (xa[k] > *x) {
+ if (xa[k] > x) {
khi = k;
} else {
klo = k;
@@ -99,15 +93,15 @@
}
h__ = xa[khi] - xa[klo];
if (h__ <= 0.) {
- s_stop("dsplint: bad xa input.", (ftnlen)22);
+ return 3;
}
- a = (xa[khi] - *x) / h__;
- b = (*x - xa[klo]) / h__;
+ a = (xa[khi] - x) / h__;
+ b = (x - xa[klo]) / h__;
*y = a * ya[klo] + b * ya[khi] + (a * (a * a - 1.) * y2a[klo] + b * (b *
b - 1.) * y2a[khi]) * (h__ * h__) / 6.;
*dy = (-ya[klo] + ya[khi] + (-(a * 3. * a - 1.) * y2a[klo] + (b * 3. * b
- 1.) * y2a[khi]) * (h__ * h__) / 6.) / h__;
*d2y = a * y2a[klo] + b * y2a[khi];
return 0;
-} /* dsplint_ */
+}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|