VERB_code_2.3
spline.cpp
Go to the documentation of this file.
1 
6 //#include <iostream>
7 #include <math.h>
8 #include "../../Exceptions/error.h"
9 
10 //using namespace std;
11 
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 
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 }
void splint(double *xa, double *ya, double *y2a, int n, double x, double *y)
Definition: spline.cpp:86
void spline(double *x, double *y, int n, double yp1, double ypn, double *y2)
Definition: spline.cpp:30
Error message - stack of single_errors.
Definition: error.h:54