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 <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
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
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
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
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
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 }