Radau15 Class Reference

#include <orsa_integrator.h>

Inheritance diagram for Radau15:

Inheritance graph
[legend]
Collaboration diagram for Radau15:

Collaboration graph
[legend]
List of all members.

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 *)
Integratorclone () 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

Constructor & Destructor Documentation

Radau15  ) 
 

Definition at line 40 of file orsa_integrator_ra15.cc.

Referenced by Radau15::clone().

00040                    : VariableTimestepIntegrator() {
00041     init();
00042   }

Radau15 const Radau15  ) 
 

Definition at line 44 of file orsa_integrator_ra15.cc.

00044                                     : VariableTimestepIntegrator() {
00045     type     = i.type;
00046     timestep = i.timestep;
00047     accuracy = i.accuracy;
00048     // m        = i.m;
00049     init();
00050   }

~Radau15  ) 
 

Definition at line 137 of file orsa_integrator_ra15.cc.

00137                     {
00138     
00139   }


Member Function Documentation

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.

00225 { return true; }

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   }

Here is the call graph for this function:

IntegratorType GetType  )  const [inline, inherited]
 

Definition at line 98 of file orsa_integrator.h.

References Integrator::type.

00098 { return type; }

void Step const Frame ,
Frame ,
Interaction
[virtual]
 

Implements Integrator.

Definition at line 201 of file orsa_integrator_ra15.cc.

References UniverseTypeAwareTimeStep::GetDouble(), MAX, ORSA_LOGIC_ERROR, 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   }

Here is the call graph for this function:


Member Data Documentation

double accuracy [inherited]
 

Definition at line 91 of file orsa_integrator.h.

Referenced by Evolution::GetIntegratorAccuracy(), OptimizedOrbitPositions::PropagatedOrbit(), orsa::PropagatedSky_J2000(), and Evolution::SetIntegratorAccuracy().

unsigned int m [inherited]
 

used only with variable step size integrators

Definition at line 92 of file orsa_integrator.h.

UniverseTypeAwareTimeStep timestep [inherited]
 

Definition at line 84 of file orsa_integrator.h.

Referenced by Evolution::GetIntegratorTimeStep(), Evolution::Integrate(), OptimizedOrbitPositions::PropagatedOrbit(), orsa::PropagatedSky_J2000(), and Evolution::SetIntegratorTimeStep().

UniverseTypeAwareTimeStep timestep_done [protected, inherited]
 

Definition at line 87 of file orsa_integrator.h.

IntegratorType type [protected, inherited]
 

Definition at line 101 of file orsa_integrator.h.

Referenced by Integrator::GetType().


The documentation for this class was generated from the following files:
Generated on Tue Jan 11 15:29:12 2005 for liborsa by  doxygen 1.4.0