orsa_integrator.h

Go to the documentation of this file.
00001 /* 
00002    ORSA - Orbit Reconstruction, Simulation and Analysis
00003    Copyright (C) 2002-2004 Pasquale Tricarico
00004    
00005    This program is free software; you can redistribute it and/or
00006    modify it under the terms of the GNU General Public License
00007    as published by the Free Software Foundation; either version 2
00008    of the License, or (at your option) any later version.
00009    
00010    As a special exception, Pasquale Tricarico gives permission to
00011    link this program with Qt commercial edition, and distribute the
00012    resulting executable, without including the source code for the Qt
00013    commercial edition in the source distribution.
00014    
00015    This program is distributed in the hope that it will be useful,
00016    but WITHOUT ANY WARRANTY; without even the implied warranty of
00017    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00018    GNU General Public License for more details.
00019    
00020    You should have received a copy of the GNU General Public License
00021    along with this program; if not, write to the Free Software
00022    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
00023 */
00024 
00025 #ifndef _ORSA_INTEGRATOR_H_
00026 #define _ORSA_INTEGRATOR_H_
00027 
00028 #include "orsa_frame.h"
00029 #include "orsa_error.h"
00030 #include "orsa_interaction.h"
00031 
00032 #include <string>
00033 
00034 namespace orsa {
00035   
00036   enum IntegratorType {
00037     STOER=1,
00038     BULIRSCHSTOER=2,
00039     RUNGEKUTTA=3,
00040     DISSIPATIVERUNGEKUTTA=4,
00041     RA15=5,
00042     LEAPFROG=6
00043   };
00044   
00045   inline void convert(IntegratorType &it, const unsigned int i)  {
00046     switch(i) {
00047     case 1: it = STOER;                 break;
00048     case 2: it = BULIRSCHSTOER;         break;
00049     case 3: it = RUNGEKUTTA;            break;
00050     case 4: it = DISSIPATIVERUNGEKUTTA; break;
00051     case 5: it = RA15;                  break;
00052     case 6: it = LEAPFROG;              break;
00053       //
00054     default:
00055       ORSA_ERROR("conversion problem: i = %i",i);
00056       break;       
00057     }
00058   }
00059   
00060   inline std::string label(const IntegratorType it) {
00061     std::string s;
00062     switch (it) {
00063     case STOER:                  s="Stoer";                              break;
00064     case BULIRSCHSTOER:          s="Bulirsch-Stoer";                     break;
00065     case RUNGEKUTTA:             s="Runge-Kutta 4th order";              break;
00066     case DISSIPATIVERUNGEKUTTA:  s="Dissipative Runge-Kutta 4th order";  break;
00067     case RA15:                   s="Everhart's RADAU 15th order";        break;
00068     case LEAPFROG:               s="Leapfrog 2nd order";                 break;
00069     }   
00070     return s;
00071   }
00072   
00073   //! This is the interface for all the Integrator classes.
00074   class Integrator {
00075   public:
00076     virtual void Step(const Frame &, Frame &, Interaction *) = 0;
00077     virtual ~Integrator() { };
00078     
00079   public:
00080     virtual Integrator * clone() const = 0;
00081     
00082   public:
00083     // double timestep,timestep_done;
00084     UniverseTypeAwareTimeStep timestep;
00085     
00086   protected:
00087     UniverseTypeAwareTimeStep timestep_done;
00088     
00089   public:
00090     // these vars should be moved inside each single class
00091     double accuracy; //! used only with variable step size integrators
00092     unsigned int m;  //! substeps for multisteps integrators
00093     
00094   public:
00095     virtual bool can_handle_velocity_dependant_interactions() const { return false; }
00096     
00097   public:
00098     inline IntegratorType GetType() const { return type; }
00099     
00100   protected: 
00101     IntegratorType type;
00102   };
00103   
00104   void make_new_integrator(Integrator**, const IntegratorType);
00105   
00106   // 
00107   // First classification: fixed timestep and variable timestep
00108   //
00109   
00110   class FixedTimestepIntegrator : public Integrator {
00111     
00112   };
00113   
00114   class MultistepIntegrator : public FixedTimestepIntegrator {
00115     
00116   };
00117   
00118   class VariableTimestepIntegrator : public Integrator {
00119     
00120   };
00121   
00122   //
00123   // derived classes
00124   //
00125   
00126   //! Advances using a number of substeps (midpoints). For conservative and non-conservative Interaction.
00127   class ModifiedMidpoint : public MultistepIntegrator {
00128     
00129   };
00130   
00131   //! Like the ModifiedMidpoint, but for conservative Interaction.
00132   class Stoer : public MultistepIntegrator {
00133   public:
00134     Stoer();
00135     Stoer(int);
00136     Stoer(const Stoer &);
00137     ~Stoer();
00138     
00139   public:
00140     void Step(const Frame&, Frame&, Interaction*);
00141     
00142   public:
00143     Integrator * clone() const;
00144   };
00145   
00146   /* 
00147      class BulirschStoer : public VariableTimestepIntegrator {
00148      public:
00149      BulirschStoer();
00150      // BulirschStoer(Interaction*);
00151      
00152      ~BulirschStoer();
00153      
00154      private:
00155      void init();
00156      
00157      public:
00158      void Step(const Frame&, Frame&, Interaction*);
00159      
00160      
00161      // for internal computations only
00162      private:
00163      double h, h_done, h_next, h_try;
00164      Frame main_frame, frame_save, frame_seq, frame_out_stoerm;
00165      std::vector<unsigned int> midpoint_sequence; // {2,4,6,8,10,12,14,16,18,20};
00166      const unsigned int max_counter_midpoint; // = 8; // ==> midpoint_sequence[]={2,4,...,14,16}
00167      const unsigned int max_n_points; //= max_counter_midpoint;  // to delete in future
00168      Stoer *stoer;
00169      //
00170      bool flag_first;
00171      unsigned int k_optimum, k_max;
00172      double old_accuracy;
00173      std::vector<double> a;
00174      std::vector< std::vector <double> > alpha;
00175      double x_new;
00176      };
00177   */
00178   
00179   class RungeKutta : public FixedTimestepIntegrator {
00180    public:
00181     RungeKutta();
00182     RungeKutta(const RungeKutta &);
00183     ~RungeKutta();
00184     
00185   public:
00186     void Step(const Frame&, Frame&, Interaction*);
00187     
00188   public:
00189     Integrator * clone() const;
00190   };
00191   
00192   class DissipativeRungeKutta : public FixedTimestepIntegrator {
00193   public:
00194     DissipativeRungeKutta();
00195     DissipativeRungeKutta(const DissipativeRungeKutta &);
00196     ~DissipativeRungeKutta();
00197     
00198   public:
00199     bool can_handle_velocity_dependant_interactions() const { return true; }
00200     
00201   public:
00202     void Step(const Frame&, Frame&, Interaction*);
00203     
00204   public:
00205     Integrator * clone() const;
00206   };
00207   
00208   /* 
00209      class Symplectic : public FixedTimestepIntegrator {
00210      
00211      };
00212   */
00213   
00214   class Radau15 : public VariableTimestepIntegrator {
00215   public:
00216     Radau15();
00217     Radau15(const Radau15 &);
00218     ~Radau15();
00219     
00220   private:
00221     void init();
00222     void Bodies_Mass_or_N_Bodies_Changed(const Frame&);
00223     
00224   public:
00225     bool can_handle_velocity_dependant_interactions() const { return true; }
00226     
00227   public:
00228     void Step(const Frame&, Frame&, Interaction*);
00229     
00230   public:
00231     Integrator * clone() const;
00232     
00233   private:
00234     // Gauss-Radau spacings for substeps within a sequence, for the 15th order 
00235     // integrator. The sum of the H values should be 3.733333333333333 
00236     // const vector<double> h;
00237     // vector<double> h;
00238     double h[8];
00239     
00240     // Constant coefficients used in series expansions for X and V
00241     //  XC: 1/2,  1/6,  1/12, 1/20, 1/30, 1/42, 1/56, 1/72
00242     //  VC: 1/2,  1/3,  1/4,  1/5,  1/6,  1/7,  1/8
00243     // vector<double> xc,vc;
00244     double xc[8],vc[7];
00245     
00246     // vector<double> r,c,d,s;
00247     double r[28],c[21],d[21],s[9];
00248     
00249     std::vector< std::vector<double> > g,b,e;
00250     // double g[7][3000],b[7][3000],e[7][3000];
00251     
00252     unsigned int nv,niter;
00253     
00254     std::vector<double> x,v,a,x1,v1,a1;
00255     // double x[3000],v[3000],a[3000],x1[3000],v1[3000],a1[3000];
00256     
00257     std::vector<double> mass;
00258     // double mass[1000];
00259     
00260     std::vector<Vector> acc;
00261     // Vector acc[1000];
00262     
00263     unsigned int size;
00264   };  
00265   
00266   class Leapfrog : public FixedTimestepIntegrator {
00267   public:
00268     Leapfrog(); 
00269     Leapfrog(const Leapfrog &);
00270     ~Leapfrog();
00271     
00272   public:
00273     void Step(const Frame&, Frame&, Interaction*);
00274     
00275   public:
00276     Integrator * clone() const;
00277   };
00278   
00279 } // namespace orsa
00280 
00281 #endif // _ORSA_INTEGRATOR_H_
00282 

Generated on Fri Nov 3 20:37:42 2006 for liborsa by  doxygen 1.4.7