AdamsStoermer.hpp

00001 
00002 //
00003 // $Id: AdamsStoermer.hpp,v 1.1 2004/08/12 16:21:33 werner Exp $
00004 //
00005 // $Log: AdamsStoermer.hpp,v $
00006 // Revision 1.1  2004/08/12 16:21:33  werner
00007 // Integration of numerical geodesics now compiles. Working is not yet satisfying.
00008 //
00009 // Revision 1.8  2004/05/12 14:36:22  werner
00010 // Separation of chart definitions on a manifold and the tangential space.
00011 // Introduced convenient coordinate transformations.
00012 //
00013 // Revision 1.7  2004/05/11 18:05:10  werner
00014 //
00015 // Revision 1.6  2004/05/11 16:53:47  werner
00016 // Introduction of coordinate systems for improved type-safety.
00017 // Will support easy coordinate transformations soon.
00018 //
00019 // Revision 1.5  2004/05/06 22:42:16  werner
00020 // Towards a specification of a spacetime via the Acceleration structure.
00021 //
00022 // Revision 1.4  2004/05/05 15:56:52  werner
00023 // Separation of DOP core routines with dynamic size into the ODE library.
00024 //
00025 // Revision 1.3  2004/05/03 13:33:33  werner
00026 // integration improved
00027 //
00028 // Revision 1.2  2004/03/29 11:51:02  werner
00029 // Common interface among simple integrators, DiffMe and Vecal, and preliminiary work on integrating dop853.
00030 //
00031 // Revision 1.1  2004/03/22 11:55:02  werner
00032 // Schwarzschild geodesic integration.
00033 //
00034 // Revision 1.1  2004/02/13 16:36:21  werner
00035 // Initial preliminiary version of the Vector Algebra Library.
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 } /* namespace VecAl */ 
00127 
00128 #endif /* __AdamsStoermer_Geodesic_HPP */