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