orsa_interaction.cc

Go to the documentation of this file.
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_interaction.h"
00026 #include "orsa_secure_math.h"
00027 #include "orsa_universe.h"
00028 #include "orsa_error.h"
00029 
00030 #include <cmath>
00031 #include <iostream>
00032 
00033 using namespace std;
00034 
00035 namespace orsa {
00036   
00037   std::string label(const InteractionType it) {
00038     
00039     std::string s("");
00040     
00041     switch (it) {
00042     case NEWTON:                                s = "Newton";                               break;
00043     case NEWTON_MPI:                            s = "Newton (MPI)";                         break;
00044     case ARMONICOSCILLATOR:                     s = "Armonic Oscillator";                   break;
00045     case GALACTIC_POTENTIAL_ALLEN:              s = "Galactic Potential (Allen)";           break; 
00046     case GALACTIC_POTENTIAL_ALLEN_PLUS_NEWTON:  s = "Galactic Potential (Allen) + Newton";  break;
00047     case JPL_PLANETS_NEWTON:                    s = "JPL planets + Newton";                 break; 
00048     case GRAVITATIONALTREE:                     s = "Gravitational TreeCode";               break;
00049     case RELATIVISTIC:                          s = "Newton + Relativistic effects";        break;
00050     }   
00051     
00052     return s;
00053   }
00054   
00055   void make_new_interaction(Interaction ** i, const InteractionType type) {
00056     
00057     delete (*i);
00058     (*i) = 0;
00059     
00060     switch (type) {
00061     case NEWTON:                               (*i) = new Newton;                           break;
00062     case NEWTON_MPI:                        
00063 #ifdef HAVE_MPI
00064       (*i) = new Newton_MPI;                  
00065 #else
00066       ORSA_WARNING("read NEWTON_MPI interaction from application without MPI support.");
00067 #endif
00068       break;
00069     case ARMONICOSCILLATOR:                    (*i) = new ArmonicOscillator(1,1);           break;
00070     case GALACTIC_POTENTIAL_ALLEN:             (*i) = new GalacticPotentialAllen;           break;
00071     case GALACTIC_POTENTIAL_ALLEN_PLUS_NEWTON: (*i) = new GalacticPotentialAllenPlusNewton; break;
00072     case JPL_PLANETS_NEWTON:                   /*************************************/      break;
00073     case GRAVITATIONALTREE:                    (*i) = new GravitationalTree;                break;
00074     case RELATIVISTIC:                         (*i) = new Relativistic;                     break;
00075     }    
00076   }
00077   
00078   // MappedTable 
00079   
00080   void MappedTable::load(const std::vector<Body> & f, const bool skip_JPL_planets) {
00081     const unsigned int f_size = f.size();
00082     N = f_size;
00083     mapping.resize(N);
00084     M = 0;
00085     for (unsigned int k=0; k<f_size; ++k) {
00086       mapping[k] = k; 
00087       if (!f[k].has_zero_mass()) {
00088         mapping[k] = mapping[M];
00089         mapping[M] = k;
00090         ++M;
00091       }
00092     }
00093     //
00094     MN = M*N;
00095     //
00096     if (MN != distance_vector.size()) {
00097       distance_vector.resize(MN);
00098       d1.resize(MN);
00099       d2.resize(MN);
00100       d3.resize(MN);
00101       d4.resize(MN);
00102       one_over_distance.resize(MN);
00103       one_over_distance_square.resize(MN);
00104       one_over_distance_cube.resize(MN);
00105       distance_vector_over_distance_cube.resize(MN);
00106     }
00107     //
00108     Vector d;
00109     Vector v;
00110     double l;
00111     double one_over_d;
00112     double one_over_d2;
00113     double one_over_d3;
00114     unsigned int index;
00115     //
00116     for (unsigned int i=0; i<(N-1); ++i) {
00117       for (unsigned int j=i+1; j<N; ++j) {
00118         // if (i == j) continue;
00119         if (f[i].has_zero_mass() && f[j].has_zero_mass()) continue;
00120         if (skip_JPL_planets && (f[i].JPLPlanet() != NONE) && (f[i].JPLPlanet() != NONE)) continue;
00121         //
00122         index = ij_to_index(i,j);
00123         //
00124         // unsing: inline Vector Body::distanceVector(const Body & b) const { return b.position()-position(); }
00125         // so: d = f[j].position() - f[i].position();
00126         d           = f[i].DistanceVector(f[j]);
00127         //
00128         if (distance_vector[index] == d) continue; // everything already computed!
00129         //
00130         l           = d.Length();
00131         one_over_d  = 1.0/l;
00132         one_over_d2 = one_over_d*one_over_d;
00133         one_over_d3 = one_over_d*one_over_d2;
00134         //
00135         // index = ij_to_index(i,j);
00136         //
00137         //
00138         distance_vector[index]        = d;
00139         d1[index]                     = l;
00140         d2[index]                     = l*l;
00141         d3[index]                     = d2[index]*l;
00142         d4[index]                     = d3[index]*l;
00143         one_over_distance[index]        = one_over_d;
00144         one_over_distance_square[index] = one_over_d2;
00145         one_over_distance_cube[index]   = one_over_d3;
00146         distance_vector_over_distance_cube[index] = d*one_over_d3;
00147       }
00148     }
00149   }
00150   
00151   // Newton
00152   
00153   Newton::Newton() : Interaction(), include_multipole_moments(false), include_relativistic_effects(false), include_fast_relativistic_effects(false), one_over_c2(1.0/(GetC()*GetC())) {
00154     skip_JPL_planets = false;
00155   }
00156   
00157   Newton::Newton(const Newton & n) : Interaction(), include_multipole_moments(n.include_multipole_moments), include_relativistic_effects(n.include_relativistic_effects), include_fast_relativistic_effects(n.include_fast_relativistic_effects), one_over_c2(1.0/(GetC()*GetC())) {
00158     skip_JPL_planets = n.skip_JPL_planets;
00159   }
00160   
00161   Interaction * Newton::clone() const {
00162     return new Newton(*this);
00163   }
00164   
00165   void Newton::fast_newton_acc(const Frame & f, std::vector<Vector> & a) {
00166     
00167     const unsigned int size = f.size();
00168     
00169     Vector d;
00170     
00171     double l;
00172     
00173     for (unsigned int i=0;i<size-1;++i) {
00174       for (unsigned int j=i+1;j<size;++j) {
00175         
00176         if (zero_mass[i] && zero_mass[j]) continue;
00177 
00178         if (skip[i] && skip[j]) continue;
00179         
00180         d = f[i].DistanceVector(f[j]);
00181         
00182         l = d.Length();
00183         
00184         if (d.IsZero()) {
00185           ORSA_WARNING("two objects in the same position! (%s and %s)",f[i].name().c_str(),f[j].name().c_str());
00186           continue;
00187         }
00188         
00189         d /= l*l*l;
00190         
00191         a[i] += d * f[j].mu();
00192         a[j] -= d * f[i].mu();
00193         
00194         // fprintf(stderr,"body %s: acc. due to %s = %g    force = %g\n",f[i].name().c_str(),f[j].name().c_str(),a[i].Length(),a[i].Length()*f[i].mass());
00195       } 
00196     } 
00197   }
00198   
00199   void Newton::Acceleration(const Frame & f, vector<Vector> & a) {
00200     
00201     const unsigned int size = f.size();
00202     
00203     if (size < 2) return;
00204     
00205     a.resize(size);
00206     /* 
00207        for (unsigned int i=0;i<size;++i)
00208        a[i].Set(0.0,0.0,0.0);
00209     */
00210     
00211     zero_mass.resize(size);
00212     for (unsigned int i=0;i<size;++i) {
00213       zero_mass[i] = f[i].has_zero_mass();
00214     }
00215     
00216     skip.resize(size);
00217     if (skip_JPL_planets) {
00218       for (unsigned int k=0;k<size;++k) {
00219         skip[k] = (f[k].JPLPlanet() != NONE);
00220       }
00221     } else {
00222       for (unsigned int k=0;k<size;++k) {
00223         skip[k] = false;
00224       }
00225     }
00226     
00227     const bool only_newton = !(include_multipole_moments || include_relativistic_effects || include_fast_relativistic_effects);
00228     
00229     if (only_newton) {
00230       for (unsigned int i=0;i<size;++i) {
00231         a[i].Set(0.0,0.0,0.0);
00232       }
00233       //
00234       fast_newton_acc(f,a);
00235       //
00236       return;
00237     }
00238     
00239     if (include_relativistic_effects && include_fast_relativistic_effects) {
00240       ORSA_WARNING("Both the accurate and the fast version of the Relativistic corrections are activated!");
00241     }
00242     
00243     if (a_newton.size() != size) {
00244       a_newton.resize(size);
00245       a_multipoles.resize(size);
00246       a_relativity.resize(size);
00247     }
00248     //
00249     for (unsigned int i=0;i<size;++i) {
00250       a_newton[i].Set(0.0,0.0,0.0);
00251       a_multipoles[i].Set(0.0,0.0,0.0);
00252       a_relativity[i].Set(0.0,0.0,0.0);
00253     }
00254     
00255     // mod_f initialization needed!
00256     Frame mod_f = f;
00257     //
00258     if (include_relativistic_effects) {
00259       // a_newton.resize(size);
00260       //
00261       // The RelativisticBarycenter is _very_ close to the normal Barycenter,
00262       // and the accuracy of a "double" is often barely enough to distinguish them.
00263       // I'll use the normal Barycenter...
00264       //
00265       const Vector barycenter          = f.Barycenter();
00266       const Vector barycenter_velocity = f.BarycenterVelocity();
00267       // 
00268       // mod_f already equalt to f
00269       // mod_f = f;
00270       //
00271       for (unsigned int k=0; k<f.size(); ++k) {
00272         mod_f[k].AddToPosition(-barycenter);
00273         mod_f[k].AddToVelocity(-barycenter_velocity);
00274       }
00275       //
00276       mapped_table.load(mod_f,skip_JPL_planets);
00277     } else {
00278       mapped_table.load(f,skip_JPL_planets);
00279     }
00280     
00281     // note: mapped(i,j) = something(i) - something(j)
00282     
00283     /* 
00284        {
00285        // test        
00286        for (unsigned int i=0;i<size-1;++i) {
00287        for (unsigned int j=i+1;j<size;++j) {
00288        if (zero_mass[i] && zero_mass[j]) continue;
00289        if (skip[i] && skip[j]) continue;
00290        if (fabs((mapped_table.Distance(i,j)-f[i].DistanceVector(f[j]).Length())/(mapped_table.Distance(i,j))) > 1.0e-10) {
00291        fprintf(stderr,
00292        "PROBLEM-1!! d(%10s,%10s): mapped=%20.12f num=%20.12f\n",
00293        f[i].name().c_str(),
00294        f[j].name().c_str(),
00295        mapped_table.Distance(i,j),
00296        f[i].DistanceVector(f[j]).Length());
00297        }
00298        if (fabs((mapped_table.DistanceVector(i,j).x+f[i].DistanceVector(f[j]).x)/mapped_table.DistanceVector(i,j).x) > 1.0e-10) {
00299        fprintf(stderr,
00300        "PROBLEM-2!! d(%10s,%10s): mapped=%20.12f num=%20.12f\n",
00301        f[i].name().c_str(),
00302        f[j].name().c_str(),
00303        mapped_table.DistanceVector(i,j).x,
00304        f[i].DistanceVector(f[j]).x);
00305        }
00306        }
00307        }
00308        }
00309     */
00310     
00311     if (include_multipole_moments) {
00312       axis.resize(size);
00313       x_axis.resize(size);
00314       R1.resize(size);
00315       R2.resize(size);
00316       R3.resize(size);
00317       R4.resize(size);
00318       //
00319       {
00320         Angle alpha,delta,meridian;
00321         Vector planet_axis;
00322         Vector planet_x_axis;
00323         for (unsigned int i=0;i<size;++i) {
00324           if (f[i].JPLPlanet() == NONE) continue;
00325           alpha_delta_meridian(f[i].JPLPlanet(),f,alpha,delta,meridian);
00326           //
00327           planet_axis.Set(0,0,1);
00328           planet_x_axis.Set(1,0,0);
00329           // glRotated(90.0+(180.0/pi)*W.GetRad(),0,0,1);
00330           // no halfpi needed!
00331           planet_axis.rotate(meridian.GetRad(),0,0);
00332           planet_x_axis.rotate(meridian.GetRad(),0,0);
00333           // glRotated(90.0-(180.0/pi)*delta.GetRad(),1,0,0);
00334           planet_axis.rotate(0,halfpi-delta.GetRad(),0);
00335           planet_x_axis.rotate(0,halfpi-delta.GetRad(),0);
00336           // glRotated(90.0+(180.0/pi)*alpha.GetRad(),0,0,1);
00337           planet_axis.rotate(halfpi+alpha.GetRad(),0,0);
00338           planet_x_axis.rotate(halfpi+alpha.GetRad(),0,0);
00339           //
00340           if (universe->GetReferenceSystem() == ECLIPTIC) {
00341             // glRotated(-(180.0/pi)*obleq_J2000().GetRad(),1,0,0);
00342             planet_axis.rotate(0,-obleq_J2000().GetRad(),0);
00343             planet_x_axis.rotate(0,-obleq_J2000().GetRad(),0);
00344           }
00345           //
00346           axis[i] = planet_axis;
00347           x_axis[i] = planet_x_axis;
00348         }
00349       }
00350       //
00351       {
00352         double R;
00353         for (unsigned int i=0;i<size;++i) {
00354           R = f[i].radius();
00355           if (R1[i] != R) {
00356             R1[i] = R;
00357             R2[i] = R*R;
00358             R3[i] = R2[i]*R;
00359             R4[i] = R3[i]*R;
00360           }
00361         }
00362       }
00363     }
00364     
00365     for (unsigned int i=0;i<size-1;++i) {
00366       for (unsigned int j=i+1;j<size;++j) {
00367         
00368         if (zero_mass[i] && zero_mass[j]) continue;
00369         
00370         if (skip[i] && skip[j]) continue;
00371         
00372         if (mapped_table.Distance(i,j) < (std::numeric_limits<double>::min() * 1.0e3)) {
00373           ORSA_WARNING("two objects in the same position! (%s and %s)",f[i].name().c_str(),f[j].name().c_str());
00374           continue;
00375         }
00376         
00377         a_newton[i] += mapped_table.DistanceVectorOverDistanceCube(j,i) * f[j].mu();
00378         a_newton[j] += mapped_table.DistanceVectorOverDistanceCube(i,j) * f[i].mu();
00379       }
00380     }
00381     
00382     if (include_relativistic_effects) {
00383       
00384       for (unsigned int i=0;i<size-1;++i) {     
00385         for (unsigned int j=i+1;j<size;++j) {
00386           
00387           /* 
00388              for (unsigned int i=0;i<size;++i) {        
00389              for (unsigned int j=0;j<size;++j) {
00390           */
00391           
00392           // if (i == j) continue;
00393           
00394           if (zero_mass[i] && zero_mass[j]) continue;
00395           
00396           if (skip[i] && skip[j]) continue;
00397           
00398           // Now mod_f should be the Frame in the relativistic barycenter reference frame
00399           
00400           // Relativistic interaction, from:
00401           // "Explanatory Supplement to the Astronomical Almanac",
00402           // Edited by P. Kenneth Seidelmann, U.S. Naval Observatory,
00403           // 2nd edition, 1992, sect. 5.211, pag. 281.
00404           
00405           double sum_ik = 0.0;
00406           double sum_jk = 0.0;
00407           for (unsigned int k=0; k<size; ++k) {
00408             if (i != k) sum_ik += f[k].mu()*mapped_table.OneOverDistance(i,k);
00409             if (j != k) sum_jk += f[k].mu()*mapped_table.OneOverDistance(j,k);
00410           }
00411           
00412           a_relativity[i] += one_over_c2*f[j].mu()*mapped_table.DistanceVectorOverDistanceCube(j,i)*( // missing Newton term...
00413                                                                                                      - 4.0*sum_ik
00414                                                                                                      - 1.0*sum_jk
00415                                                                                                      + mod_f[i].velocity().LengthSquared()
00416                                                                                                      + 2.0*mod_f[j].velocity().LengthSquared()
00417                                                                                                      - 4.0*mod_f[i].velocity()*mod_f[j].velocity()
00418                                                                                                      - 1.5*pow(mapped_table.DistanceVector(i,j)*mod_f[j].velocity()*mapped_table.OneOverDistance(i,j),2)
00419                                                                                                      + 0.5*mapped_table.DistanceVector(j,i)*a_newton[j]
00420                                                                                                      )
00421             + one_over_c2*f[j].mu()*mapped_table.OneOverDistanceCube(i,j)*((mapped_table.DistanceVector(i,j)*(4.0*mod_f[i].velocity()-3.0*mod_f[j].velocity()))*(mod_f[i].velocity()-mod_f[j].velocity()))
00422             + 3.5*one_over_c2*f[j].mu()*mapped_table.OneOverDistance(i,j)*a_newton[j];
00423           
00424           a_relativity[j] += one_over_c2*f[i].mu()*mapped_table.DistanceVectorOverDistanceCube(i,j)*( // missing Newton term...
00425                                                                                                      - 4.0*sum_jk
00426                                                                                                      - 1.0*sum_ik
00427                                                                                                      + mod_f[j].velocity().LengthSquared()
00428                                                                                                      + 2.0*mod_f[i].velocity().LengthSquared()
00429                                                                                                      - 4.0*mod_f[j].velocity()*mod_f[i].velocity()
00430                                                                                                      - 1.5*pow(mapped_table.DistanceVector(j,i)*mod_f[i].velocity()*mapped_table.OneOverDistance(j,i),2)
00431                                                                                                      + 0.5*mapped_table.DistanceVector(i,j)*a_newton[i]
00432                                                                                                      )
00433             + one_over_c2*f[i].mu()*mapped_table.OneOverDistanceCube(j,i)*((mapped_table.DistanceVector(j,i)*(4.0*mod_f[j].velocity()-3.0*mod_f[i].velocity()))*(mod_f[j].velocity()-mod_f[i].velocity()))
00434             + 3.5*one_over_c2*f[i].mu()*mapped_table.OneOverDistance(j,i)*a_newton[i];
00435         }
00436       }
00437     }
00438     
00439     if (include_fast_relativistic_effects) {
00440       for (unsigned int i=0;i<size-1;++i) {     
00441         for (unsigned int j=i+1;j<size;++j) {
00442           
00443           // if (i == j) continue;
00444           
00445           if (zero_mass[i] && zero_mass[j]) continue;
00446           
00447           if (skip[i] && skip[j]) continue;
00448           
00449           // tests...
00450           // good old code...
00451           
00452           // const Vector r = mod_f[i].position() - mod_f[j].position();
00453           const Vector r = mapped_table.DistanceVector(i,j);
00454           //
00455           const Vector v  = mod_f[i].velocity() - mod_f[j].velocity();
00456           const double v2 = v.LengthSquared();
00457           
00458           // see http://www.projectpluto.com/relativi.htm
00459           //
00460           // acc = acc_{newton} - \frac{G M}{r^3 c^2} ( (\frac{4 G M}{r} - v^2) \vec(r) + 4 (v \cdot r) \vec(v) )
00461           //
00462           
00463           a_relativity[i] += one_over_c2*f[j].mu()*mapped_table.OneOverDistanceCube(i,j)*((4.0*f[j].mu()*mapped_table.OneOverDistance(i,j)-v2)*r + 
00464                                                                                           4.0*(r*v)*v);
00465           a_relativity[j] -= one_over_c2*f[i].mu()*mapped_table.OneOverDistanceCube(i,j)*((4.0*f[i].mu()*mapped_table.OneOverDistance(i,j)-v2)*r +
00466                                                                                           4.0*(r*v)*v);
00467         }
00468       }
00469     }
00470     
00471     if (include_multipole_moments) {
00472       
00473       /* for (unsigned int i=0;i<size-1;++i) {
00474          for (unsigned int j=i+1;j<size;++j) {
00475       */
00476       // scan all i's and j's, becaues i is always the extended body,
00477       // while j is always the point mass body
00478       for (unsigned int i=0;i<size;++i) {
00479         for (unsigned int j=0;j<size;++j) {
00480           
00481           if (i == j) continue;
00482           
00483           if (zero_mass[i] && zero_mass[j]) continue;
00484           
00485           if (skip[i] && skip[j]) continue;
00486           
00487           /* 
00488              const Vector      d = mapped_table.DistanceVector(i,j);
00489              const Vector unit_d = d.Normalized();
00490              
00491              const double l2 = mapped_table.Distance2(i,j);
00492              const double l3 = mapped_table.Distance3(i,j);
00493              const double l4 = mapped_table.Distance4(i,j);
00494           */
00495           
00496           if ((R1[i] > 0.0) && (fabs(f[i].J2()) > 0.0) && (f[i].JPLPlanet() != NONE)) {
00497             
00498             // cerr << "figure: extended body: " << f[i].name() << "  point-mass body: " << f[j].name() << endl;
00499             
00500             const Vector unit_x_local = mapped_table.DistanceVector(j,i).Normalized();
00501             const Vector unit_y_local = ExternalProduct(axis[i],unit_x_local);
00502             const Vector unit_z_local = ExternalProduct(unit_x_local,unit_y_local);
00503             
00504             const Vector tmp_y_axis = ExternalProduct(axis[i],x_axis[i]);
00505             const double tmp_uy_y   = unit_y_local*tmp_y_axis;
00506             const double tmp_uy_mx  = unit_y_local*x_axis[i];
00507             const double lambda = atan2(tmp_uy_mx,tmp_uy_y);
00508             //
00509             double s1l,c1l; orsa::sincos(lambda,s1l,c1l);
00510             double s2l,c2l; orsa::sincos(lambda,s2l,c2l);
00511             double s3l,c3l; orsa::sincos(lambda,s3l,c3l);
00512             double s4l,c4l; orsa::sincos(lambda,s4l,c4l);
00513             
00514             const double theta = acos(unit_x_local*axis[i]);
00515             //
00516             double st,ct; orsa::sincos(theta,st,ct);
00517             const double cosec_theta = 1.0/st;      
00518             
00519             const double ar  = R1[i]*mapped_table.OneOverDistance(i,j);
00520             const double ar2 = ar*ar;
00521             const double ar3 = ar*ar2;
00522             const double ar4 = ar*ar3;
00523             
00524             const Legendre l(axis[i]*unit_x_local);
00525             
00526             // cerr << "x: " << axis[i]*unit_x_local << endl;
00527             
00528             // const double base_coefficient = - f[j].mu() * mapped_table.OneOverDistanceSquare(i,j);
00529             const double base_coefficient = - mapped_table.OneOverDistanceSquare(i,j);
00530             
00531             const double acc_x_local = 
00532               base_coefficient *
00533               ( // zonal harmonics
00534                f[i].J2()*ar2*3.0*l.P2 +
00535                f[i].J3()*ar3*4.0*l.P3 + 
00536                f[i].J4()*ar4*5.0*l.P4 +
00537                // tesseral harmonics
00538                ar2*(-3.0*l.P22*(f[i].C22()*c2l) ) +
00539                ar3*(-4.0*l.P31*(f[i].C31()*c1l+f[i].S31()*s1l) +
00540                     -4.0*l.P32*(f[i].C32()*c2l+f[i].S32()*s2l) +
00541                     -4.0*l.P33*(f[i].C33()*c3l+f[i].S33()*s3l) ) +
00542                ar4*(-5.0*l.P41*(f[i].C41()*c1l+f[i].S41()*s1l) +
00543                     -5.0*l.P42*(f[i].C42()*c2l+f[i].S42()*s2l) +
00544                     -5.0*l.P43*(f[i].C43()*c3l+f[i].S43()*s3l) + 
00545                     -5.0*l.P44*(f[i].C44()*c4l+f[i].S44()*s4l) ) );
00546             
00547             const double acc_y_local = 
00548               base_coefficient *
00549               ( // zonal harmonics
00550                // no zonal harmonics along the y direction
00551                // tesseral harmonics
00552                ar2*(2.0*cosec_theta*l.P22*(-f[i].C22()*s2l) ) +
00553                ar3*(1.0*cosec_theta*l.P31*(-f[i].C31()*s1l+f[i].S31()*c1l) +
00554                     2.0*cosec_theta*l.P32*(-f[i].C32()*s2l+f[i].S32()*c2l) +
00555                     3.0*cosec_theta*l.P33*(-f[i].C33()*s3l+f[i].S33()*c3l) ) +
00556                ar4*(1.0*cosec_theta*l.P41*(-f[i].C41()*s1l+f[i].S41()*c1l) +
00557                     2.0*cosec_theta*l.P42*(-f[i].C42()*s2l+f[i].S42()*c2l) +
00558                     3.0*cosec_theta*l.P43*(-f[i].C43()*s3l+f[i].S43()*c3l) + 
00559                     4.0*cosec_theta*l.P44*(-f[i].C44()*s4l+f[i].S44()*c4l) ) );
00560             
00561             const double acc_z_local = 
00562               base_coefficient *
00563               ( // zonal harmonics
00564                f[i].J2()*ar2*(-st)*l.dP2 +
00565                f[i].J3()*ar3*(-st)*l.dP3 + 
00566                f[i].J4()*ar4*(-st)*l.dP4 +
00567                // tesseral harmonics
00568                ar2*(st*l.dP22*(f[i].C22()*c2l) ) +
00569                ar3*(st*l.dP31*(f[i].C31()*c1l+f[i].S31()*s1l) +
00570                     st*l.dP32*(f[i].C32()*c2l+f[i].S32()*s2l) +
00571                     st*l.dP33*(f[i].C33()*c3l+f[i].S33()*s3l) ) +
00572                ar4*(st*l.dP41*(f[i].C41()*c1l+f[i].S41()*s1l) +
00573                     st*l.dP42*(f[i].C42()*c2l+f[i].S42()*s2l) +
00574                     st*l.dP43*(f[i].C43()*c3l+f[i].S43()*s3l) + 
00575                     st*l.dP44*(f[i].C44()*c4l+f[i].S44()*s4l) ) );
00576             
00577             const Vector local_acc =
00578               unit_x_local * acc_x_local +
00579               unit_y_local * acc_y_local +
00580               unit_z_local * acc_z_local;
00581             
00582             // cerr << "F: i: " << f[i].name() << " j: " << f[j].name() << " ax: " << acc_x_local << " tot: " << local_acc.Length() << endl;
00583             
00584             a_multipoles[i] += f[j].mu() * local_acc;
00585             
00586             // reaction
00587             a_multipoles[j] -= f[i].mu() * local_acc;
00588             
00589           }
00590         }
00591       }
00592     }
00593     
00594     for (unsigned int i=0;i<size;++i) {
00595       a[i] = a_newton[i] + (a_multipoles[i] + a_relativity[i]);
00596     }
00597     
00598     // done
00599   }
00600   
00601   double Newton::PotentialEnergy(const Frame &f) {
00602     
00603     if (f.size() < 2) return(0.0);
00604     
00605     double energy = 0.0;
00606     
00607     unsigned int i,j;
00608     
00609     Vector d;
00610     
00611     double l;
00612     
00613     for (i=1;i<f.size();++i) {
00614       
00615       if (f[i].mu()==0) continue;
00616       
00617       for (j=0;j<i;++j) {
00618         
00619         if (f[j].mu()==0) continue;
00620         
00621         d = f[i].DistanceVector(f[j]);
00622         
00623         l = d.Length();
00624         
00625         if (d.IsZero()) {
00626           ORSA_WARNING("two objects in the same position! (%s and %s)",f[i].name().c_str(),f[j].name().c_str());
00627           continue;
00628         }
00629         
00630         energy -= f[i].mu()*f[j].mass()/l;
00631         
00632       } 
00633     }  
00634     
00635     return (energy);
00636   }
00637   
00638   // Relativistic
00639   
00640   Relativistic::Relativistic() : Interaction(), g(GetG()), c_2(GetC()*GetC()) {
00641     
00642   }
00643   
00644   Relativistic::Relativistic(const Relativistic &) : Interaction(), g(GetG()), c_2(GetC()*GetC()) {
00645     
00646   }
00647   
00648   Interaction * Relativistic::clone() const {
00649     return new Relativistic(*this);
00650   }
00651   
00652   // see http://www.projectpluto.com/relativi.htm
00653   //
00654   // acc = acc_{newton} - \frac{G M}{r^3 c^2} ( (\frac{4 G M}{r} - v^2) \vec(r) + 4 (v \cdot r) \vec(v) )
00655   //
00656   void Relativistic::Acceleration(const Frame &f, vector<Vector> &a) {
00657     
00658     if (f.size() < 2) return;
00659     
00660     a.resize(f.size());
00661     
00662     unsigned int i,j;
00663     
00664     for (i=0;i<a.size();++i)
00665       a[i].Set(0,0,0);
00666     
00667     Vector r;
00668     Vector v;
00669     
00670     double r_1, r_3;
00671     
00672     for (i=1;i<f.size();++i) {
00673       
00674       for (j=0;j<i;++j) {
00675         
00676         /*** acc. on 'i' due to 'j' ***/
00677         
00678         if ((f[i].mass()==0) && (f[j].mass()==0)) continue;
00679         
00680         // 'i' - 'j'
00681         r = f[i].position() - f[j].position();
00682         v = f[i].velocity() - f[j].velocity();
00683         
00684         if (r.IsZero()) {
00685           ORSA_WARNING("two objects in the same position! (%s and %s)",f[i].name().c_str(),f[j].name().c_str());
00686           continue;
00687         }
00688         
00689         r_1 = r.Length();
00690         r_3 = r.LengthSquared()*r_1;
00691         
00692         a[i] -= f[j].mass()/r_3*r + f[j].mass()/(r_3*c_2)*((4.0*g*f[j].mass()/r_1-v.LengthSquared())*r + 4.0*(r*v)*v);
00693         
00694         // now, acc. on 'j' due to 'i' so r -> -r and v -> -v ***/
00695         //
00696         a[j] += f[i].mass()/r_3*r + f[i].mass()/(r_3*c_2)*((4.0*g*f[i].mass()/r_1-v.LengthSquared())*r + 4.0*(r*v)*v);
00697         
00698       } 
00699       
00700     } 
00701     
00702     for (i=0;i<a.size();++i) a[i] *= g;
00703   }
00704   
00705   double Relativistic::PotentialEnergy(const Frame &f) {
00706     Newton newton;
00707     return (newton.PotentialEnergy(f));
00708   }
00709   
00710   // ArmonicOscillator
00711   
00712   ArmonicOscillator::ArmonicOscillator(const double free_length_in, const double k_in) : Interaction(), free_length(free_length_in), k(k_in) {
00713     
00714   }
00715   
00716   ArmonicOscillator::ArmonicOscillator(const ArmonicOscillator & i) : Interaction(), free_length(i.free_length), k(i.k) {
00717     
00718   }
00719   
00720   Interaction * ArmonicOscillator::clone() const {
00721     return new ArmonicOscillator(*this);
00722   }
00723   
00724   void ArmonicOscillator::Acceleration(const Frame& f, vector<Vector>& a) {
00725     
00726     if (f.size() < 2) return;
00727     
00728     a.resize(f.size());
00729     
00730     unsigned int i,j;
00731     
00732     Vector d, da;
00733 
00734     double ls;
00735     
00736     for (i=0;i<f.size();++i)
00737       a[i].Set(0,0,0);
00738     
00739     for (i=1;i<f.size();++i) {
00740       
00741       if (f[i].mass()==0) continue;
00742       
00743       for (j=0;j<i;++j) {
00744         
00745         // if ((f[i].mass==0) || (f[j].mass==0)) continue;
00746         
00747         // assert(i!=j);
00748         
00749         // d  =  f[i].position;
00750         // d -=  f[j].position;
00751         
00752         d = f[i].DistanceVector(f[j]);
00753 
00754         ls = d.Length();
00755         
00756         if (d.IsZero()) {
00757           ORSA_WARNING("two objects in the same position! (%s and %s)",f[i].name().c_str(),f[j].name().c_str());
00758           continue;
00759         }
00760         
00761         da = d*(ls-free_length)/ls;
00762         
00763         a[i] += da;
00764         a[j] -= da;
00765         
00766       } 
00767     }  
00768     
00769     for (i=0;i<a.size();++i) 
00770       if (f[i].mass() != 0) a[i] *= k/f[i].mass();
00771     
00772   }
00773   
00774   double ArmonicOscillator::PotentialEnergy(const Frame &f) {
00775     
00776     double energy=0.0;
00777     
00778     if (f.size() < 2) return(0.0);
00779     
00780     unsigned int i,j;
00781     
00782     Vector d;
00783     
00784     double ls;
00785     
00786     for (i=1;i<f.size();++i) {
00787       
00788       if (f[i].mass()==0) continue;
00789       
00790       for (j=0;j<i;++j) {
00791         
00792         // if ((f[i].mass==0) && (f[j].mass==0)) continue;
00793         
00794         // d  =  f[i].position;
00795         // d -=  f[j].position;
00796         
00797         d = f[i].DistanceVector(f[j]);
00798         
00799         if (d.IsZero()) {
00800           ORSA_WARNING("two objects in the same position! (%s and %s)",f[i].name().c_str(),f[j].name().c_str());
00801           continue;
00802         }
00803         
00804         // ls = secure_pow(d.Length()-free_length,2);
00805         ls = pow(d.Length()-free_length,2); // don't use the secure version for performance reasons
00806         
00807         energy += ls/2.0;
00808         
00809       } 
00810     }  
00811     
00812     return (energy*k);
00813     
00814   }
00815   
00816   // GalacticPotentialAllen
00817   // see: http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1991RMxAA..22..255A&db_key=AST&high=3e8aa73e3d07980
00818   /* @ARTICLE{1991RMxAA..22..255A,
00819      author = {{Allen}, C. and {Santillan}, A.},
00820      title = "{An improved model of the galactic mass distribution for orbit computations}",
00821      journal = {Revista Mexicana de Astronomia y Astrofisica},
00822      year = 1991,
00823      month = oct,
00824      volume = 22,
00825      pages = {255-263},
00826      adsurl = {http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1991RMxAA..22..255A&db_key=AST},
00827      adsnote = {Provided by the NASA Astrophysics Data System}
00828      }
00829   */
00830   
00831   GalacticPotentialAllen::GalacticPotentialAllen() : Interaction() {
00832     
00833     g = GetG();
00834     
00835     mb = FromUnits(1.40592e10,   MSUN);
00836     bb = FromUnits(0.3873e3,     PARSEC);
00837     md = FromUnits(8.5608e10,    MSUN);     
00838     ad = FromUnits(5.3178e3,     PARSEC);     
00839     bd = FromUnits(0.25e3,       PARSEC);
00840     mh = FromUnits(10.3836745e10,MSUN);
00841     ah = FromUnits(12.0e3,       PARSEC); 
00842     
00843   }
00844   
00845   
00846   GalacticPotentialAllen::GalacticPotentialAllen(const GalacticPotentialAllen &) : Interaction() {
00847     
00848     g = GetG();
00849     
00850     mb = FromUnits(1.40592e10,   MSUN);
00851     bb = FromUnits(0.3873e3,     PARSEC);
00852     md = FromUnits(8.5608e10,    MSUN);     
00853     ad = FromUnits(5.3178e3,     PARSEC);     
00854     bd = FromUnits(0.25e3,       PARSEC);
00855     mh = FromUnits(10.3836745e10,MSUN);
00856     ah = FromUnits(12.0e3,       PARSEC); 
00857     
00858   }
00859   
00860   Interaction * GalacticPotentialAllen::clone() const {
00861     return new GalacticPotentialAllen(*this);
00862   }
00863   
00864   void GalacticPotentialAllen::Acceleration(const Frame &f, vector<Vector> &a) {
00865     
00866     a.resize(f.size());
00867     
00868     unsigned int i;
00869     for (i=0;i<a.size();++i) a[i].Set(0,0,0);
00870     
00871     double r2,z2,r,rho;
00872     double fbr,fdr,fhr,fhrc;
00873     double fbz,fdz,fhz,fhzc;
00874     double fr,fz;
00875     Vector x;
00876     
00877     for (i=0;i<f.size();++i) {
00878       
00879       // if (!f[i].StillAlive()) continue;
00880       
00881       x = f[i].position();
00882       
00883       // auxiliary components
00884       r2  = x.x*x.x+x.y*x.y;
00885       z2  = x.z*x.z;
00886       r   = sqrt(r2);
00887       rho = x.Length();
00888       
00889       // accel. components along the r direction
00890       fbr  = -(mb*r/(secure_pow(r2+z2+bb*bb,1.5)));
00891       fdr  = -(md*r/(secure_pow(r2+secure_pow(ad+sqrt(z2+bd*bd),2),1.5)));
00892       fhr  = 1.02*mh*r/(ah*ah*secure_pow(1+1.0/secure_pow(rho/ah,1.02),2)*secure_pow(rho/ah,2.02)*rho);
00893       fhrc = -(mh/(1.02*ah)*(1.0404*(r*secure_pow(rho/ah,0.02))/(ah*rho*secure_pow(1.0+secure_pow(rho/ah,1.02),2))+1.02*(r*secure_pow(rho/ah,0.02))/(ah*rho*(1.0+secure_pow(rho/ah,1.02)))));
00894       
00895       // accel. components along the z direction
00896       fbz  = -(mb*x.z/(secure_pow(r2+z2+bb*bb,1.5)));
00897       fdz  = -(md*x.z*(ad+sqrt(z2+bd*bd))/(sqrt(z2+bd*bd)*secure_pow(r2+secure_pow(ad+sqrt(z2+bd*bd),2),1.5)));
00898       fhz  = 1.02*mh*x.z/(ah*ah*secure_pow(1.0+1.0/secure_pow(rho/ah,1.02),2)*secure_pow(rho/ah,2.02)*rho);
00899       fhzc = -(mh/(1.02*ah)*(1.0404*(x.z*secure_pow(rho/ah,0.02))/(ah*rho*secure_pow(1.0+secure_pow(rho/ah,1.02),2))+1.02*(x.z*secure_pow(rho/ah,0.02))/(ah*rho*(1.0+secure_pow(rho/ah,1.02)))));
00900       
00901       fr=fbr+fdr+fhr+fhrc;
00902       fz=fbz+fdz+fhz+fhzc;
00903       
00904       a[i].x = fr*x.x/r;
00905       a[i].y = fr*x.y/r;
00906       a[i].z = fz;
00907     } 
00908     
00909     for (i=0;i<a.size();++i) a[i] *= g;
00910     
00911     // for (i=0;i<a.size();++i) printf("a[%i] = %g %g %g\n",i,a[i].x,a[i].y,a[i].z);
00912     
00913   }
00914   
00915   double GalacticPotentialAllen::PotentialEnergy(const Frame &f) {
00916     
00917     // to be tested...
00918     
00919     double energy = 0.0;
00920     
00921     double pb,pd,ph;
00922     double r2,z2,r,rho;
00923     unsigned int i;
00924     Vector x;
00925     
00926     for (i=0;i<f.size();++i) {
00927       
00928       x = f[i].position();
00929       
00930       // auxiliary components
00931       r2  = x.x*x.x+x.y*x.y;
00932       z2  = x.z*x.z;
00933       r   = sqrt(r2);
00934       rho = x.Length();
00935       
00936       pb = mb/sqrt(r2+z2+bb*bb);
00937       pd = md/sqrt(r2+secure_pow(ad+sqrt(z2+bd*bd),2));
00938       ph = (mh/(rho/ah))*(secure_pow(rho/ah,2.02))/(1.0+secure_pow(rho/ah,1.02));
00939       
00940       energy -= (pb+pd+ph);
00941       
00942     } 
00943     
00944     return (energy*g);
00945   }
00946   
00947   // JPLPlanetsNewton
00948   
00949   JPLPlanetsNewton::JPLPlanetsNewton(list<JPL_planets> & l_in) : Interaction(), l(l_in) {
00950     if (universe->GetUniverseType() != Real) {
00951       cerr << "error: using the JPLPlanetsNewton interaction in a non-Real universe!" << endl;
00952       exit(0);
00953     }
00954     g = GetG();
00955   }
00956   
00957   JPLPlanetsNewton::JPLPlanetsNewton(const JPLPlanetsNewton & i) : Interaction(), l(i.l) {
00958     if (universe->GetUniverseType() != Real) {
00959       cerr << "error: using the JPLPlanetsNewton interaction in a non-Real universe!" << endl;
00960       exit(0);
00961     }
00962     g = GetG();
00963   }
00964   
00965   Interaction * JPLPlanetsNewton::clone() const {
00966     return new JPLPlanetsNewton(*this);
00967   }
00968   
00969   void JPLPlanetsNewton::Acceleration(const Frame & f, vector<Vector> & a) {
00970     
00971     a.resize(f.size());
00972     
00973     if (planets_frame.GetDate() != f.GetDate()) SetupSolarSystem(planets_frame,l,f);
00974     
00975     for (unsigned int i=0;i<a.size();++i)
00976       a[i].Set(0,0,0);
00977     
00978     for (unsigned int i=0;i<f.size();++i) {
00979       if (f[i].mass() > 0) {
00980         ORSA_ERROR("using the JPLPlanetsNewton interaction with massive objects!");
00981         // exit(0);
00982         return;
00983       }
00984     }
00985     
00986     Vector d;
00987     
00988     double l;
00989     
00990     for (unsigned int i=0;i<f.size();++i) {
00991       
00992       // if (!f[i].StillAlive()) continue;
00993       
00994       for (unsigned int j=0;j<planets_frame.size();++j) {
00995         
00996         // d  =              f[i].position;
00997         // d -=  planets_frame[j].position;
00998         
00999         d = f[i].DistanceVector(planets_frame[j]);
01000         
01001         // ls = d.LengthSquared();
01002         l = d.Length();
01003         
01004         if (d.IsZero()) {
01005           ORSA_WARNING("two objects in the same position! (%s and %s)",
01006                   f[i].name().c_str(),
01007                   planets_frame[j].name().c_str());
01008           continue;
01009         }
01010         
01011         // d *= secure_pow(ls,-1.5);
01012         d /= l*l*l;
01013         
01014         a[i] += d * planets_frame[j].mass();
01015         
01016       }
01017     }
01018     
01019     for (unsigned int i=0;i<a.size();++i) a[i] *= g;
01020   }
01021   
01022   double JPLPlanetsNewton::PotentialEnergy(const Frame & f) {
01023     SetupSolarSystem(planets_frame,l,f);
01024     return newton_itg.PotentialEnergy(planets_frame);
01025   }
01026   
01027 } // namespace orsa

Generated on Thu Jul 13 06:45:22 2006 for liborsa by  doxygen 1.4.7