orsa_integrator_runge_kutta.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 <iostream>
00026 
00027 #include "orsa_integrator.h"
00028 #include "orsa_common.h"
00029 
00030 using namespace std;
00031 
00032 namespace orsa {
00033   
00034   RungeKutta::RungeKutta() : FixedTimestepIntegrator() {
00035     type = RUNGEKUTTA;
00036   }
00037   
00038   RungeKutta::RungeKutta(const RungeKutta & i) : FixedTimestepIntegrator() {
00039     type     = i.type;
00040     timestep = i.timestep;
00041     accuracy = i.accuracy;
00042     // m        = i.m;
00043   }
00044   
00045   RungeKutta::~RungeKutta() {
00046     
00047   }
00048   
00049   Integrator * RungeKutta::clone() const {
00050     return new RungeKutta(*this);
00051   }
00052   
00053   void RungeKutta::Step(const Frame & frame_in, Frame & frame_out, Interaction * interaction) {
00054     
00055     // NON-DISSIPATIVE (velocity indipendent) version
00056     
00057     const unsigned int n =  frame_in.size();
00058     
00059     const double h = timestep.GetDouble();
00060     
00061     vector<Vector> acc(n);
00062     
00063     vector<Vector> b1(n), b2(n), b3(n), b4(n);
00064     
00065     unsigned int j;
00066     
00067     interaction->Acceleration(frame_in,acc);
00068     //
00069     b1 = acc;
00070     
00071     Frame tmp_frame = frame_in;
00072     for (j=0;j<tmp_frame.size();++j) {
00073       tmp_frame[j].AddToPosition(frame_in[j].velocity() * h / 2.0);
00074     }
00075     //
00076     if (interaction->IsSkippingJPLPlanets()) {
00077       tmp_frame.SetTime(frame_in + h / 2.0);
00078       tmp_frame.ForceJPLEphemerisData();
00079     }
00080     //
00081     interaction->Acceleration(tmp_frame,acc);
00082     //
00083     b2 = acc;
00084     
00085     tmp_frame = frame_in;
00086     for (j=0;j<tmp_frame.size();++j) {
00087       tmp_frame[j].AddToPosition(frame_in[j].velocity() * h / 2.0 + b1[j] * h * h / 4.0);
00088     }
00089     //
00090     if (interaction->IsSkippingJPLPlanets()) {
00091       tmp_frame.SetTime(frame_in + h / 2.0);
00092       tmp_frame.ForceJPLEphemerisData();
00093     }
00094     //
00095     interaction->Acceleration(tmp_frame,acc);
00096     //
00097     b3 = acc;
00098     
00099     tmp_frame = frame_in;
00100     for (j=0;j<tmp_frame.size();++j) {
00101       tmp_frame[j].AddToPosition(frame_in[j].velocity() * h + b2[j] * h * h / 2.0);
00102     }
00103     //
00104     if (interaction->IsSkippingJPLPlanets()) {
00105       tmp_frame.SetTime(frame_in + timestep);
00106       tmp_frame.ForceJPLEphemerisData();
00107     }
00108     //
00109     interaction->Acceleration(tmp_frame,acc);
00110     //
00111     b4 = acc;
00112     
00113     frame_out = frame_in;
00114     for (j=0;j<frame_out.size();++j) {
00115       frame_out[j].AddToPosition(frame_in[j].velocity() * h + h * h / 6.0 * ( b1[j] + b2[j] + b3[j] ));
00116       frame_out[j].AddToVelocity(h / 6.0 * ( b1[j] + 2 * b2[j] + 2 * b3[j] + b4[j] ));
00117     }
00118     
00119     frame_out += timestep;
00120   }
00121   
00122   // DissipativeRungeKutta
00123   
00124   DissipativeRungeKutta::DissipativeRungeKutta() : FixedTimestepIntegrator() {
00125     type = DISSIPATIVERUNGEKUTTA;
00126   }
00127   
00128   DissipativeRungeKutta::DissipativeRungeKutta(const DissipativeRungeKutta & i) : FixedTimestepIntegrator() {
00129     type     = i.type;
00130     timestep = i.timestep;
00131     accuracy = i.accuracy;
00132     // m        = i.m;
00133   }
00134   
00135   DissipativeRungeKutta::~DissipativeRungeKutta() {
00136     
00137   }
00138   
00139   Integrator * DissipativeRungeKutta::clone() const {
00140     return new DissipativeRungeKutta(*this);
00141   }
00142   
00143   void DissipativeRungeKutta::Step(const Frame & frame_in, Frame & frame_out, Interaction * interaction) {
00144     
00145     // DISSIPATIVE (velocity dependent) version
00146     
00147     const unsigned int n =  frame_in.size();
00148     
00149     const double h = timestep.GetDouble();
00150     
00151     vector<Vector> acc(n);
00152     
00153     vector<Vector> b1(n), b2(n), b3(n), b4(n);
00154     
00155     unsigned int j;
00156     
00157     interaction->Acceleration(frame_in,acc);
00158     //
00159     b1 = acc;
00160     
00161     Frame tmp_frame = frame_in;
00162     for (j=0;j<tmp_frame.size();++j) {
00163       tmp_frame[j].AddToPosition(frame_in[j].velocity() * h / 2.0 + b1[j] * h * h / 8.0);
00164       tmp_frame[j].AddToVelocity(b1[j] * h / 2.0);
00165     }
00166     //
00167     if (interaction->IsSkippingJPLPlanets()) {
00168       tmp_frame.SetTime(frame_in + h / 2.0);
00169       tmp_frame.ForceJPLEphemerisData();
00170     }
00171     //
00172     interaction->Acceleration(tmp_frame,acc);
00173     //
00174     b2 = acc;
00175     
00176     tmp_frame = frame_in;
00177     for (j=0;j<tmp_frame.size();++j) {
00178       tmp_frame[j].AddToPosition(frame_in[j].velocity() * h / 2.0 + b1[j] * h * h / 8.0);
00179       tmp_frame[j].AddToVelocity(b2[j] * h / 2.0);
00180     }
00181     //
00182     if (interaction->IsSkippingJPLPlanets()) {
00183       tmp_frame.SetTime(frame_in + h / 2.0);
00184       tmp_frame.ForceJPLEphemerisData();
00185     }
00186     //
00187     interaction->Acceleration(tmp_frame,acc);
00188     //
00189     b3 = acc;
00190     
00191     tmp_frame = frame_in;
00192     for (j=0;j<tmp_frame.size();++j) {
00193       tmp_frame[j].AddToPosition(frame_in[j].velocity() * h + b3[j] * h * h / 2.0);
00194       tmp_frame[j].AddToVelocity(b3[j] * h);
00195     }
00196     //
00197     if (interaction->IsSkippingJPLPlanets()) {
00198       tmp_frame.SetTime(frame_in + timestep);
00199       tmp_frame.ForceJPLEphemerisData();
00200     }
00201     //
00202     interaction->Acceleration(tmp_frame,acc);
00203     //
00204     b4 = acc;
00205     
00206     frame_out = frame_in;
00207     for (j=0;j<frame_out.size();++j) {
00208       frame_out[j].AddToPosition(frame_in[j].velocity() * h + h * h / 6.0 * ( b1[j] + b2[j] + b3[j] ));
00209       frame_out[j].AddToVelocity(h / 6.0 * ( b1[j] + 2 * b2[j] + 2 * b3[j] + b4[j] ));
00210     }
00211     
00212     frame_out += timestep;
00213   }
00214   
00215 } // namespace orsa

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