VERB_code_2.2  2
 All Classes Namespaces Files Functions Variables Typedefs Macros Pages
spline.cpp
Go to the documentation of this file.
1 /**
2  * Spline interpolation
3  * \file spline.cpp
4  */
5 
6 //#include <iostream>
7 #include <math.h>
8 #include "../../Exceptions/error.h"
9 
10 //using namespace std;
11 
12 /**
13  * Spline interpolation - calculation of an array with second derivatives (called only once for each interpolation of an array)
14  *
15  * Given arrays x[1..n] and y[1..n] containing a tabulated function, i.e., yi = f (xi ), with
16  * x1 < x2 < . . . < xN , and given values yp1 and ypn for the first derivative of the interpolating
17  * function at points 1 and n, respectively, this routine returns an array y2[1..n] that contains
18  * the second derivatives of the interpolating function at the tabulated points xi . If yp1 and/or
19  * ypn are equal to 1 - 10^30 or larger, the routine is signaled to set the corresponding boundary
20  * condition for a natural spline, with zero second derivative on that boundary.
21  *
22  * \param *x - old x
23  * \param *y - old function
24  * \param n - number of points
25  * \param yp1 - derivative on boundary / flag if > 10^30
26  * \param ypn - derivative on boundary / flag if > 10^30
27  * \param *y2 - array with second derivatives of the function
28  *
29  */
30 void spline(double *x, double *y, int n, double yp1, double ypn, double *y2)
31 {
32  int i,k;
33  double p,qn,sig,un,*u,*vector();
34  void free_vector();
35 
36  u = new double[n];
37  if (yp1 > 0.99e30) {
38  y2[0]=0.0;
39  u[0]=0.0;
40  } else {
41  y2[0] = -0.5;
42  u[0]=(3.0/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-yp1);
43  }
44  for (i=1;i<=n-2;i++) {
45  sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
46  p=sig*y2[i-1]+2.0;
47  y2[i]=(sig-1.0)/p;
48  u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
49  u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
50  }
51  if (ypn > 0.99e30) {
52  qn=0.0;
53  un=0.0;
54  } else {
55  qn=0.5;
56  un=(3.0/(x[n-1]-x[n-2]))*(ypn-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
57  }
58  y2[n-1]=(un-qn*u[n-2])/(qn*y2[n-2]+1.0);
59  for (k=n-2;k>=0;k--) {
60  y2[k]=y2[k]*y2[k+1]+u[k];
61 /* if (!(y2[k] >= 0 || y2[k] <= 0) || y2[k] > 0.99e30) {
62  cout << "error in spline() - inf" << endl;
63  getchar();
64  exit(-1);
65  }*/
66  }
67 
68  delete u;
69 }
70 
71 
72 /**
73  * Spline interpolation
74  *
75  * Given the arrays xa[1..n] and ya[1..n], which tabulate a function (with the xai -s in order),
76  * and given the array y2a[1..n], which is the output from spline above, and given a value of
77  * x, this routine returns a cubic-spline interpolated value y.
78  *
79  * \param *xa - old grid array
80  * \param *ya - old function array
81  * \param *y2a - old function second derivatives array
82  * \param n - size of arrays
83  * \param x - new grid point
84  * \param *y - new grid value
85  */
86 void splint(double *xa, double *ya, double *y2a, int n, double x, double *y)
87 {
88  int klo,khi,k;
89  double h,b,a;
90  void nrerror();
91 
92  klo=0;//1;
93  khi=n-1;//n;
94  while (khi-klo > 1) {
95  k=(khi+klo) >> 1;
96  if (xa[k] > x) khi=k;
97  else klo=k;
98  }
99  h=xa[khi]-xa[klo];
100  if (h == 0.0) //nrerror("Bad XA input to routine SPLINT");
101  throw error_msg("SPLINE_ERROR", "Bad XA input to routine SPLINT");
102  a=(xa[khi]-x)/h;
103  b=(x-xa[klo])/h;
104  *y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
105 /* if (!(*y >= 0 || *y <= 0) || *y > 0.99e30) {
106  cout << "error in splint() - inf" << endl;
107  getchar();
108  exit(-1);
109  }*/
110 
111 }