00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00038 #ifndef __AdamsStoermer_Geodesic_HPP
00039 #define __AdamsStoermer_Geodesic_HPP "Created 27.02.2001 21:42:27 by werner"
00040
00041 #include <vecal/VVector.hpp>
00042 #include <vecal/Matrix.hpp>
00043
00044 #include <ode/Integrator.hpp>
00045 #include <spacetime/RungeKutta.hpp>
00046
00047 namespace Traum
00048 {
00049
00055 template <class Acceleration>
00056 class AdamsStoermerGeodesic : public RungeKuttaGeodesic<Acceleration>
00057 {
00058 public:
00059 typedef typename RungeKuttaGeodesic<Acceleration>::Scalar_t Scalar_t;
00060 typedef typename RungeKuttaGeodesic<Acceleration>::Point_t Point_t;
00061 typedef typename RungeKuttaGeodesic<Acceleration>::Vector_t Vector_t;
00062
00063 using RungeKuttaGeodesic<Acceleration>::x;
00064 using RungeKuttaGeodesic<Acceleration>::v;
00065
00066 private:
00067 int step;
00068 Vector_t b0, b1, b2, b3;
00069
00070 void predict(Point_t&xp, Vector_t&vp, const Scalar_t&ds) const
00071 {
00072 xp = x + v*ds + ds*ds/360.*(323.*b0-264.*b1+159.*b2-38.*b3);
00073 vp = v + ds/24. *( 55.*b0- 59.*b1+ 37.*b2- 9.*b3);
00074 }
00075
00076 void correct(Point_t&xp, Vector_t&vp, const Scalar_t&ds) const
00077 {
00078 Vector_t bp = -Accel(xp, vp);
00079 xp = x + v*ds
00080 + ds*ds/1440.*(135.*bp+752.*b0-246.*b1+ 96.*b2-17.*b3);
00081 vp = v + ds/720. *(251.*bp+646.*b0-264.*b1+106.*b2-19.*b3);
00082 }
00083
00084 void shift()
00085 {
00086 b3 = b2;
00087 b2 = b1;
00088 b1 = b0;
00089 b0 = -Accel(x, v);
00090 }
00091
00092 public:
00093 AdamsStoermerGeodesic (const Acceleration&A)
00094 : RungeKuttaGeodesic<Acceleration>(A)
00095 , step(0)
00096 {}
00097
00098 void reset_step()
00099 {
00100 step = 0;
00101 }
00102
00103 success_code advance(bool bk=false)
00104 {
00105 Scalar_t ds = bk?-this->ds:this->ds;
00106
00107 if (step++ < 4)
00108 {
00109 RungeKuttaGeodesic<Acceleration>::advance(bk);
00110 shift();
00111 return StepOk;
00112 }
00113 Point_t xp;
00114 Vector_t vp;
00115 predict(xp, vp, ds);
00116 correct(xp, vp, ds);
00117 correct(xp, vp, ds);
00118 shift();
00119 x = xp;
00120 v = vp;
00121
00122 return StepOk;
00123 }
00124 };
00125
00126 }
00127
00128 #endif