00001 /* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 #include "cddefines.h" 00004 #include "thirdparty.h" 00005 00006 STATIC double *d3_np_fs( long n, double a[], const double b[], double x[] ); 00007 00008 //****************************************************************************80 00009 // 00010 // The routines d3_np_fs, spline_cubic_set, spline_cubic_val where written 00011 // by John Burkardt (Computer Science Department, Florida State University) 00012 // and have been slightly modified and adapted for use in Cloudy by Peter 00013 // van Hoof (Royal Observatory of Belgium). 00014 // 00015 // The original sources can be found at 00016 // http://www.scs.fsu.edu/~burkardt/cpp_src/spline/spline.html 00017 // 00018 //****************************************************************************80 00019 00020 STATIC double *d3_np_fs( long n, double a[], const double b[], double x[] ) 00021 00022 //****************************************************************************80 00023 // 00024 // Purpose: 00025 // 00026 // D3_NP_FS factors and solves a D3 system. 00027 // 00028 // Discussion: 00029 // 00030 // The D3 storage format is used for a tridiagonal matrix. 00031 // The superdiagonal is stored in entries (1,2:N), the diagonal in 00032 // entries (2,1:N), and the subdiagonal in (3,1:N-1). Thus, the 00033 // original matrix is "collapsed" vertically into the array. 00034 // 00035 // This algorithm requires that each diagonal entry be nonzero. 00036 // It does not use pivoting, and so can fail on systems that 00037 // are actually nonsingular. 00038 // 00039 // Example: 00040 // 00041 // Here is how a D3 matrix of order 5 would be stored: 00042 // 00043 // * A12 A23 A34 A45 00044 // A11 A22 A33 A44 A55 00045 // A21 A32 A43 A54 * 00046 // 00047 // Modified: 00048 // 00049 // 15 November 2003 00050 // 00051 // Author: 00052 // 00053 // John Burkardt 00054 // 00055 // Parameters: 00056 // 00057 // Input, long N, the order of the linear system. 00058 // 00059 // Input/output, double A[3*N]. 00060 // On input, the nonzero diagonals of the linear system. 00061 // On output, the data in these vectors has been overwritten 00062 // by factorization information. 00063 // 00064 // Input, double B[N], the right hand side. 00065 // 00066 // Output, double X[N] and D3_NP_FS[N], the solution of the linear system. 00067 // D3_NP_FS returns NULL if there was an error because one of the diagonal 00068 // entries was zero. 00069 // 00070 { 00071 long i; 00072 double xmult; 00073 00074 DEBUG_ENTRY( "d3_np_fs()" ); 00075 // 00076 // Check. 00077 // 00078 for( i = 0; i < n; i++ ) 00079 { 00080 if( a[1+i*3] == 0.0 ) 00081 { 00082 return NULL; 00083 } 00084 } 00085 00086 x[0] = b[0]; 00087 00088 for( i = 1; i < n; i++ ) 00089 { 00090 xmult = a[2+(i-1)*3] / a[1+(i-1)*3]; 00091 a[1+i*3] = a[1+i*3] - xmult * a[0+i*3]; 00092 x[i] = b[i] - xmult * x[i-1]; 00093 } 00094 00095 x[n-1] = x[n-1] / a[1+(n-1)*3]; 00096 for( i = n-2; 0 <= i; i-- ) 00097 { 00098 x[i] = ( x[i] - a[0+(i+1)*3] * x[i+1] ) / a[1+i*3]; 00099 } 00100 return x; 00101 } 00102 00103 //****************************************************************************80 00104 00105 void spline_cubic_set( long n, const double t[], const double y[], double ypp[], 00106 int ibcbeg, double ybcbeg, int ibcend, double ybcend ) 00107 00108 //****************************************************************************80 00109 // 00110 // Purpose: 00111 // 00112 // SPLINE_CUBIC_SET computes the second derivatives of a piecewise cubic spline. 00113 // 00114 // Discussion: 00115 // 00116 // For data interpolation, the user must call SPLINE_SET to determine 00117 // the second derivative data, passing in the data to be interpolated, 00118 // and the desired boundary conditions. 00119 // 00120 // The data to be interpolated, plus the SPLINE_SET output, defines 00121 // the spline. The user may then call SPLINE_VAL to evaluate the 00122 // spline at any point. 00123 // 00124 // The cubic spline is a piecewise cubic polynomial. The intervals 00125 // are determined by the "knots" or abscissas of the data to be 00126 // interpolated. The cubic spline has continous first and second 00127 // derivatives over the entire interval of interpolation. 00128 // 00129 // For any point T in the interval T(IVAL), T(IVAL+1), the form of 00130 // the spline is 00131 // 00132 // SPL(T) = A(IVAL) 00133 // + B(IVAL) * ( T - T(IVAL) ) 00134 // + C(IVAL) * ( T - T(IVAL) )**2 00135 // + D(IVAL) * ( T - T(IVAL) )**3 00136 // 00137 // If we assume that we know the values Y(*) and YPP(*), which represent 00138 // the values and second derivatives of the spline at each knot, then 00139 // the coefficients can be computed as: 00140 // 00141 // A(IVAL) = Y(IVAL) 00142 // B(IVAL) = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) ) 00143 // - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6 00144 // C(IVAL) = YPP(IVAL) / 2 00145 // D(IVAL) = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) ) 00146 // 00147 // Since the first derivative of the spline is 00148 // 00149 // SPL'(T) = B(IVAL) 00150 // + 2 * C(IVAL) * ( T - T(IVAL) ) 00151 // + 3 * D(IVAL) * ( T - T(IVAL) )**2, 00152 // 00153 // the requirement that the first derivative be continuous at interior 00154 // knot I results in a total of N-2 equations, of the form: 00155 // 00156 // B(IVAL-1) + 2 C(IVAL-1) * (T(IVAL)-T(IVAL-1)) 00157 // + 3 * D(IVAL-1) * (T(IVAL) - T(IVAL-1))**2 = B(IVAL) 00158 // 00159 // or, setting H(IVAL) = T(IVAL+1) - T(IVAL) 00160 // 00161 // ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1) 00162 // - ( YPP(IVAL) + 2 * YPP(IVAL-1) ) * H(IVAL-1) / 6 00163 // + YPP(IVAL-1) * H(IVAL-1) 00164 // + ( YPP(IVAL) - YPP(IVAL-1) ) * H(IVAL-1) / 2 00165 // = 00166 // ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL) 00167 // - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * H(IVAL) / 6 00168 // 00169 // or 00170 // 00171 // YPP(IVAL-1) * H(IVAL-1) + 2 * YPP(IVAL) * ( H(IVAL-1) + H(IVAL) ) 00172 // + YPP(IVAL) * H(IVAL) 00173 // = 00174 // 6 * ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL) 00175 // - 6 * ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1) 00176 // 00177 // Boundary conditions must be applied at the first and last knots. 00178 // The resulting tridiagonal system can be solved for the YPP values. 00179 // 00180 // Modified: 00181 // 00182 // 06 February 2004 00183 // 00184 // Author: 00185 // 00186 // John Burkardt 00187 // 00188 // Parameters: 00189 // 00190 // Input, long N, the number of data points. N must be at least 2. 00191 // In the special case where N = 2 and IBCBEG = IBCEND = 0, the 00192 // spline will actually be linear. 00193 // 00194 // Input, double T[N], the knot values, that is, the points were data is 00195 // specified. The knot values should be distinct, and increasing. 00196 // 00197 // Input, double Y[N], the data values to be interpolated. 00198 // 00199 // Input, int IBCBEG, left boundary condition flag: 00200 // 0: the cubic spline should be a quadratic over the first interval; 00201 // 1: the first derivative at the left endpoint should be YBCBEG; 00202 // 2: the second derivative at the left endpoint should be YBCBEG. 00203 // 00204 // Input, double YBCBEG, the values to be used in the boundary 00205 // conditions if IBCBEG is equal to 1 or 2. 00206 // 00207 // Input, int IBCEND, right boundary condition flag: 00208 // 0: the cubic spline should be a quadratic over the last interval; 00209 // 1: the first derivative at the right endpoint should be YBCEND; 00210 // 2: the second derivative at the right endpoint should be YBCEND. 00211 // 00212 // Input, double YBCEND, the values to be used in the boundary 00213 // conditions if IBCEND is equal to 1 or 2. 00214 // 00215 // Output, double YPP[N] and SPLINE_CUBIC_SET[N], the second derivatives 00216 // of the cubic spline. 00217 // 00218 { 00219 long i; 00220 00221 DEBUG_ENTRY( "spline_cubic_set()" ); 00222 // 00223 // Check. 00224 // 00225 ASSERT( n >= 2 ); 00226 00227 # ifndef NDEBUG 00228 for( i = 0; i < n - 1; i++ ) 00229 { 00230 if( t[i+1] <= t[i] ) 00231 { 00232 fprintf( ioQQQ, "SPLINE_CUBIC_SET - Fatal error!\n" ); 00233 fprintf( ioQQQ, " The knots must be strictly increasing\n" ); 00234 cdEXIT(EXIT_FAILURE); 00235 } 00236 } 00237 # endif 00238 00239 double *a = (double*)MALLOC((size_t)(3*n)*sizeof(double)); 00240 double *b = (double*)MALLOC((size_t)n*sizeof(double)); 00241 // 00242 // Set up the first equation. 00243 // 00244 if( ibcbeg == 0 ) 00245 { 00246 b[0] = 0.0; 00247 a[1+0*3] = 1.0; 00248 a[0+1*3] = -1.0; 00249 } 00250 else if( ibcbeg == 1 ) 00251 { 00252 b[0] = ( y[1] - y[0] ) / ( t[1] - t[0] ) - ybcbeg; 00253 a[1+0*3] = ( t[1] - t[0] ) / 3.0; 00254 a[0+1*3] = ( t[1] - t[0] ) / 6.0; 00255 } 00256 else if( ibcbeg == 2 ) 00257 { 00258 b[0] = ybcbeg; 00259 a[1+0*3] = 1.0; 00260 a[0+1*3] = 0.0; 00261 } 00262 else 00263 { 00264 fprintf( ioQQQ, "SPLINE_CUBIC_SET - Fatal error!\n" ); 00265 fprintf( ioQQQ, " IBCBEG must be 0, 1 or 2, but I found %d.\n", ibcbeg ); 00266 cdEXIT(EXIT_FAILURE); 00267 } 00268 // 00269 // Set up the intermediate equations. 00270 // 00271 for( i = 1; i < n-1; i++ ) 00272 { 00273 b[i] = ( y[i+1] - y[i] ) / ( t[i+1] - t[i] ) 00274 - ( y[i] - y[i-1] ) / ( t[i] - t[i-1] ); 00275 a[2+(i-1)*3] = ( t[i] - t[i-1] ) / 6.0; 00276 a[1+ i *3] = ( t[i+1] - t[i-1] ) / 3.0; 00277 a[0+(i+1)*3] = ( t[i+1] - t[i] ) / 6.0; 00278 } 00279 // 00280 // Set up the last equation. 00281 // 00282 if( ibcend == 0 ) 00283 { 00284 b[n-1] = 0.0; 00285 a[2+(n-2)*3] = -1.0; 00286 a[1+(n-1)*3] = 1.0; 00287 } 00288 else if( ibcend == 1 ) 00289 { 00290 b[n-1] = ybcend - ( y[n-1] - y[n-2] ) / ( t[n-1] - t[n-2] ); 00291 a[2+(n-2)*3] = ( t[n-1] - t[n-2] ) / 6.0; 00292 a[1+(n-1)*3] = ( t[n-1] - t[n-2] ) / 3.0; 00293 } 00294 else if( ibcend == 2 ) 00295 { 00296 b[n-1] = ybcend; 00297 a[2+(n-2)*3] = 0.0; 00298 a[1+(n-1)*3] = 1.0; 00299 } 00300 else 00301 { 00302 fprintf( ioQQQ, "SPLINE_CUBIC_SET - Fatal error!\n" ); 00303 fprintf( ioQQQ, " IBCEND must be 0, 1 or 2, but I found %d.\n", ibcend ); 00304 cdEXIT(EXIT_FAILURE); 00305 } 00306 // 00307 // Solve the linear system. 00308 // 00309 if( n == 2 && ibcbeg == 0 && ibcend == 0 ) 00310 { 00311 ypp[0] = 0.0; 00312 ypp[1] = 0.0; 00313 } 00314 else 00315 { 00316 if( d3_np_fs( n, a, b, ypp ) == NULL ) 00317 { 00318 fprintf( ioQQQ, "SPLINE_CUBIC_SET - Fatal error!\n" ); 00319 fprintf( ioQQQ, " The linear system could not be solved.\n" ); 00320 cdEXIT(EXIT_FAILURE); 00321 } 00322 } 00323 00324 free(b); 00325 free(a); 00326 return; 00327 } 00328 00329 //****************************************************************************80 00330 00331 void spline_cubic_val( long n, const double t[], double tval, const double y[], const double ypp[], 00332 double *yval, double *ypval, double *yppval ) 00333 00334 //****************************************************************************80 00335 // 00336 // Purpose: 00337 // 00338 // SPLINE_CUBIC_VAL evaluates a piecewise cubic spline at a point. 00339 // 00340 // Discussion: 00341 // 00342 // SPLINE_CUBIC_SET must have already been called to define the values of YPP. 00343 // 00344 // For any point T in the interval T(IVAL), T(IVAL+1), the form of 00345 // the spline is 00346 // 00347 // SPL(T) = A 00348 // + B * ( T - T(IVAL) ) 00349 // + C * ( T - T(IVAL) )**2 00350 // + D * ( T - T(IVAL) )**3 00351 // 00352 // Here: 00353 // A = Y(IVAL) 00354 // B = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) ) 00355 // - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6 00356 // C = YPP(IVAL) / 2 00357 // D = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) ) 00358 // 00359 // Modified: 00360 // 00361 // 04 February 1999 00362 // 00363 // Author: 00364 // 00365 // John Burkardt 00366 // 00367 // Parameters: 00368 // 00369 // Input, long N, the number of knots. 00370 // 00371 // Input, double Y[N], the data values at the knots. 00372 // 00373 // Input, double T[N], the knot values. 00374 // 00375 // Input, double TVAL, a point, typically between T[0] and T[N-1], at 00376 // which the spline is to be evalulated. If TVAL lies outside 00377 // this range, extrapolation is used. 00378 // 00379 // Input, double Y[N], the data values at the knots. 00380 // 00381 // Input, double YPP[N], the second derivatives of the spline at 00382 // the knots. 00383 // 00384 // Output, double *YPVAL, the derivative of the spline at TVAL. 00385 // 00386 // Output, double *YPPVAL, the second derivative of the spline at TVAL. 00387 // 00388 // Output, double SPLINE_VAL, the value of the spline at TVAL. 00389 // 00390 { 00391 DEBUG_ENTRY( "spline_cubic_val()" ); 00392 // 00393 // Determine the interval [ T(I), T(I+1) ] that contains TVAL. 00394 // Values below T[0] or above T[N-1] use extrapolation. 00395 // 00396 long ival = hunt_bisect( t, n, tval ); 00397 // 00398 // In the interval I, the polynomial is in terms of a normalized 00399 // coordinate between 0 and 1. 00400 // 00401 double dt = tval - t[ival]; 00402 double h = t[ival+1] - t[ival]; 00403 00404 if( yval != NULL ) 00405 { 00406 *yval = y[ival] 00407 + dt * ( ( y[ival+1] - y[ival] ) / h 00408 - ( ypp[ival+1] / 6.0 + ypp[ival] / 3.0 ) * h 00409 + dt * ( 0.5 * ypp[ival] 00410 + dt * ( ( ypp[ival+1] - ypp[ival] ) / ( 6.0 * h ) ) ) ); 00411 } 00412 if( ypval != NULL ) 00413 { 00414 *ypval = ( y[ival+1] - y[ival] ) / h 00415 - ( ypp[ival+1] / 6.0 + ypp[ival] / 3.0 ) * h 00416 + dt * ( ypp[ival] 00417 + dt * ( 0.5 * ( ypp[ival+1] - ypp[ival] ) / h ) ); 00418 } 00419 if( yppval != NULL ) 00420 { 00421 *yppval = ypp[ival] + dt * ( ypp[ival+1] - ypp[ival] ) / h; 00422 } 00423 return; 00424 } 00425 00426 //****************************************************************************80 00427 // 00428 // the routines lagrange and linint where written by Peter van Hoof (ROB) 00429 // 00430 //****************************************************************************80 00431 00433 double lagrange(const double x[], /* x[n] */ 00434 const double y[], /* y[n] */ 00435 long n, 00436 double xval) 00437 { 00438 double yval = 0.; 00439 00440 DEBUG_ENTRY( "lagrange()" ); 00441 00442 for( long i=0; i < n; i++ ) 00443 { 00444 double l = 1.; 00445 for( long j=0; j < n; j++ ) 00446 { 00447 if( i != j ) 00448 l *= (xval-x[j])/(x[i]-x[j]); 00449 } 00450 yval += y[i]*l; 00451 } 00452 return yval; 00453 } 00454 00456 double linint(const double x[], /* x[n] */ 00457 const double y[], /* y[n] */ 00458 long n, 00459 double xval) 00460 { 00461 double yval; 00462 00463 DEBUG_ENTRY( "linint()" ); 00464 00465 ASSERT( n >= 2 ); 00466 00467 if( xval <= x[0] ) 00468 yval = y[0]; 00469 else if( xval >= x[n-1] ) 00470 yval = y[n-1]; 00471 else 00472 { 00473 long ilo = hunt_bisect( x, n, xval ); 00474 double deriv = (y[ilo+1]-y[ilo])/(x[ilo+1]-x[ilo]); 00475 yval = y[ilo] + deriv*(xval-x[ilo]); 00476 } 00477 return yval; 00478 }