cloudy trunk
Loading...
Searching...
No Matches
thirdparty_interpolate.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3#include "cddefines.h"
4#include "thirdparty.h"
5
6STATIC double *d3_np_fs( long n, double a[], const double b[], double x[] );
7
8//****************************************************************************80
9//
10// The routines d3_np_fs, spline_cubic_set, spline_cubic_val where written
11// by John Burkardt (Computer Science Department, Florida State University)
12// and have been slightly modified and adapted for use in Cloudy by Peter
13// van Hoof (Royal Observatory of Belgium).
14//
15// The original sources can be found at
16// http://www.scs.fsu.edu/~burkardt/cpp_src/spline/spline.html
17//
18//****************************************************************************80
19
20STATIC double *d3_np_fs( long n, double a[], const double b[], double x[] )
21
22//****************************************************************************80
23//
24// Purpose:
25//
26// D3_NP_FS factors and solves a D3 system.
27//
28// Discussion:
29//
30// The D3 storage format is used for a tridiagonal matrix.
31// The superdiagonal is stored in entries (1,2:N), the diagonal in
32// entries (2,1:N), and the subdiagonal in (3,1:N-1). Thus, the
33// original matrix is "collapsed" vertically into the array.
34//
35// This algorithm requires that each diagonal entry be nonzero.
36// It does not use pivoting, and so can fail on systems that
37// are actually nonsingular.
38//
39// Example:
40//
41// Here is how a D3 matrix of order 5 would be stored:
42//
43// * A12 A23 A34 A45
44// A11 A22 A33 A44 A55
45// A21 A32 A43 A54 *
46//
47// Modified:
48//
49// 15 November 2003
50//
51// Author:
52//
53// John Burkardt
54//
55// Parameters:
56//
57// Input, long N, the order of the linear system.
58//
59// Input/output, double A[3*N].
60// On input, the nonzero diagonals of the linear system.
61// On output, the data in these vectors has been overwritten
62// by factorization information.
63//
64// Input, double B[N], the right hand side.
65//
66// Output, double X[N] and D3_NP_FS[N], the solution of the linear system.
67// D3_NP_FS returns NULL if there was an error because one of the diagonal
68// entries was zero.
69//
70{
71 long i;
72 double xmult;
73
74 DEBUG_ENTRY( "d3_np_fs()" );
75//
76// Check.
77//
78 for( i = 0; i < n; i++ )
79 {
80 if( a[1+i*3] == 0.0 )
81 {
82 return NULL;
83 }
84 }
85
86 x[0] = b[0];
87
88 for( i = 1; i < n; i++ )
89 {
90 xmult = a[2+(i-1)*3] / a[1+(i-1)*3];
91 a[1+i*3] = a[1+i*3] - xmult * a[0+i*3];
92 x[i] = b[i] - xmult * x[i-1];
93 }
94
95 x[n-1] = x[n-1] / a[1+(n-1)*3];
96 for( i = n-2; 0 <= i; i-- )
97 {
98 x[i] = ( x[i] - a[0+(i+1)*3] * x[i+1] ) / a[1+i*3];
99 }
100 return x;
101}
102
103//****************************************************************************80
104
105void spline_cubic_set( long n, const double t[], const double y[], double ypp[],
106 int ibcbeg, double ybcbeg, int ibcend, double ybcend )
107
108//****************************************************************************80
109//
110// Purpose:
111//
112// SPLINE_CUBIC_SET computes the second derivatives of a piecewise cubic spline.
113//
114// Discussion:
115//
116// For data interpolation, the user must call SPLINE_SET to determine
117// the second derivative data, passing in the data to be interpolated,
118// and the desired boundary conditions.
119//
120// The data to be interpolated, plus the SPLINE_SET output, defines
121// the spline. The user may then call SPLINE_VAL to evaluate the
122// spline at any point.
123//
124// The cubic spline is a piecewise cubic polynomial. The intervals
125// are determined by the "knots" or abscissas of the data to be
126// interpolated. The cubic spline has continous first and second
127// derivatives over the entire interval of interpolation.
128//
129// For any point T in the interval T(IVAL), T(IVAL+1), the form of
130// the spline is
131//
132// SPL(T) = A(IVAL)
133// + B(IVAL) * ( T - T(IVAL) )
134// + C(IVAL) * ( T - T(IVAL) )**2
135// + D(IVAL) * ( T - T(IVAL) )**3
136//
137// If we assume that we know the values Y(*) and YPP(*), which represent
138// the values and second derivatives of the spline at each knot, then
139// the coefficients can be computed as:
140//
141// A(IVAL) = Y(IVAL)
142// B(IVAL) = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
143// - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
144// C(IVAL) = YPP(IVAL) / 2
145// D(IVAL) = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
146//
147// Since the first derivative of the spline is
148//
149// SPL'(T) = B(IVAL)
150// + 2 * C(IVAL) * ( T - T(IVAL) )
151// + 3 * D(IVAL) * ( T - T(IVAL) )**2,
152//
153// the requirement that the first derivative be continuous at interior
154// knot I results in a total of N-2 equations, of the form:
155//
156// B(IVAL-1) + 2 C(IVAL-1) * (T(IVAL)-T(IVAL-1))
157// + 3 * D(IVAL-1) * (T(IVAL) - T(IVAL-1))**2 = B(IVAL)
158//
159// or, setting H(IVAL) = T(IVAL+1) - T(IVAL)
160//
161// ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
162// - ( YPP(IVAL) + 2 * YPP(IVAL-1) ) * H(IVAL-1) / 6
163// + YPP(IVAL-1) * H(IVAL-1)
164// + ( YPP(IVAL) - YPP(IVAL-1) ) * H(IVAL-1) / 2
165// =
166// ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
167// - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * H(IVAL) / 6
168//
169// or
170//
171// YPP(IVAL-1) * H(IVAL-1) + 2 * YPP(IVAL) * ( H(IVAL-1) + H(IVAL) )
172// + YPP(IVAL) * H(IVAL)
173// =
174// 6 * ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
175// - 6 * ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
176//
177// Boundary conditions must be applied at the first and last knots.
178// The resulting tridiagonal system can be solved for the YPP values.
179//
180// Modified:
181//
182// 06 February 2004
183//
184// Author:
185//
186// John Burkardt
187//
188// Parameters:
189//
190// Input, long N, the number of data points. N must be at least 2.
191// In the special case where N = 2 and IBCBEG = IBCEND = 0, the
192// spline will actually be linear.
193//
194// Input, double T[N], the knot values, that is, the points were data is
195// specified. The knot values should be distinct, and increasing.
196//
197// Input, double Y[N], the data values to be interpolated.
198//
199// Input, int IBCBEG, left boundary condition flag:
200// 0: the cubic spline should be a quadratic over the first interval;
201// 1: the first derivative at the left endpoint should be YBCBEG;
202// 2: the second derivative at the left endpoint should be YBCBEG.
203//
204// Input, double YBCBEG, the values to be used in the boundary
205// conditions if IBCBEG is equal to 1 or 2.
206//
207// Input, int IBCEND, right boundary condition flag:
208// 0: the cubic spline should be a quadratic over the last interval;
209// 1: the first derivative at the right endpoint should be YBCEND;
210// 2: the second derivative at the right endpoint should be YBCEND.
211//
212// Input, double YBCEND, the values to be used in the boundary
213// conditions if IBCEND is equal to 1 or 2.
214//
215// Output, double YPP[N] and SPLINE_CUBIC_SET[N], the second derivatives
216// of the cubic spline.
217//
218{
219 long i;
220
221 DEBUG_ENTRY( "spline_cubic_set()" );
222//
223// Check.
224//
225 ASSERT( n >= 2 );
226
227# ifndef NDEBUG
228 for( i = 0; i < n - 1; i++ )
229 {
230 if( t[i+1] <= t[i] )
231 {
232 fprintf( ioQQQ, "SPLINE_CUBIC_SET - Fatal error!\n" );
233 fprintf( ioQQQ, " The knots must be strictly increasing\n" );
235 }
236 }
237# endif
238
239 double *a = (double*)MALLOC((size_t)(3*n)*sizeof(double));
240 double *b = (double*)MALLOC((size_t)n*sizeof(double));
241//
242// Set up the first equation.
243//
244 if( ibcbeg == 0 )
245 {
246 b[0] = 0.0;
247 a[1+0*3] = 1.0;
248 a[0+1*3] = -1.0;
249 }
250 else if( ibcbeg == 1 )
251 {
252 b[0] = ( y[1] - y[0] ) / ( t[1] - t[0] ) - ybcbeg;
253 a[1+0*3] = ( t[1] - t[0] ) / 3.0;
254 a[0+1*3] = ( t[1] - t[0] ) / 6.0;
255 }
256 else if( ibcbeg == 2 )
257 {
258 b[0] = ybcbeg;
259 a[1+0*3] = 1.0;
260 a[0+1*3] = 0.0;
261 }
262 else
263 {
264 fprintf( ioQQQ, "SPLINE_CUBIC_SET - Fatal error!\n" );
265 fprintf( ioQQQ, " IBCBEG must be 0, 1 or 2, but I found %d.\n", ibcbeg );
267 }
268//
269// Set up the intermediate equations.
270//
271 for( i = 1; i < n-1; i++ )
272 {
273 b[i] = ( y[i+1] - y[i] ) / ( t[i+1] - t[i] )
274 - ( y[i] - y[i-1] ) / ( t[i] - t[i-1] );
275 a[2+(i-1)*3] = ( t[i] - t[i-1] ) / 6.0;
276 a[1+ i *3] = ( t[i+1] - t[i-1] ) / 3.0;
277 a[0+(i+1)*3] = ( t[i+1] - t[i] ) / 6.0;
278 }
279//
280// Set up the last equation.
281//
282 if( ibcend == 0 )
283 {
284 b[n-1] = 0.0;
285 a[2+(n-2)*3] = -1.0;
286 a[1+(n-1)*3] = 1.0;
287 }
288 else if( ibcend == 1 )
289 {
290 b[n-1] = ybcend - ( y[n-1] - y[n-2] ) / ( t[n-1] - t[n-2] );
291 a[2+(n-2)*3] = ( t[n-1] - t[n-2] ) / 6.0;
292 a[1+(n-1)*3] = ( t[n-1] - t[n-2] ) / 3.0;
293 }
294 else if( ibcend == 2 )
295 {
296 b[n-1] = ybcend;
297 a[2+(n-2)*3] = 0.0;
298 a[1+(n-1)*3] = 1.0;
299 }
300 else
301 {
302 fprintf( ioQQQ, "SPLINE_CUBIC_SET - Fatal error!\n" );
303 fprintf( ioQQQ, " IBCEND must be 0, 1 or 2, but I found %d.\n", ibcend );
305 }
306//
307// Solve the linear system.
308//
309 if( n == 2 && ibcbeg == 0 && ibcend == 0 )
310 {
311 ypp[0] = 0.0;
312 ypp[1] = 0.0;
313 }
314 else
315 {
316 if( d3_np_fs( n, a, b, ypp ) == NULL )
317 {
318 fprintf( ioQQQ, "SPLINE_CUBIC_SET - Fatal error!\n" );
319 fprintf( ioQQQ, " The linear system could not be solved.\n" );
321 }
322 }
323
324 free(b);
325 free(a);
326 return;
327}
328
329//****************************************************************************80
330
331void spline_cubic_val( long n, const double t[], double tval, const double y[], const double ypp[],
332 double *yval, double *ypval, double *yppval )
333
334//****************************************************************************80
335//
336// Purpose:
337//
338// SPLINE_CUBIC_VAL evaluates a piecewise cubic spline at a point.
339//
340// Discussion:
341//
342// SPLINE_CUBIC_SET must have already been called to define the values of YPP.
343//
344// For any point T in the interval T(IVAL), T(IVAL+1), the form of
345// the spline is
346//
347// SPL(T) = A
348// + B * ( T - T(IVAL) )
349// + C * ( T - T(IVAL) )**2
350// + D * ( T - T(IVAL) )**3
351//
352// Here:
353// A = Y(IVAL)
354// B = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
355// - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
356// C = YPP(IVAL) / 2
357// D = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
358//
359// Modified:
360//
361// 04 February 1999
362//
363// Author:
364//
365// John Burkardt
366//
367// Parameters:
368//
369// Input, long N, the number of knots.
370//
371// Input, double Y[N], the data values at the knots.
372//
373// Input, double T[N], the knot values.
374//
375// Input, double TVAL, a point, typically between T[0] and T[N-1], at
376// which the spline is to be evalulated. If TVAL lies outside
377// this range, extrapolation is used.
378//
379// Input, double Y[N], the data values at the knots.
380//
381// Input, double YPP[N], the second derivatives of the spline at
382// the knots.
383//
384// Output, double *YPVAL, the derivative of the spline at TVAL.
385//
386// Output, double *YPPVAL, the second derivative of the spline at TVAL.
387//
388// Output, double SPLINE_VAL, the value of the spline at TVAL.
389//
390{
391 DEBUG_ENTRY( "spline_cubic_val()" );
392//
393// Determine the interval [ T(I), T(I+1) ] that contains TVAL.
394// Values below T[0] or above T[N-1] use extrapolation.
395//
396 long ival = hunt_bisect( t, n, tval );
397//
398// In the interval I, the polynomial is in terms of a normalized
399// coordinate between 0 and 1.
400//
401 double dt = tval - t[ival];
402 double h = t[ival+1] - t[ival];
403
404 if( yval != NULL )
405 {
406 *yval = y[ival]
407 + dt * ( ( y[ival+1] - y[ival] ) / h
408 - ( ypp[ival+1] / 6.0 + ypp[ival] / 3.0 ) * h
409 + dt * ( 0.5 * ypp[ival]
410 + dt * ( ( ypp[ival+1] - ypp[ival] ) / ( 6.0 * h ) ) ) );
411 }
412 if( ypval != NULL )
413 {
414 *ypval = ( y[ival+1] - y[ival] ) / h
415 - ( ypp[ival+1] / 6.0 + ypp[ival] / 3.0 ) * h
416 + dt * ( ypp[ival]
417 + dt * ( 0.5 * ( ypp[ival+1] - ypp[ival] ) / h ) );
418 }
419 if( yppval != NULL )
420 {
421 *yppval = ypp[ival] + dt * ( ypp[ival+1] - ypp[ival] ) / h;
422 }
423 return;
424}
425
426//****************************************************************************80
427//
428// the routines lagrange and linint where written by Peter van Hoof (ROB)
429//
430//****************************************************************************80
431
433double lagrange(const double x[], /* x[n] */
434 const double y[], /* y[n] */
435 long n,
436 double xval)
437{
438 double yval = 0.;
439
440 DEBUG_ENTRY( "lagrange()" );
441
442 for( long i=0; i < n; i++ )
443 {
444 double l = 1.;
445 for( long j=0; j < n; j++ )
446 {
447 if( i != j )
448 l *= (xval-x[j])/(x[i]-x[j]);
449 }
450 yval += y[i]*l;
451 }
452 return yval;
453}
454
456double linint(const double x[], /* x[n] */
457 const double y[], /* y[n] */
458 long n,
459 double xval)
460{
461 double yval;
462
463 DEBUG_ENTRY( "linint()" );
464
465 ASSERT( n >= 2 );
466
467 if( xval <= x[0] )
468 yval = y[0];
469 else if( xval >= x[n-1] )
470 yval = y[n-1];
471 else
472 {
473 long ilo = hunt_bisect( x, n, xval );
474 double deriv = (y[ilo+1]-y[ilo])/(x[ilo+1]-x[ilo]);
475 yval = y[ilo] + deriv*(xval-x[ilo]);
476 }
477 return yval;
478}
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define MALLOC(exp)
Definition cddefines.h:501
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long hunt_bisect(const T x[], long n, T xval)
Definition thirdparty.h:270
STATIC double * d3_np_fs(long n, double a[], const double b[], double x[])
double linint(const double x[], const double y[], long n, double xval)
double lagrange(const double x[], const double y[], long n, double xval)
void spline_cubic_set(long n, const double t[], const double y[], double ypp[], int ibcbeg, double ybcbeg, int ibcend, double ybcend)
void spline_cubic_val(long n, const double t[], double tval, const double y[], const double ypp[], double *yval, double *ypval, double *yppval)