$search
00001 00010 /***************************************************************************** 00011 ** Includes 00012 *****************************************************************************/ 00013 00014 #include <cmath> 00015 #include <ecl/containers/array.hpp> 00016 #include <ecl/exceptions/standard_exception.hpp> 00017 #include "../../include/ecl/geometry/tension_spline.hpp" 00018 00019 /***************************************************************************** 00020 ** Namespaces 00021 *****************************************************************************/ 00022 00023 namespace ecl { 00024 namespace blueprints { 00025 00026 /***************************************************************************** 00027 ** Using 00028 *****************************************************************************/ 00029 00030 using ecl::Array; 00031 using ecl::OutOfRangeError; 00032 using ecl::StandardException; 00033 00034 /***************************************************************************** 00035 ** Implementation [C2CubicSpline] 00036 *****************************************************************************/ 00037 /* 00038 * This algorithm can be found on page 323 of "Numerical Analysis" [Kincaid & 00039 * Cheney]. 00040 */ 00041 C2TensionSpline::C2TensionSpline(const ecl::Array<double>& x_set, const ecl::Array<double>& y_set, 00042 const double &tau) ecl_assert_throw_decl(ecl::StandardException) : 00043 x_data(x_set), 00044 y_data(y_set) 00045 { 00046 ecl_assert_throw( x_data.size() == y_data.size(), StandardException(LOC,OutOfRangeError) ); 00047 00048 /********************* 00049 ** Basic Parameters 00050 **********************/ 00051 unsigned int n = x_data.size()-1; // No. of data points (x_0...x_n) 00052 yddot_data.resize(n+1); 00053 yddot_data[0] = 0; 00054 tension = tau; 00055 00056 /****************************************** 00057 ** Tridiagonal parameters 00058 *******************************************/ 00059 // These equations are derived from the c2 constraints - it results in a tridiagonal system 00060 // of equations. See p.319 of kincaid and cheney for the tridiagonal solution (note this 00061 // is only valid for diagonally dominant tridiagonal matrices! Is this a problem for us? 00062 Array<double> h(n), a(n), beta(n), gamma(n), u(n), v(n); 00063 h[0] = x_set[1] - x_set[0]; 00064 for ( unsigned int i = 0; i < n; ++i ) { 00065 h[i] = x_set[i+1]-x_set[i]; 00066 a[i] = 1/h[i] - tension/sinh(tension*h[i]); 00067 beta[i] = tension*(cosh(tension*h[i])/sinh(tension*h[i])) - 1/h[i]; 00068 gamma[i] = tension*tension*(y_data[i+1]-y_data[i])/h[i]; 00069 } 00070 /* 00071 * [ u1 a1 0 0 0 ] [z1] [ v1 ] 00072 * [ a1 u2 a2 0 0 ] [z2] [ v2 ] 00073 * [ 0 a2 u3 a3 0 ] [z3] [ v3 ] 00074 * 00075 * but with back substitution here, 00076 * 00077 * [ u1 a1 0 0 0 ] [z1] [ v1 ] 00078 * [ 0 (u2-a1^2/u1) a2 0 0 ] [z2] [ v2-a1*v1/u1 ] 00079 * 00080 * etc 00081 */ 00082 u[1] = beta[1] + beta[0]; 00083 v[1] = gamma[1]-gamma[0]; 00084 for ( unsigned int i = 2; i < n; ++i ) { 00085 u[i] = beta[i] + beta[i-1] - a[i-1]*a[i-1]/u[i-1]; 00086 v[i] = gamma[i]-gamma[i-1] - a[i-1]*v[i-1]/u[i-1]; 00087 } 00088 /* 00089 * Now we can roll the upper triangular matrix backwards 00090 */ 00091 yddot_data[n] = 0; 00092 for ( unsigned int i = n-1; i > 0; --i ) { 00093 yddot_data[i] = (v[i] - a[i]*yddot_data[i+1])/u[i]; 00094 } 00095 yddot_data[0] = 0; 00096 } 00097 00098 TensionSpline C2TensionSpline::instantiate() { 00099 TensionSpline cubic; 00100 apply(cubic); 00101 return cubic; 00102 }; 00103 00104 void C2TensionSpline::apply(TensionSpline& spline) const { 00105 00106 spline.discretised_domain = x_data; 00107 spline.tension = tension; 00108 spline.functions.resize(x_data.size()-1); // One less polynomials than there is points 00109 for (unsigned int i = 0; i < spline.functions.size(); ++i ) { 00110 spline.functions[i] = TensionFunction::Interpolation( 00111 x_data[i], y_data[i], yddot_data[i], 00112 x_data[i+1], y_data[i+1], yddot_data[i+1] ); 00113 } 00114 }; 00115 00116 } // namespace blueprints 00117 00118 /***************************************************************************** 00119 ** Implementation [BluePrintFactory][TensionSpline] 00120 *****************************************************************************/ 00121 00122 using blueprints::C2TensionSpline; 00123 00124 00125 C2TensionSpline BluePrintFactory< TensionSpline >::Natural(const Array<double>& x_set, 00126 const Array<double>& y_set, const double &tau) ecl_assert_throw_decl(StandardException) { 00127 return C2TensionSpline(x_set, y_set, tau); 00128 }; 00129 00130 } // namespace ecl