Revision 795 branches/replacement/src/linad99/cspline.cpp

cspline.cpp (revision 795)
4 4
 * Author: David Fournier
5 5
 * Copyright (c) 2009 ADMB Foundation
6 6
 */
7
/**
8
 * \file
9
 * Contains routines for cubic spline interpolation
10
 * for constant types.
11
 */
12
/**
13
 * \defgroup cub_spline
14
 */
15

  
7 16
#include <fvar.hpp>
8 17

  
9 18
dvector spline(BOR_CONST dvector &x,BOR_CONST dvector&y,double yp1,double ypn);
10
dvector spline_cubic_set (int n,_CONST dvector& t,_CONST dvector& y, int ibcbeg,
19
dvector spline_cubic_set(int n,_CONST dvector& t,_CONST dvector& y, int ibcbeg,
11 20
  double ybcbeg, int ibcend, double ybcend );
12 21
double spline_cubic_val(int n, _CONST dvector& t, double tval,
13 22
  _CONST dvector& y, _CONST dvector& ypp);
......
50 59
 *        end point
51 60
 * \return an array containing the second derivatives
52 61
*/
53
dvector spline(BOR_CONST dvector &_x,BOR_CONST dvector&_y,double yp1,double ypn)
62
dvector spline(BOR_CONST dvector& _x,BOR_CONST dvector& _y,double yp1,double ypn)
54 63
{
55 64
  int ibcbeg, ibcend;
56 65
  double ybcbeg, ybcend;
......
83 92
  return ret;
84 93
}
85 94

  
86

  
87
/*dvector spline(BOR_CONST dvector &_x,BOR_CONST dvector&_y,double yp1,double ypn)
88
{
89
  dvector& x=(dvector&) _x;
90
  dvector& y=(dvector&) _y;
91
  int orig_min=x.indexmin();
92
  x.shift(1);
93
  y.shift(1);
94
  // need to check that x is monotone increasing;
95
  if  ( x.indexmax() != y.indexmax() )
96
  {
97
    cerr << " Incompatible bounds in input to spline" << endl;
98
  }
99
  int n=x.indexmax();
100
  dvector y2(1,n);
101
  int i,k;
102
  double  p,qn,sig,un;
103

  
104
  dvector u(1,n-1);
105
  if (yp1 > 0.99e30)
106
  {
107
    y2[1]=u[1]=0.0;
108
  }
109
  else
110
  {
111
    y2[1] = -0.5;
112
    u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
113
  }
114
  for (i=2;i<=n-1;i++)
115
  {
116
    sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
117
    p=sig*y2[i-1]+2.0;
118
    y2[i]=(sig-1.0)/p;
119
    u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
120
    u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
121
  }
122
  if (ypn > 0.99e30)
123
  {
124
    qn=un=0.0;
125
  }
126
  else 
127
  {
128
    qn=0.5;
129
    un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
130
  }
131
  y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
132
  for (k=n-1;k>=1;k--)
133
  {
134
    y2[k]=y2[k]*y2[k+1]+u[k];
135
  }
136
  x.shift(orig_min);
137
  y.shift(orig_min);
138
  y2.shift(orig_min);
139
  return y2;
140
}*/
141

  
142

  
143

  
144

  
145
/** \ingroup cub_spline
95
/**
146 96
 *  Cubic spline interpolation.
147
 *
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff