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,n1); 

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<=n1;i++) 

115 
{ 

116 
sig=(x[i]x[i1])/(x[i+1]x[i1]); 

117 
p=sig*y2[i1]+2.0; 

118 
y2[i]=(sig1.0)/p; 

119 
u[i]=(y[i+1]y[i])/(x[i+1]x[i])  (y[i]y[i1])/(x[i]x[i1]); 

120 
u[i]=(6.0*u[i]/(x[i+1]x[i1])sig*u[i1])/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[n1]))*(ypn(y[n]y[n1])/(x[n]x[n1])); 

130 
} 

131 
y2[n]=(unqn*u[n1])/(qn*y2[n1]+1.0); 

132 
for (k=n1;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 
* 
Also available in: Unified diff