orsa_integrator_stoer.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 namespace orsa {
00031   
00032   Stoer::Stoer() : MultistepIntegrator() {
00033     type = STOER;
00034     m = 8;
00035   }
00036   
00037   Stoer::Stoer(int m_) : MultistepIntegrator() {
00038     type = STOER;
00039     m = m_;
00040   }
00041   
00042   Stoer::Stoer(const Stoer & i) : MultistepIntegrator() {
00043     type     = i.type;
00044     timestep = i.timestep;
00045     accuracy = i.accuracy;
00046     m        = i.m;
00047   }
00048   
00049   Stoer::~Stoer() {
00050     
00051   };
00052   
00053   Integrator * Stoer::clone() const {
00054     return new Stoer(*this);
00055   }
00056   
00057   void Stoer::Step(const Frame &frame_in, Frame &frame_out, Interaction *interaction) {
00058     
00059     const unsigned int n = frame_in.size();
00060     double h, h2, hh; // local_x;
00061     
00062     unsigned int i, k;
00063     
00064     // h  = h_tot / m ;
00065     h  = timestep.GetDouble() / m ;
00066     h2 = h * 0.5   ;
00067     hh = h * h     ;
00068     
00069     std::vector<Vector> temp(n),delta(n);
00070     
00071     frame_out = frame_in;
00072     
00073     // initial values
00074     interaction->Acceleration(frame_out,temp);
00075     
00076     for (i=0;i<frame_out.size();++i) {
00077       delta[i] +=  h * ( frame_out[i].velocity() + h2 * temp[i] );  
00078       frame_out[i].AddToPosition(delta[i]); 
00079     }
00080     
00081     for ( k=1; k <= (m-1) ; ++k) {
00082       if (interaction->IsSkippingJPLPlanets()) {
00083         frame_out.SetTime(frame_in + h*k);
00084         frame_out.ForceJPLEphemerisData();
00085       }
00086       for (i=0;i<frame_out.size();++i) {     
00087         interaction->Acceleration(frame_out,temp);
00088         delta[i] += hh * temp[i];
00089         frame_out[i].AddToPosition(delta[i]);
00090       }
00091     }
00092     
00093     for (i=0;i<frame_out.size();++i) 
00094       frame_out[i].SetVelocity(delta[i]/h + h2 * temp[i]);   
00095     
00096     // frame_out.time += timestep;
00097     // frame_out.SetTime(frame_in.Time() + timestep);
00098     // frame_out.AddTimeStep(timestep);
00099     // frame_out += timestep;
00100     frame_out.SetTime(frame_in + timestep);
00101   }
00102   
00103 } // namespace orsa

Generated on Thu Jul 13 06:45:22 2006 for liborsa by  doxygen 1.4.7