00001 /* This file is part of Cloudy and is copyright (C)1978-2007 by Gary J. Ferland 00002 * 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 DEBUG_EXIT( "d3_np_fs()" ); 00083 return NULL; 00084 } 00085 } 00086 00087 x[0] = b[0]; 00088 00089 for( i = 1; i < n; i++ ) 00090 { 00091 xmult = a[2+(i-1)*3] / a[1+(i-1)*3]; 00092 a[1+i*3] = a[1+i*3] - xmult * a[0+i*3]; 00093 x[i] = b[i] - xmult * x[i-1]; 00094 } 00095 00096 x[n-1] = x[n-1] / a[1+(n-1)*3]; 00097 for( i = n-2; 0 <= i; i-- ) 00098 { 00099 x[i] = ( x[i] - a[0+(i+1)*3] * x[i+1] ) / a[1+i*3]; 00100 } 00101 00102 DEBUG_EXIT( "d3_np_fs()" ); 00103 return x; 00104 } 00105 00106 //****************************************************************************80 00107 00108 void spline_cubic_set( long n, const double t[], const double y[], double ypp[], 00109 int ibcbeg, double ybcbeg, int ibcend, double ybcend ) 00110 00111 //****************************************************************************80 00112 // 00113 // Purpose: 00114 // 00115 // SPLINE_CUBIC_SET computes the second derivatives of a piecewise cubic spline. 00116 // 00117 // Discussion: 00118 // 00119 // For data interpolation, the user must call SPLINE_SET to determine 00120 // the second derivative data, passing in the data to be interpolated, 00121 // and the desired boundary conditions. 00122 // 00123 // The data to be interpolated, plus the SPLINE_SET output, defines 00124 // the spline. The user may then call SPLINE_VAL to evaluate the 00125 // spline at any point. 00126 // 00127 // The cubic spline is a piecewise cubic polynomial. The intervals 00128 // are determined by the "knots" or abscissas of the data to be 00129 // interpolated. The cubic spline has continous first and second 00130 // derivatives over the entire interval of interpolation. 00131 // 00132 // For any point T in the interval T(IVAL), T(IVAL+1), the form of 00133 // the spline is 00134 // 00135 // SPL(T) = A(IVAL) 00136 // + B(IVAL) * ( T - T(IVAL) ) 00137 // + C(IVAL) * ( T - T(IVAL) )**2 00138 // + D(IVAL) * ( T - T(IVAL) )**3 00139 // 00140 // If we assume that we know the values Y(*) and YPP(*), which represent 00141 // the values and second derivatives of the spline at each knot, then 00142 // the coefficients can be computed as: 00143 // 00144 // A(IVAL) = Y(IVAL) 00145 // B(IVAL) = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) ) 00146 // - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6 00147 // C(IVAL) = YPP(IVAL) / 2 00148 // D(IVAL) = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) ) 00149 // 00150 // Since the first derivative of the spline is 00151 // 00152 // SPL'(T) = B(IVAL) 00153 // + 2 * C(IVAL) * ( T - T(IVAL) ) 00154 // + 3 * D(IVAL) * ( T - T(IVAL) )**2, 00155 // 00156 // the requirement that the first derivative be continuous at interior 00157 // knot I results in a total of N-2 equations, of the form: 00158 // 00159 // B(IVAL-1) + 2 C(IVAL-1) * (T(IVAL)-T(IVAL-1)) 00160 // + 3 * D(IVAL-1) * (T(IVAL) - T(IVAL-1))**2 = B(IVAL) 00161 // 00162 // or, setting H(IVAL) = T(IVAL+1) - T(IVAL) 00163 // 00164 // ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1) 00165 // - ( YPP(IVAL) + 2 * YPP(IVAL-1) ) * H(IVAL-1) / 6 00166 // + YPP(IVAL-1) * H(IVAL-1) 00167 // + ( YPP(IVAL) - YPP(IVAL-1) ) * H(IVAL-1) / 2 00168 // = 00169 // ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL) 00170 // - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * H(IVAL) / 6 00171 // 00172 // or 00173 // 00174 // YPP(IVAL-1) * H(IVAL-1) + 2 * YPP(IVAL) * ( H(IVAL-1) + H(IVAL) ) 00175 // + YPP(IVAL) * H(IVAL) 00176 // = 00177 // 6 * ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL) 00178 // - 6 * ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1) 00179 // 00180 // Boundary conditions must be applied at the first and last knots. 00181 // The resulting tridiagonal system can be solved for the YPP values. 00182 // 00183 // Modified: 00184 // 00185 // 06 February 2004 00186 // 00187 // Author: 00188 // 00189 // John Burkardt 00190 // 00191 // Parameters: 00192 // 00193 // Input, long N, the number of data points. N must be at least 2. 00194 // In the special case where N = 2 and IBCBEG = IBCEND = 0, the 00195 // spline will actually be linear. 00196 // 00197 // Input, double T[N], the knot values, that is, the points were data is 00198 // specified. The knot values should be distinct, and increasing. 00199 // 00200 // Input, double Y[N], the data values to be interpolated. 00201 // 00202 // Input, int IBCBEG, left boundary condition flag: 00203 // 0: the cubic spline should be a quadratic over the first interval; 00204 // 1: the first derivative at the left endpoint should be YBCBEG; 00205 // 2: the second derivative at the left endpoint should be YBCBEG. 00206 // 00207 // Input, double YBCBEG, the values to be used in the boundary 00208 // conditions if IBCBEG is equal to 1 or 2. 00209 // 00210 // Input, int IBCEND, right boundary condition flag: 00211 // 0: the cubic spline should be a quadratic over the last interval; 00212 // 1: the first derivative at the right endpoint should be YBCEND; 00213 // 2: the second derivative at the right endpoint should be YBCEND. 00214 // 00215 // Input, double YBCEND, the values to be used in the boundary 00216 // conditions if IBCEND is equal to 1 or 2. 00217 // 00218 // Output, double YPP[N] and SPLINE_CUBIC_SET[N], the second derivatives 00219 // of the cubic spline. 00220 // 00221 { 00222 long i; 00223 00224 DEBUG_ENTRY( "spline_cubic_set()" ); 00225 // 00226 // Check. 00227 // 00228 ASSERT( n >= 2 ); 00229 00230 # ifndef NDEBUG 00231 for( i = 0; i < n - 1; i++ ) 00232 { 00233 if( t[i+1] <= t[i] ) 00234 { 00235 fprintf( ioQQQ, "SPLINE_CUBIC_SET - Fatal error!\n" ); 00236 fprintf( ioQQQ, " The knots must be strictly increasing\n" ); 00237 puts( "[Stop in spline_cubic_set]" ); 00238 cdEXIT(EXIT_FAILURE); 00239 } 00240 } 00241 # endif 00242 00243 double *a = (double*)MALLOC((size_t)(3*n)*sizeof(double)); 00244 double *b = (double*)MALLOC((size_t)n*sizeof(double)); 00245 // 00246 // Set up the first equation. 00247 // 00248 if( ibcbeg == 0 ) 00249 { 00250 b[0] = 0.0; 00251 a[1+0*3] = 1.0; 00252 a[0+1*3] = -1.0; 00253 } 00254 else if( ibcbeg == 1 ) 00255 { 00256 b[0] = ( y[1] - y[0] ) / ( t[1] - t[0] ) - ybcbeg; 00257 a[1+0*3] = ( t[1] - t[0] ) / 3.0; 00258 a[0+1*3] = ( t[1] - t[0] ) / 6.0; 00259 } 00260 else if( ibcbeg == 2 ) 00261 { 00262 b[0] = ybcbeg; 00263 a[1+0*3] = 1.0; 00264 a[0+1*3] = 0.0; 00265 } 00266 else 00267 { 00268 fprintf( ioQQQ, "SPLINE_CUBIC_SET - Fatal error!\n" ); 00269 fprintf( ioQQQ, " IBCBEG must be 0, 1 or 2, but I found %d.\n", ibcbeg ); 00270 puts( "[Stop in spline_cubic_set]" ); 00271 cdEXIT(EXIT_FAILURE); 00272 } 00273 // 00274 // Set up the intermediate equations. 00275 // 00276 for( i = 1; i < n-1; i++ ) 00277 { 00278 b[i] = ( y[i+1] - y[i] ) / ( t[i+1] - t[i] ) 00279 - ( y[i] - y[i-1] ) / ( t[i] - t[i-1] ); 00280 a[2+(i-1)*3] = ( t[i] - t[i-1] ) / 6.0; 00281 a[1+ i *3] = ( t[i+1] - t[i-1] ) / 3.0; 00282 a[0+(i+1)*3] = ( t[i+1] - t[i] ) / 6.0; 00283 } 00284 // 00285 // Set up the last equation. 00286 // 00287 if( ibcend == 0 ) 00288 { 00289 b[n-1] = 0.0; 00290 a[2+(n-2)*3] = -1.0; 00291 a[1+(n-1)*3] = 1.0; 00292 } 00293 else if( ibcend == 1 ) 00294 { 00295 b[n-1] = ybcend - ( y[n-1] - y[n-2] ) / ( t[n-1] - t[n-2] ); 00296 a[2+(n-2)*3] = ( t[n-1] - t[n-2] ) / 6.0; 00297 a[1+(n-1)*3] = ( t[n-1] - t[n-2] ) / 3.0; 00298 } 00299 else if( ibcend == 2 ) 00300 { 00301 b[n-1] = ybcend; 00302 a[2+(n-2)*3] = 0.0; 00303 a[1+(n-1)*3] = 1.0; 00304 } 00305 else 00306 { 00307 fprintf( ioQQQ, "SPLINE_CUBIC_SET - Fatal error!\n" ); 00308 fprintf( ioQQQ, " IBCEND must be 0, 1 or 2, but I found %d.\n", ibcend ); 00309 puts( "[Stop in spline_cubic_set]" ); 00310 cdEXIT(EXIT_FAILURE); 00311 } 00312 // 00313 // Solve the linear system. 00314 // 00315 if( n == 2 && ibcbeg == 0 && ibcend == 0 ) 00316 { 00317 ypp[0] = 0.0; 00318 ypp[1] = 0.0; 00319 } 00320 else 00321 { 00322 if( d3_np_fs( n, a, b, ypp ) == NULL ) 00323 { 00324 fprintf( ioQQQ, "SPLINE_CUBIC_SET - Fatal error!\n" ); 00325 fprintf( ioQQQ, " The linear system could not be solved.\n" ); 00326 puts( "[Stop in spline_cubic_set]" ); 00327 cdEXIT(EXIT_FAILURE); 00328 } 00329 } 00330 00331 free(b); 00332 free(a); 00333 00334 DEBUG_EXIT( "spline_cubic_set()" ); 00335 return; 00336 } 00337 00338 //****************************************************************************80 00339 00340 void spline_cubic_val( long n, const double t[], double tval, const double y[], const double ypp[], 00341 double *yval, double *ypval, double *yppval ) 00342 00343 //****************************************************************************80 00344 // 00345 // Purpose: 00346 // 00347 // SPLINE_CUBIC_VAL evaluates a piecewise cubic spline at a point. 00348 // 00349 // Discussion: 00350 // 00351 // SPLINE_CUBIC_SET must have already been called to define the values of YPP. 00352 // 00353 // For any point T in the interval T(IVAL), T(IVAL+1), the form of 00354 // the spline is 00355 // 00356 // SPL(T) = A 00357 // + B * ( T - T(IVAL) ) 00358 // + C * ( T - T(IVAL) )**2 00359 // + D * ( T - T(IVAL) )**3 00360 // 00361 // Here: 00362 // A = Y(IVAL) 00363 // B = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) ) 00364 // - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6 00365 // C = YPP(IVAL) / 2 00366 // D = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) ) 00367 // 00368 // Modified: 00369 // 00370 // 04 February 1999 00371 // 00372 // Author: 00373 // 00374 // John Burkardt 00375 // 00376 // Parameters: 00377 // 00378 // Input, long N, the number of knots. 00379 // 00380 // Input, double Y[N], the data values at the knots. 00381 // 00382 // Input, double T[N], the knot values. 00383 // 00384 // Input, double TVAL, a point, typically between T[0] and T[N-1], at 00385 // which the spline is to be evalulated. If TVAL lies outside 00386 // this range, extrapolation is used. 00387 // 00388 // Input, double Y[N], the data values at the knots. 00389 // 00390 // Input, double YPP[N], the second derivatives of the spline at 00391 // the knots. 00392 // 00393 // Output, double *YPVAL, the derivative of the spline at TVAL. 00394 // 00395 // Output, double *YPPVAL, the second derivative of the spline at TVAL. 00396 // 00397 // Output, double SPLINE_VAL, the value of the spline at TVAL. 00398 // 00399 { 00400 DEBUG_ENTRY( "spline_cubic_val()" ); 00401 // 00402 // Determine the interval [ T(I), T(I+1) ] that contains TVAL. 00403 // Values below T[0] or above T[N-1] use extrapolation. 00404 // 00405 long ival = hunt_bisect( t, n, tval ); 00406 // 00407 // In the interval I, the polynomial is in terms of a normalized 00408 // coordinate between 0 and 1. 00409 // 00410 double dt = tval - t[ival]; 00411 double h = t[ival+1] - t[ival]; 00412 00413 if( yval != NULL ) 00414 { 00415 *yval = y[ival] 00416 + dt * ( ( y[ival+1] - y[ival] ) / h 00417 - ( ypp[ival+1] / 6.0 + ypp[ival] / 3.0 ) * h 00418 + dt * ( 0.5 * ypp[ival] 00419 + dt * ( ( ypp[ival+1] - ypp[ival] ) / ( 6.0 * h ) ) ) ); 00420 } 00421 if( ypval != NULL ) 00422 { 00423 *ypval = ( y[ival+1] - y[ival] ) / h 00424 - ( ypp[ival+1] / 6.0 + ypp[ival] / 3.0 ) * h 00425 + dt * ( ypp[ival] 00426 + dt * ( 0.5 * ( ypp[ival+1] - ypp[ival] ) / h ) ); 00427 } 00428 if( yppval != NULL ) 00429 { 00430 *yppval = ypp[ival] + dt * ( ypp[ival+1] - ypp[ival] ) / h; 00431 } 00432 00433 DEBUG_EXIT( "spline_cubic_val()" ); 00434 return; 00435 } 00436 00437 //****************************************************************************80 00438 // 00439 // the routines lagrange and linint where written by Peter van Hoof (ROB) 00440 // 00441 //****************************************************************************80 00442 00444 double lagrange(const double x[], /* x[n] */ 00445 const double y[], /* y[n] */ 00446 long n, 00447 double xval) 00448 { 00449 double yval = 0.; 00450 00451 DEBUG_ENTRY( "lagrange()" ); 00452 00453 for( long i=0; i < n; i++ ) 00454 { 00455 double l = 1.; 00456 for( long j=0; j < n; j++ ) 00457 { 00458 if( i != j ) 00459 l *= (xval-x[j])/(x[i]-x[j]); 00460 } 00461 yval += y[i]*l; 00462 } 00463 00464 DEBUG_EXIT( "lagrange()" ); 00465 return yval; 00466 } 00467 00469 double linint(const double x[], /* x[n] */ 00470 const double y[], /* y[n] */ 00471 long n, 00472 double xval) 00473 { 00474 double yval; 00475 00476 DEBUG_ENTRY( "linint()" ); 00477 00478 ASSERT( n >= 2 ); 00479 00480 if( xval <= x[0] ) 00481 yval = y[0]; 00482 else if( xval >= x[n-1] ) 00483 yval = y[n-1]; 00484 else 00485 { 00486 long ilo = hunt_bisect( x, n, xval ); 00487 double deriv = (y[ilo+1]-y[ilo])/(x[ilo+1]-x[ilo]); 00488 yval = y[ilo] + deriv*(xval-x[ilo]); 00489 } 00490 00491 DEBUG_EXIT( "linint()" ); 00492 return yval; 00493 }