00001
00006
00007 #include <math.h>
00008 #include "../../Exceptions/error.h"
00009
00010
00011
00030 void spline(double *x, double *y, int n, double yp1, double ypn, double *y2)
00031 {
00032 int i,k;
00033 double p,qn,sig,un,*u,*vector();
00034 void free_vector();
00035
00036 u = new double[n];
00037 if (yp1 > 0.99e30) {
00038 y2[0]=0.0;
00039 u[0]=0.0;
00040 } else {
00041 y2[0] = -0.5;
00042 u[0]=(3.0/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-yp1);
00043 }
00044 for (i=1;i<=n-2;i++) {
00045 sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
00046 p=sig*y2[i-1]+2.0;
00047 y2[i]=(sig-1.0)/p;
00048 u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
00049 u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
00050 }
00051 if (ypn > 0.99e30) {
00052 qn=0.0;
00053 un=0.0;
00054 } else {
00055 qn=0.5;
00056 un=(3.0/(x[n-1]-x[n-2]))*(ypn-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
00057 }
00058 y2[n-1]=(un-qn*u[n-2])/(qn*y2[n-2]+1.0);
00059 for (k=n-2;k>=0;k--) {
00060 y2[k]=y2[k]*y2[k+1]+u[k];
00061
00062
00063
00064
00065
00066 }
00067
00068 delete u;
00069 }
00070
00071
00086 void splint(double *xa, double *ya, double *y2a, int n, double x, double *y)
00087 {
00088 int klo,khi,k;
00089 double h,b,a;
00090 void nrerror();
00091
00092 klo=0;
00093 khi=n-1;
00094 while (khi-klo > 1) {
00095 k=(khi+klo) >> 1;
00096 if (xa[k] > x) khi=k;
00097 else klo=k;
00098 }
00099 h=xa[khi]-xa[klo];
00100 if (h == 0.0)
00101 throw error_msg("SPLINE_ERROR", "Bad XA input to routine SPLINT");
00102 a=(xa[khi]-x)/h;
00103 b=(x-xa[klo])/h;
00104 *y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
00105
00106
00107
00108
00109
00110
00111 }