00001 #include "TrackPropagation/RungeKutta/interface/RK4OneStep.h" 00002 #include "TrackPropagation/RungeKutta/interface/RKCartesianDerivative.h" 00003 00004 #include <iostream> 00005 00006 CartesianState 00007 RK4OneStep::operator()( const CartesianState& start, const RKCartesianDerivative& deriv, 00008 double step) const 00009 { 00010 double s0 = 0; // derivatives do't depend on absolute value of the integration variable 00011 CartesianState k1 = step * deriv( s0, start); 00012 CartesianState k2 = step * deriv( s0+step/2, start+k1/2); 00013 CartesianState k3 = step * deriv( s0+step/2, start+k2/2); 00014 CartesianState k4 = step * deriv( s0+step, start+k3); 00015 /* 00016 std::cout << "k1 = " << k1.position() << k1.momentum() << std::endl; 00017 std::cout << "k2 = " << k2.position() << k2.momentum() << std::endl; 00018 std::cout << "k3 = " << k3.position() << k3.momentum() << std::endl; 00019 std::cout << "k4 = " << k4.position() << k4.momentum() << std::endl; 00020 */ 00021 CartesianState result = start + k1/6 + k2/3 + k3/3 + k4/6; 00022 return result; 00023 }