00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
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
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
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
00055 }
00056
00057 Integrator * Leapfrog::clone() const {
00058 return new Leapfrog(*this);
00059 }
00060
00061 Leapfrog::~Leapfrog() {
00062
00063 }
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107 void Leapfrog::Step(const Frame & frame_in, Frame & frame_out, Interaction * interaction) {
00108
00109
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 }