cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
thirdparty_interpolate.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 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 
6 STATIC 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 
20 STATIC 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 
105 void 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  vector<double> a(3*n), b(n);
240 //
241 // Set up the first equation.
242 //
243  if( ibcbeg == 0 )
244  {
245  b[0] = 0.0;
246  a[1+0*3] = 1.0;
247  a[0+1*3] = -1.0;
248  }
249  else if( ibcbeg == 1 )
250  {
251  b[0] = ( y[1] - y[0] ) / ( t[1] - t[0] ) - ybcbeg;
252  a[1+0*3] = ( t[1] - t[0] ) / 3.0;
253  a[0+1*3] = ( t[1] - t[0] ) / 6.0;
254  }
255  else if( ibcbeg == 2 )
256  {
257  b[0] = ybcbeg;
258  a[1+0*3] = 1.0;
259  a[0+1*3] = 0.0;
260  }
261  else
262  {
263  fprintf( ioQQQ, "SPLINE_CUBIC_SET - Fatal error!\n" );
264  fprintf( ioQQQ, " IBCBEG must be 0, 1 or 2, but I found %d.\n", ibcbeg );
266  }
267 //
268 // Set up the intermediate equations.
269 //
270  for( i = 1; i < n-1; i++ )
271  {
272  b[i] = ( y[i+1] - y[i] ) / ( t[i+1] - t[i] )
273  - ( y[i] - y[i-1] ) / ( t[i] - t[i-1] );
274  a[2+(i-1)*3] = ( t[i] - t[i-1] ) / 6.0;
275  a[1+ i *3] = ( t[i+1] - t[i-1] ) / 3.0;
276  a[0+(i+1)*3] = ( t[i+1] - t[i] ) / 6.0;
277  }
278 //
279 // Set up the last equation.
280 //
281  if( ibcend == 0 )
282  {
283  b[n-1] = 0.0;
284  a[2+(n-2)*3] = -1.0;
285  a[1+(n-1)*3] = 1.0;
286  }
287  else if( ibcend == 1 )
288  {
289  b[n-1] = ybcend - ( y[n-1] - y[n-2] ) / ( t[n-1] - t[n-2] );
290  a[2+(n-2)*3] = ( t[n-1] - t[n-2] ) / 6.0;
291  a[1+(n-1)*3] = ( t[n-1] - t[n-2] ) / 3.0;
292  }
293  else if( ibcend == 2 )
294  {
295  b[n-1] = ybcend;
296  a[2+(n-2)*3] = 0.0;
297  a[1+(n-1)*3] = 1.0;
298  }
299  else
300  {
301  fprintf( ioQQQ, "SPLINE_CUBIC_SET - Fatal error!\n" );
302  fprintf( ioQQQ, " IBCEND must be 0, 1 or 2, but I found %d.\n", ibcend );
304  }
305 //
306 // Solve the linear system.
307 //
308  if( n == 2 && ibcbeg == 0 && ibcend == 0 )
309  {
310  ypp[0] = 0.0;
311  ypp[1] = 0.0;
312  }
313  else
314  {
315  if( d3_np_fs( n, &a[0], &b[0], ypp ) == NULL )
316  {
317  fprintf( ioQQQ, "SPLINE_CUBIC_SET - Fatal error!\n" );
318  fprintf( ioQQQ, " The linear system could not be solved.\n" );
320  }
321  }
322 
323  return;
324 }
325 
326 //****************************************************************************80
327 
328 void spline_cubic_val( long n, const double t[], double tval, const double y[], const double ypp[],
329  double *yval, double *ypval, double *yppval )
330 
331 //****************************************************************************80
332 //
333 // Purpose:
334 //
335 // SPLINE_CUBIC_VAL evaluates a piecewise cubic spline at a point.
336 //
337 // Discussion:
338 //
339 // SPLINE_CUBIC_SET must have already been called to define the values of YPP.
340 //
341 // For any point T in the interval T(IVAL), T(IVAL+1), the form of
342 // the spline is
343 //
344 // SPL(T) = A
345 // + B * ( T - T(IVAL) )
346 // + C * ( T - T(IVAL) )**2
347 // + D * ( T - T(IVAL) )**3
348 //
349 // Here:
350 // A = Y(IVAL)
351 // B = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
352 // - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
353 // C = YPP(IVAL) / 2
354 // D = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
355 //
356 // Modified:
357 //
358 // 04 February 1999
359 //
360 // Author:
361 //
362 // John Burkardt
363 //
364 // Parameters:
365 //
366 // Input, long N, the number of knots.
367 //
368 // Input, double Y[N], the data values at the knots.
369 //
370 // Input, double T[N], the knot values.
371 //
372 // Input, double TVAL, a point, typically between T[0] and T[N-1], at
373 // which the spline is to be evalulated. If TVAL lies outside
374 // this range, extrapolation is used.
375 //
376 // Input, double Y[N], the data values at the knots.
377 //
378 // Input, double YPP[N], the second derivatives of the spline at
379 // the knots.
380 //
381 // Output, double *YPVAL, the derivative of the spline at TVAL.
382 //
383 // Output, double *YPPVAL, the second derivative of the spline at TVAL.
384 //
385 // Output, double SPLINE_VAL, the value of the spline at TVAL.
386 //
387 {
388  DEBUG_ENTRY( "spline_cubic_val()" );
389 //
390 // Determine the interval [ T(I), T(I+1) ] that contains TVAL.
391 // Values below T[0] or above T[N-1] use extrapolation.
392 //
393  long ival = hunt_bisect( t, n, tval );
394 //
395 // In the interval I, the polynomial is in terms of a normalized
396 // coordinate between 0 and 1.
397 //
398  double dt = tval - t[ival];
399  double h = t[ival+1] - t[ival];
400 
401  if( yval != NULL )
402  {
403  *yval = y[ival]
404  + dt * ( ( y[ival+1] - y[ival] ) / h
405  - ( ypp[ival+1] / 6.0 + ypp[ival] / 3.0 ) * h
406  + dt * ( 0.5 * ypp[ival]
407  + dt * ( ( ypp[ival+1] - ypp[ival] ) / ( 6.0 * h ) ) ) );
408  }
409  if( ypval != NULL )
410  {
411  *ypval = ( y[ival+1] - y[ival] ) / h
412  - ( ypp[ival+1] / 6.0 + ypp[ival] / 3.0 ) * h
413  + dt * ( ypp[ival]
414  + dt * ( 0.5 * ( ypp[ival+1] - ypp[ival] ) / h ) );
415  }
416  if( yppval != NULL )
417  {
418  *yppval = ypp[ival] + dt * ( ypp[ival+1] - ypp[ival] ) / h;
419  }
420  return;
421 }
422 
423 //****************************************************************************80
424 //
425 // the routines lagrange and linint where written by Peter van Hoof (ROB)
426 //
427 //****************************************************************************80
428 
430 double lagrange(const double x[], /* x[n] */
431  const double y[], /* y[n] */
432  long n,
433  double xval)
434 {
435  double yval = 0.;
436 
437  DEBUG_ENTRY( "lagrange()" );
438 
439  for( long i=0; i < n; i++ )
440  {
441  double l = 1.;
442  for( long j=0; j < n; j++ )
443  {
444  if( i != j )
445  l *= (xval-x[j])/(x[i]-x[j]);
446  }
447  yval += y[i]*l;
448  }
449  return yval;
450 }
451 
453 double linint(const double x[], /* x[n] */
454  const double y[], /* y[n] */
455  long n,
456  double xval)
457 {
458  double yval;
459 
460  DEBUG_ENTRY( "linint()" );
461 
462  ASSERT( n >= 2 );
463 
464  if( xval <= x[0] )
465  yval = y[0];
466  else if( xval >= x[n-1] )
467  yval = y[n-1];
468  else
469  {
470  long ilo = hunt_bisect( x, n, xval );
471  double deriv = (y[ilo+1]-y[ilo])/(x[ilo+1]-x[ilo]);
472  yval = y[ilo] + deriv*(xval-x[ilo]);
473  }
474  return yval;
475 }
void spline_cubic_val(long n, const double t[], double tval, const double y[], const double ypp[], double *yval, double *ypval, double *yppval)
void spline_cubic_set(long n, const double t[], const double y[], double ypp[], int ibcbeg, double ybcbeg, int ibcend, double ybcend)
FILE * ioQQQ
Definition: cddefines.cpp:7
long hunt_bisect(const T x[], long n, T xval)
Definition: thirdparty.h:303
STATIC double * d3_np_fs(long n, double a[], const double b[], double x[])
#define STATIC
Definition: cddefines.h:118
#define EXIT_FAILURE
Definition: cddefines.h:168
#define cdEXIT(FAIL)
Definition: cddefines.h:482
double lagrange(const double x[], const double y[], long n, double xval)
#define ASSERT(exp)
Definition: cddefines.h:613
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
double linint(const double x[], const double y[], long n, double xval)
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121