ADMB Documentation  11.6rc.3396
 All Classes Files Functions Variables Typedefs Friends Defines
cspline.cpp
Go to the documentation of this file.
00001 
00006 #include <fvar.hpp>
00007 
00008 dvector spline(const dvector &x,const dvector&y,double yp1,double ypn);
00009 
00010 double splint(const dvector& xa,const dvector& ya,const dvector& y2a,double x);
00011 
00012 cubic_spline_function::cubic_spline_function(const dvector & _x,
00013   const dvector& _y,double yp1,double ypn) : x(_x) , y(_y)
00014 {
00015   y2.allocate(x);
00016   y2=spline(x,y,yp1,ypn);
00017 }
00018 
00019 double cubic_spline_function::operator () (double u)
00020 {
00021   // need to deal with u<x(1) or u>x(2)
00022   return splint(x,y,y2,u);
00023 }
00024 
00025 dvector cubic_spline_function::operator () (const dvector& u)
00026 {
00027   int mmin=u.indexmin();
00028   int mmax=u.indexmax();
00029   dvector z(mmin,mmax);
00030   for (int i=mmin;i<=mmax;i++)
00031   {
00032     z(i)=splint(x,y,y2,u(i));
00033   }
00034   return z;
00035 }
00036 
00037 dvector spline(const dvector &_x,const dvector&_y,double yp1,double ypn)
00038 {
00039   dvector& x=(dvector&) _x;
00040   dvector& y=(dvector&) _y;
00041   int orig_min=x.indexmin();
00042   x.shift(1);
00043   y.shift(1);
00044   // need to check that x is monotone increasing;
00045   if  ( x.indexmax() != y.indexmax() )
00046   {
00047     cerr << " Incompatible bounds in input to spline" << endl;
00048   }
00049   int n=x.indexmax();
00050   dvector y2(1,n);
00051   int i,k;
00052   double  p,qn,sig,un;
00053 
00054   dvector u(1,n-1);
00055   if (yp1 > 0.99e30)
00056   {
00057     y2[1]=u[1]=0.0;
00058   }
00059   else
00060   {
00061     y2[1] = -0.5;
00062     u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
00063   }
00064   for (i=2;i<=n-1;i++)
00065   {
00066     sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
00067     p=sig*y2[i-1]+2.0;
00068     y2[i]=(sig-1.0)/p;
00069     u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
00070     u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
00071   }
00072   if (ypn > 0.99e30)
00073   {
00074     qn=un=0.0;
00075   }
00076   else
00077   {
00078     qn=0.5;
00079     un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
00080   }
00081   y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
00082   for (k=n-1;k>=1;k--)
00083   {
00084     y2[k]=y2[k]*y2[k+1]+u[k];
00085   }
00086   x.shift(orig_min);
00087   y.shift(orig_min);
00088   y2.shift(orig_min);
00089   return y2;
00090 }
00091 
00092 double splint(const dvector& _xa,const dvector& _ya,const dvector& _y2a,
00093   double x)
00094 {
00095   dvector& xa=(dvector&) _xa;
00096   dvector& ya=(dvector&) _ya;
00097   dvector& y2a=(dvector&) _y2a;
00098   int orig_min=xa.indexmin();
00099   xa.shift(1);
00100   ya.shift(1);
00101   y2a.shift(1);
00102   double y;
00103   int n = xa.indexmax();
00104   int klo,khi,k;
00105   double h,b,a;
00106 
00107   klo=1;
00108   khi=n;
00109   while (khi-klo > 1)
00110   {
00111     k=(khi+klo) >> 1;
00112     if (xa[k] > x) khi=k;
00113     else klo=k;
00114   }
00115   h=xa[khi]-xa[klo];
00116   a=(xa[khi]-x)/h;
00117   b=(x-xa[klo])/h;
00118   y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
00119   xa.shift(orig_min);
00120   ya.shift(orig_min);
00121   y2a.shift(orig_min);
00122   return y;
00123 }