orsa_integrator.cc

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 #include "orsa_integrator.h"
00026 
00027 namespace orsa {
00028   
00029   void make_new_integrator(Integrator ** i, const IntegratorType type) {
00030     
00031     delete (*i);
00032     (*i) = 0;
00033     
00034     switch (type) {
00035     case STOER:                 (*i) = new Stoer;                  break;
00036       // case BULIRSCHSTOER:         (*i) = new BulirschStoer;          break;
00037     case RUNGEKUTTA:            (*i) = new RungeKutta;             break;
00038     case DISSIPATIVERUNGEKUTTA: (*i) = new DissipativeRungeKutta;  break;
00039     case RA15:                  (*i) = new Radau15;                break;
00040     case LEAPFROG:              (*i) = new Leapfrog;               break;
00041     }
00042   }
00043   
00044   // Leapfrog
00045   
00046   Leapfrog::Leapfrog() : FixedTimestepIntegrator() {
00047     type = LEAPFROG;
00048   }
00049   
00050   Leapfrog::Leapfrog(const Leapfrog & i) : FixedTimestepIntegrator() {
00051     type     = i.type;
00052     timestep = i.timestep;
00053     accuracy = i.accuracy;
00054     // m        = i.m;
00055   }
00056   
00057   Integrator * Leapfrog::clone() const {
00058     return new Leapfrog(*this);
00059   }
00060   
00061   Leapfrog::~Leapfrog() {
00062     
00063   }
00064   
00065   /* 
00066      // OLD method: kick-drift-kick
00067      void Leapfrog::Step(const Frame & frame_in, Frame & frame_out, Interaction * interaction) {
00068      
00069      // NON-DISSIPATIVE (velocity indipendent) version
00070      
00071      const unsigned int n = frame_in.size();
00072      
00073      const double h  = timestep.GetDouble();
00074      const double h2 = 0.5*h;
00075      
00076      acc.resize(n);
00077      
00078      if (acc_time != frame_in) {
00079      interaction->Acceleration(frame_in,acc);
00080      acc_time = frame_in;
00081      }
00082      
00083      std::vector<Vector> vel_h2(n);
00084      
00085      for (unsigned int j=0;j<n;++j) {
00086      vel_h2[j] = frame_in[j].velocity()+acc[j]*h2;
00087      }
00088      
00089      frame_out = frame_in;
00090      frame_out += timestep;
00091      
00092      for (unsigned int j=0;j<n;++j) {
00093      frame_out[j].AddToPosition(vel_h2[j]*h);
00094      }
00095      
00096      interaction->Acceleration(frame_out,acc);
00097      acc_time = frame_out;
00098      
00099      for (unsigned int j=0;j<n;++j) {
00100      frame_out[j].SetVelocity(vel_h2[j]+acc[j]*h2);
00101      }
00102      
00103      }
00104   */
00105   
00106   // new method: drift-kick-drift
00107   void Leapfrog::Step(const Frame & frame_in, Frame & frame_out, Interaction * interaction) {
00108     
00109     // NON-DISSIPATIVE (velocity indipendent) version
00110     
00111     const unsigned int n = frame_in.size();
00112     
00113     const double h  = timestep.GetDouble();
00114     const double h2 = 0.5*h;
00115     
00116     frame_out = frame_in;
00117     
00118     frame_out += 0.5*timestep;
00119     for (unsigned int j=0;j<n;++j) {
00120       frame_out[j].AddToPosition(frame_out[j].velocity()*h2);
00121     }
00122     
00123     std::vector<Vector> acc(n);
00124     //
00125     if (interaction->IsSkippingJPLPlanets()) {
00126       frame_out.ForceJPLEphemerisData();
00127     }
00128     //
00129     interaction->Acceleration(frame_out,acc);
00130     
00131     for (unsigned int j=0;j<n;++j) {
00132       frame_out[j].AddToVelocity(acc[j]*h);
00133     }
00134     
00135     frame_out += 0.5*timestep;
00136     for (unsigned int j=0;j<n;++j) {
00137       frame_out[j].AddToPosition(frame_out[j].velocity()*h2);
00138     }
00139     
00140   }
00141   
00142 } // namespace orsa

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