#include <orsa_integrator.h>
Public Member Functions | |
Radau15 () | |
Radau15 (const Radau15 &) | |
~Radau15 () | |
bool | can_handle_velocity_dependant_interactions () const |
substeps for multisteps integrators | |
void | Step (const Frame &, Frame &, Interaction *) |
Integrator * | clone () const |
IntegratorType | GetType () const |
Public Attributes | |
UniverseTypeAwareTimeStep | timestep |
double | accuracy |
unsigned int | m |
used only with variable step size integrators | |
Protected Attributes | |
UniverseTypeAwareTimeStep | timestep_done |
IntegratorType | type |
Definition at line 214 of file orsa_integrator.h.
Radau15 | ( | ) |
Definition at line 44 of file orsa_integrator_ra15.cc.
References Integrator::accuracy, Integrator::timestep, and Integrator::type.
00044 : VariableTimestepIntegrator() { 00045 type = i.type; 00046 timestep = i.timestep; 00047 accuracy = i.accuracy; 00048 // m = i.m; 00049 init(); 00050 }
~Radau15 | ( | ) |
bool can_handle_velocity_dependant_interactions | ( | ) | const [inline, virtual] |
substeps for multisteps integrators
Reimplemented from Integrator.
Definition at line 225 of file orsa_integrator.h.
Integrator * clone | ( | ) | const [virtual] |
Implements Integrator.
Definition at line 52 of file orsa_integrator_ra15.cc.
References Radau15::Radau15().
00052 { 00053 return new Radau15(*this); 00054 }
IntegratorType GetType | ( | ) | const [inline, inherited] |
Definition at line 98 of file orsa_integrator.h.
References Integrator::type.
Referenced by OrsaFile::Write().
00098 { return type; }
void Step | ( | const Frame & | frame_in, | |
Frame & | frame_out, | |||
Interaction * | interaction | |||
) | [virtual] |
Implements Integrator.
Definition at line 201 of file orsa_integrator_ra15.cc.
References Interaction::Acceleration(), Integrator::accuracy, Interaction::depends_on_velocity(), Frame::ForceJPLEphemerisData(), UniverseTypeAwareTimeStep::GetDouble(), Interaction::IsSkippingJPLPlanets(), MAX, ORSA_LOGIC_ERROR, UniverseTypeAwareTime::SetTime(), Frame::size(), Integrator::timestep, Integrator::timestep_done, Vector::x, Vector::y, and Vector::z.
00201 { 00202 00203 // cerr << "-> inside Radau15::Step()..." << endl; 00204 00205 // cerr << "RA15: initial timestep: " << timestep << endl; 00206 00207 // N.B. Input/output must be in coordinates with respect to the central body. 00208 00209 // frames... 00210 frame_out = frame_in; 00211 00212 niter = 2; 00213 // if (frame_out.size() != mass.size()) { 00214 if (frame_out.size() != size) { 00215 Bodies_Mass_or_N_Bodies_Changed(frame_out); 00216 niter = 6; 00217 } else { 00218 unsigned int l=0; 00219 while (l != frame_out.size()) { 00220 if (frame_out[l].mass() != mass[l]) { 00221 Bodies_Mass_or_N_Bodies_Changed(frame_out); 00222 niter = 6; 00223 break; 00224 } 00225 ++l; 00226 } 00227 } 00228 00229 // cerr << "niter: " << niter << endl; 00230 00231 interaction->Acceleration(frame_out,acc); 00232 00233 unsigned int j,k; 00234 for(k=0;k<frame_in.size();++k) { 00235 x1[3*k] = frame_in[k].position().x; 00236 x1[3*k+1] = frame_in[k].position().y; 00237 x1[3*k+2] = frame_in[k].position().z; 00238 // 00239 v1[3*k] = frame_in[k].velocity().x; 00240 v1[3*k+1] = frame_in[k].velocity().y; 00241 v1[3*k+2] = frame_in[k].velocity().z; 00242 // 00243 a1[3*k] = acc[k].x; 00244 a1[3*k+1] = acc[k].y; 00245 a1[3*k+2] = acc[k].z; 00246 } 00247 00248 for(k=0;k<nv;++k) { 00249 g[0][k] = b[6][k]*d[15] + b[5][k]*d[10] + b[4][k]*d[6] + b[3][k]*d[3] + b[2][k]*d[1] + b[1][k]*d[0] + b[0][k]; 00250 g[1][k] = b[6][k]*d[16] + b[5][k]*d[11] + b[4][k]*d[7] + b[3][k]*d[4] + b[2][k]*d[2] + b[1][k]; 00251 g[2][k] = b[6][k]*d[17] + b[5][k]*d[12] + b[4][k]*d[8] + b[3][k]*d[5] + b[2][k]; 00252 g[3][k] = b[6][k]*d[18] + b[5][k]*d[13] + b[4][k]*d[9] + b[3][k]; 00253 g[4][k] = b[6][k]*d[19] + b[5][k]*d[14] + b[4][k]; 00254 g[5][k] = b[6][k]*d[20] + b[5][k]; 00255 g[6][k] = b[6][k]; 00256 } 00257 00258 double tmp,gk; 00259 double q1,q2,q3,q4,q5,q6,q7; 00260 00261 unsigned int main_loop_counter; 00262 for (main_loop_counter=0;main_loop_counter<niter;++main_loop_counter) { 00263 for(j=1;j<8;++j) { 00264 00265 // s[0] = timestep * h[j]; 00266 s[0] = timestep.GetDouble() * h[j]; 00267 s[1] = s[0] * s[0] * 0.5; 00268 s[2] = s[1] * h[j] * 0.3333333333333333; 00269 s[3] = s[2] * h[j] * 0.5; 00270 s[4] = s[3] * h[j] * 0.6; 00271 s[5] = s[4] * h[j] * 0.6666666666666667; 00272 s[6] = s[5] * h[j] * 0.7142857142857143; 00273 s[7] = s[6] * h[j] * 0.75; 00274 s[8] = s[7] * h[j] * 0.7777777777777778; 00275 00276 for(k=0;k<nv;++k) { 00277 x[k] = ( s[8]*b[6][k] + 00278 s[7]*b[5][k] + 00279 s[6]*b[4][k] + 00280 s[5]*b[3][k] + 00281 s[4]*b[2][k] + 00282 s[3]*b[1][k] + 00283 s[2]*b[0][k] ) + 00284 s[1]*a1[k] + 00285 s[0]*v1[k] + 00286 x1[k]; 00287 } 00288 00289 // needed only if using a velocity-dependent interaction... 00290 if (interaction->depends_on_velocity()) { 00291 // s[0] = timestep * h[j]; 00292 s[0] = timestep.GetDouble() * h[j]; 00293 s[1] = s[0] * h[j] * 0.5; 00294 s[2] = s[1] * h[j] * 0.6666666666666667; 00295 s[3] = s[2] * h[j] * 0.75; 00296 s[4] = s[3] * h[j] * 0.8; 00297 s[5] = s[4] * h[j] * 0.8333333333333333; 00298 s[6] = s[5] * h[j] * 0.8571428571428571; 00299 s[7] = s[6] * h[j] * 0.875; 00300 00301 for(k=0;k<nv;++k) { 00302 v[k] = ( s[7]*b[6][k] + 00303 s[6]*b[5][k] + 00304 s[5]*b[4][k] + 00305 s[4]*b[3][k] + 00306 s[3]*b[2][k] + 00307 s[2]*b[1][k] + 00308 s[1]*b[0][k] ) + 00309 s[0]*a1[k] + 00310 v1[k]; 00311 } 00312 } 00313 00314 { 00315 Vector rr,vv,drr,dvv; 00316 for(k=0;k<frame_out.size();++k) { 00317 00318 frame_out[k] = frame_in[k]; 00319 00320 rr.x = x[3*k]; 00321 rr.y = x[3*k+1]; 00322 rr.z = x[3*k+2]; 00323 00324 drr = rr - frame_in[k].position(); 00325 frame_out[k].AddToPosition(drr); 00326 00327 vv.x = v[3*k]; 00328 vv.y = v[3*k+1]; 00329 vv.z = v[3*k+2]; 00330 00331 dvv = vv - frame_in[k].velocity(); 00332 frame_out[k].AddToVelocity(dvv); 00333 } 00334 } 00335 00336 if (interaction->IsSkippingJPLPlanets()) { 00337 frame_out.SetTime(frame_in+timestep*h[j]); 00338 frame_out.ForceJPLEphemerisData(); 00339 } 00340 // 00341 interaction->Acceleration(frame_out,acc); 00342 00343 for(k=0;k<frame_out.size();++k) { 00344 a[3*k] = acc[k].x; 00345 a[3*k+1] = acc[k].y; 00346 a[3*k+2] = acc[k].z; 00347 } 00348 00349 switch (j) { 00350 case 1: 00351 for(k=0;k<nv;++k) { 00352 tmp = g[0][k]; 00353 g[0][k] = (a[k] - a1[k]) * r[0]; 00354 b[0][k] += g[0][k] - tmp; 00355 } break; 00356 case 2: 00357 for(k=0;k<nv;++k) { 00358 tmp = g[1][k]; 00359 gk = a[k] - a1[k]; 00360 g[1][k] = (gk*r[1] - g[0][k])*r[2]; 00361 tmp = g[1][k] - tmp; 00362 b[0][k] += tmp * c[0]; 00363 b[1][k] += tmp; 00364 } break; 00365 case 3: 00366 for(k=0;k<nv;++k) { 00367 tmp = g[2][k]; 00368 gk = a[k] - a1[k]; 00369 g[2][k] = ((gk*r[3] - g[0][k])*r[4] - g[1][k])*r[5]; 00370 tmp = g[2][k] - tmp; 00371 b[0][k] += tmp * c[1]; 00372 b[1][k] += tmp * c[2]; 00373 b[2][k] += tmp; 00374 } break; 00375 case 4: 00376 for(k=0;k<nv;++k) { 00377 tmp = g[3][k]; 00378 gk = a[k] - a1[k]; 00379 g[3][k] = (((gk*r[6] - g[0][k])*r[7] - g[1][k])*r[8] - g[2][k])*r[9]; 00380 tmp = g[3][k] - tmp; 00381 b[0][k] += tmp * c[3]; 00382 b[1][k] += tmp * c[4]; 00383 b[2][k] += tmp * c[5]; 00384 b[3][k] += tmp; 00385 } break; 00386 case 5: 00387 for(k=0;k<nv;++k) { 00388 tmp = g[4][k]; 00389 gk = a[k] - a1[k]; 00390 g[4][k] = ((((gk*r[10] - g[0][k])*r[11] - g[1][k])*r[12] - g[2][k])*r[13] - g[3][k])*r[14]; 00391 tmp = g[4][k] - tmp; 00392 b[0][k] += tmp * c[6]; 00393 b[1][k] += tmp * c[7]; 00394 b[2][k] += tmp * c[8]; 00395 b[3][k] += tmp * c[9]; 00396 b[4][k] += tmp; 00397 } break; 00398 case 6: 00399 for(k=0;k<nv;++k) { 00400 tmp = g[5][k]; 00401 gk = a[k] - a1[k]; 00402 g[5][k] = (((((gk*r[15] - g[0][k])*r[16] - g[1][k])*r[17] - g[2][k])*r[18] - g[3][k])*r[19] - g[4][k])*r[20]; 00403 tmp = g[5][k] - tmp; 00404 b[0][k] += tmp * c[10]; 00405 b[1][k] += tmp * c[11]; 00406 b[2][k] += tmp * c[12]; 00407 b[3][k] += tmp * c[13]; 00408 b[4][k] += tmp * c[14]; 00409 b[5][k] += tmp; 00410 } break; 00411 case 7: 00412 for(k=0;k<nv;++k) { 00413 tmp = g[6][k]; 00414 gk = a[k] - a1[k]; 00415 g[6][k] = ((((((gk*r[21] - g[0][k])*r[22] - g[1][k])*r[23] - g[2][k])*r[24] - g[3][k])*r[25] - g[4][k])*r[26] - g[5][k])*r[27]; 00416 tmp = g[6][k] - tmp; 00417 b[0][k] += tmp * c[15]; 00418 b[1][k] += tmp * c[16]; 00419 b[2][k] += tmp * c[17]; 00420 b[3][k] += tmp * c[18]; 00421 b[4][k] += tmp * c[19]; 00422 b[5][k] += tmp * c[20]; 00423 b[6][k] += tmp; 00424 } break; 00425 default: 00426 ORSA_LOGIC_ERROR("aieeee!!!"); 00427 } 00428 00429 } 00430 } 00431 00432 timestep_done = timestep; 00433 00434 // Estimate suitable sequence size for the next call 00435 tmp = 0.0; 00436 for(k=0;k<nv;++k) { 00437 tmp = MAX(tmp,fabs(b[6][k])); 00438 } 00439 // if (tmp!=0.0) tmp /= (72.0 * secure_pow(fabs(timestep),7)); 00440 // if (tmp!=0.0) tmp /= (72.0 * secure_pow(fabs(timestep.GetDouble()),7)); 00441 if (tmp!=0.0) tmp /= (72.0 * pow(fabs(timestep.GetDouble()),7)); 00442 00443 if (tmp < 1.0e-50) { // is equal to zero? 00444 // timestep = timestep_done * 1.4; 00445 timestep = timestep_done * 1.4; 00446 } else { 00447 00448 // old rule... 00449 // timestep = copysign(secure_pow(accuracy/tmp,0.1111111111111111),timestep_done); // 1/9=0.111... 00450 // timestep = copysign(secure_pow(accuracy/tmp,0.1111111111111111),timestep_done.GetDouble()); // 1/9=0.111... 00451 timestep = copysign(pow(accuracy/tmp,0.1111111111111111),timestep_done.GetDouble()); // 1/9=0.111... 00452 00453 } 00454 // 00455 // if (fabs(timestep/timestep_done) < 1.0) { 00456 if (fabs(timestep.GetDouble()/timestep_done.GetDouble()) < 1.0) { 00457 timestep = timestep_done * 0.8; 00458 // std::cerr << "Radau: step rejected! New proposed timestep: " << timestep.GetDouble() << std::endl; 00459 frame_out = frame_in; 00460 niter = 6; 00461 return; 00462 } 00463 // 00464 if (fabs(timestep.GetDouble()/timestep_done.GetDouble()) > 1.4) timestep = timestep_done * 1.4; 00465 00466 // std::cerr << "RA15: new timestep: " << timestep.GetDouble() << std::endl; 00467 00468 // Find new position and velocity values at end of the sequence 00469 tmp = timestep_done.GetDouble() * timestep_done.GetDouble(); 00470 for(k=0;k<nv;++k) { 00471 x1[k] = ( xc[7]*b[6][k] + 00472 xc[6]*b[5][k] + 00473 xc[5]*b[4][k] + 00474 xc[4]*b[3][k] + 00475 xc[3]*b[2][k] + 00476 xc[2]*b[1][k] + 00477 xc[1]*b[0][k] + 00478 xc[0]*a1[k] ) * tmp + v1[k]*timestep_done.GetDouble() + x1[k]; 00479 00480 v1[k] = ( vc[6]*b[6][k] + 00481 vc[5]*b[5][k] + 00482 vc[4]*b[4][k] + 00483 vc[3]*b[3][k] + 00484 vc[2]*b[2][k] + 00485 vc[1]*b[1][k] + 00486 vc[0]*b[0][k] + 00487 a1[k]) * timestep_done.GetDouble() + v1[k]; 00488 } 00489 00490 { 00491 Vector rr,vv,drr,dvv; 00492 for(k=0;k<frame_out.size();++k) { 00493 00494 frame_out[k] = frame_in[k]; 00495 00496 rr.x = x1[3*k]; 00497 rr.y = x1[3*k+1]; 00498 rr.z = x1[3*k+2]; 00499 00500 drr = rr - frame_in[k].position(); 00501 frame_out[k].AddToPosition(drr); 00502 00503 vv.x = v1[3*k]; 00504 vv.y = v1[3*k+1]; 00505 vv.z = v1[3*k+2]; 00506 00507 dvv = vv - frame_in[k].velocity(); 00508 frame_out[k].AddToVelocity(dvv); 00509 } 00510 } 00511 00512 // frame_out += timestep_done; 00513 // frame_out.SetTime(frame_in + timestep_done); 00514 frame_out.SetTime(frame_in); 00515 frame_out += timestep_done; 00516 00517 // Predict new B values to use at the start of the next sequence. The predicted 00518 // values from the last call are saved as E. The correction, BD, between the 00519 // actual and predicted values of B is applied in advance as a correction. 00520 // 00521 q1 = timestep.GetDouble() / timestep_done.GetDouble(); 00522 q2 = q1 * q1; 00523 q3 = q1 * q2; 00524 q4 = q2 * q2; 00525 q5 = q2 * q3; 00526 q6 = q3 * q3; 00527 q7 = q3 * q4; 00528 00529 for(k=0;k<nv;++k) { 00530 00531 s[0] = b[0][k] - e[0][k]; 00532 s[1] = b[1][k] - e[1][k]; 00533 s[2] = b[2][k] - e[2][k]; 00534 s[3] = b[3][k] - e[3][k]; 00535 s[4] = b[4][k] - e[4][k]; 00536 s[5] = b[5][k] - e[5][k]; 00537 s[6] = b[6][k] - e[6][k]; 00538 00539 // Estimate B values for the next sequence 00540 00541 e[0][k] = q1*(b[6][k]* 7.0 + b[5][k]* 6.0 + b[4][k]* 5.0 + b[3][k]* 4.0 + b[2][k]* 3.0 + b[1][k]*2.0 + b[0][k]); 00542 e[1][k] = q2*(b[6][k]*21.0 + b[5][k]*15.0 + b[4][k]*10.0 + b[3][k]* 6.0 + b[2][k]* 3.0 + b[1][k]); 00543 e[2][k] = q3*(b[6][k]*35.0 + b[5][k]*20.0 + b[4][k]*10.0 + b[3][k]* 4.0 + b[2][k]); 00544 e[3][k] = q4*(b[6][k]*35.0 + b[5][k]*15.0 + b[4][k]* 5.0 + b[3][k]); 00545 e[4][k] = q5*(b[6][k]*21.0 + b[5][k]* 6.0 + b[4][k]); 00546 e[5][k] = q6*(b[6][k]* 7.0 + b[5][k]); 00547 e[6][k] = q7* b[6][k]; 00548 00549 b[0][k] = e[0][k] + s[0]; 00550 b[1][k] = e[1][k] + s[1]; 00551 b[2][k] = e[2][k] + s[2]; 00552 b[3][k] = e[3][k] + s[3]; 00553 b[4][k] = e[4][k] + s[4]; 00554 b[5][k] = e[5][k] + s[5]; 00555 b[6][k] = e[6][k] + s[6]; 00556 00557 } 00558 00559 // cerr << "-> out of Radau15::Step()..." << endl; 00560 00561 }
double accuracy [inherited] |
Definition at line 91 of file orsa_integrator.h.
Referenced by DissipativeRungeKutta::DissipativeRungeKutta(), Evolution::GetIntegratorAccuracy(), Leapfrog::Leapfrog(), OptimizedOrbitPositions::PropagatedOrbit(), orsa::PropagatedSky_J2000(), Radau15::Radau15(), RungeKutta::RungeKutta(), Evolution::SetIntegratorAccuracy(), Radau15::Step(), Stoer::Stoer(), and OrsaFile::Write().
unsigned int m [inherited] |
used only with variable step size integrators
Definition at line 92 of file orsa_integrator.h.
Referenced by Stoer::Step(), Stoer::Stoer(), and OrsaFile::Write().
UniverseTypeAwareTimeStep timestep [inherited] |
Definition at line 84 of file orsa_integrator.h.
Referenced by DissipativeRungeKutta::DissipativeRungeKutta(), Evolution::GetIntegratorTimeStep(), Evolution::Integrate(), Leapfrog::Leapfrog(), OptimizedOrbitPositions::PropagatedOrbit(), orsa::PropagatedSky_J2000(), Radau15::Radau15(), RungeKutta::RungeKutta(), Evolution::SetIntegratorTimeStep(), orsa::StartFrame(), Stoer::Step(), DissipativeRungeKutta::Step(), RungeKutta::Step(), Radau15::Step(), Leapfrog::Step(), Stoer::Stoer(), and OrsaFile::Write().
UniverseTypeAwareTimeStep timestep_done [protected, inherited] |
IntegratorType type [protected, inherited] |
Definition at line 101 of file orsa_integrator.h.
Referenced by DissipativeRungeKutta::DissipativeRungeKutta(), Integrator::GetType(), Leapfrog::Leapfrog(), Radau15::Radau15(), RungeKutta::RungeKutta(), and Stoer::Stoer().