1 //
2 // Copyright (C) 2009 Alan W. Irwin
3 //
4 // This file is part of PLplot.
5 //
6 // PLplot is free software; you can redistribute it and/or modify
7 // it under the terms of the GNU Library General Public License as published
8 // by the Free Software Foundation; either version 2 of the License, or
9 // (at your option) any later version.
10 //
11 // PLplot is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU Library General Public License for more details.
15 //
16 // You should have received a copy of the GNU Library General Public License
17 // along with PLplot; if not, write to the Free Software
18 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
19 //
20 // Provenance: This code was originally developed under the GPL as part of
21 // the FreeEOS project (revision 121). This code has been converted from
22 // Fortran to C with the aid of f2c and relicensed for PLplot under the LGPL
23 // with the permission of the FreeEOS copyright holder (Alan W. Irwin).
24 //
25
26 #include "dspline.h"
27
dspline(double * x,double * y,int n,int if1,double cond1,int ifn,double condn,double * y2)28 int dspline( double *x, double *y, int n,
29 int if1, double cond1, int ifn, double condn, double *y2 )
30 {
31 int i__1, i__, k;
32 double p, u[2000], qn, un, sig;
33
34 // input parameters:
35 // x(n) are the spline knot points
36 // y(n) are the function values at the knot points
37 // if1 = 1 specifies cond1 is the first derivative at the
38 // first knot point.
39 // if1 = 2 specifies cond1 is the second derivative at the
40 // first knot point.
41 // ifn = 1 specifies condn is the first derivative at the
42 // nth knot point.
43 // ifn = 2 specifies condn is the second derivative at the
44 // nth knot point.
45 // output values:
46 // y2(n) is the second derivative of the spline evaluated at
47 // the knot points.
48 // Parameter adjustments
49 --y2;
50 --y;
51 --x;
52
53 // Function Body
54 if ( n > 2000 )
55 {
56 return 1;
57 }
58 // y2(i) = u(i) + d(i)*y2(i+1), where
59 // d(i) is temporarily stored in y2(i) (see below).
60 if ( if1 == 2 )
61 {
62 // cond1 is second derivative at first point.
63 // these two values assure that for above equation with d(i) temporarily
64 // stored in y2(i)
65 y2[1] = 0.;
66 u[0] = cond1;
67 }
68 else if ( if1 == 1 )
69 {
70 // cond1 is first derivative at first point.
71 // special case (Press et al 3.3.5 with A = 1, and B=0)
72 // of equations below where
73 // a_j = 0
74 // b_j = -(x_j+1 - x_j)/3
75 // c_j = -(x_j+1 - x_j)/6
76 // r_j = cond1 - (y_j+1 - y_j)/(x_j+1 - x_j)
77 // u(i) = r(i)/b(i)
78 // d(i) = -c(i)/b(i)
79 // N.B. d(i) is temporarily stored in y2.
80 y2[1] = -.5;
81 u[0] = 3. / ( x[2] - x[1] ) * ( ( y[2] - y[1] ) / ( x[2] - x[1] ) - cond1 );
82 }
83 else
84 {
85 return 2;
86 }
87 // if original tri-diagonal system is characterized as
88 // a_j y2_j-1 + b_j y2_j + c_j y2_j+1 = r_j
89 // Then from Press et al. 3.3.7, we have the unscaled result:
90 // a_j = (x_j - x_j-1)/6
91 // b_j = (x_j+1 - x_j-1)/3
92 // c_j = (x_j+1 - x_j)/6
93 // r_j = (y_j+1 - y_j)/(x_j+1 - x_j) - (y_j - y_j-1)/(x_j - x_j-1)
94 // In practice, all these values are divided through by b_j/2 to scale
95 // them, and from now on we will use these scaled values.
96
97 // forward elimination step: assume y2(i-1) = u(i-1) + d(i-1)*y2(i).
98 // When this is substituted into above tridiagonal equation ==>
99 // y2(i) = u(i) + d(i)*y2(i+1), where
100 // u(i) = [r(i) - a(i) u(i-1)]/[b(i) + a(i) d(i-1)]
101 // d(i) = -c(i)/[b(i) + a(i) d(i-1)]
102 // N.B. d(i) is temporarily stored in y2.
103 i__1 = n - 1;
104 for ( i__ = 2; i__ <= i__1; ++i__ )
105 {
106 // sig is scaled a(i)
107 sig = ( x[i__] - x[i__ - 1] ) / ( x[i__ + 1] - x[i__ - 1] );
108 // p is denominator = scaled a(i) d(i-1) + scaled b(i), where scaled
109 // b(i) is 2.
110 p = sig * y2[i__ - 1] + 2.;
111 // propagate d(i) equation above. Note sig-1 = -c(i)
112 y2[i__] = ( sig - 1. ) / p;
113 // propagate scaled u(i) equation above
114 u[i__ - 1] = ( ( ( y[i__ + 1] - y[i__] ) / ( x[i__ + 1] - x[i__] ) - ( y[i__]
115 - y[i__ - 1] ) / ( x[i__] - x[i__ - 1] ) ) * 6. / ( x[i__ + 1] -
116 x[i__ - 1] ) - sig * u[i__ - 2] ) / p;
117 }
118 if ( ifn == 2 )
119 {
120 // condn is second derivative at nth point.
121 // These two values assure that in the equation below.
122 qn = 0.;
123 un = condn;
124 }
125 else if ( ifn == 1 )
126 {
127 // specify condn is first derivative at nth point.
128 // special case (Press et al 3.3.5 with A = 0, and B=1)
129 // implies a_n y2(n-1) + b_n y2(n) = r_n, where
130 // a_n = (x_n - x_n-1)/6
131 // b_n = (x_n - x_n-1)/3
132 // r_n = cond1 - (y_n - y_n-1)/(x_n - x_n-1)
133 // use same propagation equation as above, only with c_n = 0
134 // ==> d_n = 0 ==> y2(n) = u(n) =>
135 // y(n) = [r(n) - a(n) u(n-1)]/[b(n) + a(n) d(n-1)]
136 // qn is scaled a_n
137 qn = .5;
138 // un is scaled r_n (N.B. un is not u(n))! Sorry for the mixed notation.
139 un = 3. / ( x[n] - x[n - 1] ) * ( condn - ( y[n] - y[n - 1] ) / ( x[n]
140 - x[n - 1] ) );
141 }
142 else
143 {
144 return 3;
145 }
146 // N.B. d(i) is temporarily stored in y2, and everything is
147 // scaled by b_n.
148 // qn is scaled a_n, 1.d0 is scaled b_n, and un is scaled r_n.
149 y2[n] = ( un - qn * u[n - 2] ) / ( qn * y2[n - 1] + 1. );
150 // back substitution.
151 for ( k = n - 1; k >= 1; --k )
152 {
153 y2[k] = y2[k] * y2[k + 1] + u[k - 1];
154 }
155 return 0;
156 }
157
158