cloudy  trunk
 All Data Structures 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-2008 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" );
234  cdEXIT(EXIT_FAILURE);
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 );
266  cdEXIT(EXIT_FAILURE);
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 );
304  cdEXIT(EXIT_FAILURE);
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" );
320  cdEXIT(EXIT_FAILURE);
321  }
322  }
323 
324  free(b);
325  free(a);
326  return;
327 }
328 
329 //****************************************************************************80
330 
331 void 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 
433 double 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 
456 double 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 }

Generated for cloudy by doxygen 1.8.1.1