PLplot
5.10.0
|
00001 // 00002 // Copyright (C) 2009 Alan W. Irwin 00003 // 00004 // This file is part of PLplot. 00005 // 00006 // PLplot is free software; you can redistribute it and/or modify 00007 // it under the terms of the GNU Library General Public License as published 00008 // by the Free Software Foundation; either version 2 of the License, or 00009 // (at your option) any later version. 00010 // 00011 // PLplot is distributed in the hope that it will be useful, 00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 // GNU Library General Public License for more details. 00015 // 00016 // You should have received a copy of the GNU Library General Public License 00017 // along with PLplot; if not, write to the Free Software 00018 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00019 // 00020 // Provenance: This code was originally developed under the GPL as part of 00021 // the FreeEOS project (revision 121). This code has been converted from 00022 // Fortran to C with the aid of f2c and relicensed for PLplot under the LGPL 00023 // with the permission of the FreeEOS copyright holder (Alan W. Irwin). 00024 // 00025 00026 #include "dspline.h" 00027 00028 int dspline( double *x, double *y, int n, 00029 int if1, double cond1, int ifn, double condn, double *y2 ) 00030 { 00031 int i__1, i__, k; 00032 double p, u[2000], qn, un, sig; 00033 00034 // input parameters: 00035 // x(n) are the spline knot points 00036 // y(n) are the function values at the knot points 00037 // if1 = 1 specifies cond1 is the first derivative at the 00038 // first knot point. 00039 // if1 = 2 specifies cond1 is the second derivative at the 00040 // first knot point. 00041 // ifn = 1 specifies condn is the first derivative at the 00042 // nth knot point. 00043 // ifn = 2 specifies condn is the second derivative at the 00044 // nth knot point. 00045 // output values: 00046 // y2(n) is the second derivative of the spline evaluated at 00047 // the knot points. 00048 // Parameter adjustments 00049 --y2; 00050 --y; 00051 --x; 00052 00053 // Function Body 00054 if ( n > 2000 ) 00055 { 00056 return 1; 00057 } 00058 // y2(i) = u(i) + d(i)*y2(i+1), where 00059 // d(i) is temporarily stored in y2(i) (see below). 00060 if ( if1 == 2 ) 00061 { 00062 // cond1 is second derivative at first point. 00063 // these two values assure that for above equation with d(i) temporarily 00064 // stored in y2(i) 00065 y2[1] = 0.; 00066 u[0] = cond1; 00067 } 00068 else if ( if1 == 1 ) 00069 { 00070 // cond1 is first derivative at first point. 00071 // special case (Press et al 3.3.5 with A = 1, and B=0) 00072 // of equations below where 00073 // a_j = 0 00074 // b_j = -(x_j+1 - x_j)/3 00075 // c_j = -(x_j+1 - x_j)/6 00076 // r_j = cond1 - (y_j+1 - y_j)/(x_j+1 - x_j) 00077 // u(i) = r(i)/b(i) 00078 // d(i) = -c(i)/b(i) 00079 // N.B. d(i) is temporarily stored in y2. 00080 y2[1] = -.5; 00081 u[0] = 3. / ( x[2] - x[1] ) * ( ( y[2] - y[1] ) / ( x[2] - x[1] ) - cond1 ); 00082 } 00083 else 00084 { 00085 return 2; 00086 } 00087 // if original tri-diagonal system is characterized as 00088 // a_j y2_j-1 + b_j y2_j + c_j y2_j+1 = r_j 00089 // Then from Press et al. 3.3.7, we have the unscaled result: 00090 // a_j = (x_j - x_j-1)/6 00091 // b_j = (x_j+1 - x_j-1)/3 00092 // c_j = (x_j+1 - x_j)/6 00093 // r_j = (y_j+1 - y_j)/(x_j+1 - x_j) - (y_j - y_j-1)/(x_j - x_j-1) 00094 // In practice, all these values are divided through by b_j/2 to scale 00095 // them, and from now on we will use these scaled values. 00096 00097 // forward elimination step: assume y2(i-1) = u(i-1) + d(i-1)*y2(i). 00098 // When this is substituted into above tridiagonal equation ==> 00099 // y2(i) = u(i) + d(i)*y2(i+1), where 00100 // u(i) = [r(i) - a(i) u(i-1)]/[b(i) + a(i) d(i-1)] 00101 // d(i) = -c(i)/[b(i) + a(i) d(i-1)] 00102 // N.B. d(i) is temporarily stored in y2. 00103 i__1 = n - 1; 00104 for ( i__ = 2; i__ <= i__1; ++i__ ) 00105 { 00106 // sig is scaled a(i) 00107 sig = ( x[i__] - x[i__ - 1] ) / ( x[i__ + 1] - x[i__ - 1] ); 00108 // p is denominator = scaled a(i) d(i-1) + scaled b(i), where scaled 00109 // b(i) is 2. 00110 p = sig * y2[i__ - 1] + 2.; 00111 // propagate d(i) equation above. Note sig-1 = -c(i) 00112 y2[i__] = ( sig - 1. ) / p; 00113 // propagate scaled u(i) equation above 00114 u[i__ - 1] = ( ( ( y[i__ + 1] - y[i__] ) / ( x[i__ + 1] - x[i__] ) - ( y[i__] 00115 - y[i__ - 1] ) / ( x[i__] - x[i__ - 1] ) ) * 6. / ( x[i__ + 1] - 00116 x[i__ - 1] ) - sig * u[i__ - 2] ) / p; 00117 } 00118 if ( ifn == 2 ) 00119 { 00120 // condn is second derivative at nth point. 00121 // These two values assure that in the equation below. 00122 qn = 0.; 00123 un = condn; 00124 } 00125 else if ( ifn == 1 ) 00126 { 00127 // specify condn is first derivative at nth point. 00128 // special case (Press et al 3.3.5 with A = 0, and B=1) 00129 // implies a_n y2(n-1) + b_n y2(n) = r_n, where 00130 // a_n = (x_n - x_n-1)/6 00131 // b_n = (x_n - x_n-1)/3 00132 // r_n = cond1 - (y_n - y_n-1)/(x_n - x_n-1) 00133 // use same propagation equation as above, only with c_n = 0 00134 // ==> d_n = 0 ==> y2(n) = u(n) => 00135 // y(n) = [r(n) - a(n) u(n-1)]/[b(n) + a(n) d(n-1)] 00136 // qn is scaled a_n 00137 qn = .5; 00138 // un is scaled r_n (N.B. un is not u(n))! Sorry for the mixed notation. 00139 un = 3. / ( x[n] - x[n - 1] ) * ( condn - ( y[n] - y[n - 1] ) / ( x[n] 00140 - x[n - 1] ) ); 00141 } 00142 else 00143 { 00144 return 3; 00145 } 00146 // N.B. d(i) is temporarily stored in y2, and everything is 00147 // scaled by b_n. 00148 // qn is scaled a_n, 1.d0 is scaled b_n, and un is scaled r_n. 00149 y2[n] = ( un - qn * u[n - 2] ) / ( qn * y2[n - 1] + 1. ); 00150 // back substitution. 00151 for ( k = n - 1; k >= 1; --k ) 00152 { 00153 y2[k] = y2[k] * y2[k + 1] + u[k - 1]; 00154 } 00155 return 0; 00156 } 00157