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