AstroToolbox/Pl_Eph_An.cpp

   1 // ------------------------------------------------------------------------ //
   2 // This source file is part of the 'ESA Advanced Concepts Team's                        //
   3 // Space Mechanics Toolbox' software.                                       //
   4 //                                                                          //
   5 // The source files are for research use only,                              //
   6 // and are distributed WITHOUT ANY WARRANTY. Use them on your own risk.     //
   7 //                                                                          //
   8 // Copyright (c) 2004-2007 European Space Agency                            //
   9 // ------------------------------------------------------------------------ //
  10
  11
  12 #include <cmath>
  13
  14 #include "Pl_Eph_An.h"
  15 #include "Astro_Functions.h"
  16
  17 // TODO: place constants in separate class and share them.
  18
  19 using namespace std;
  20
  21 void Planet_Ephemerides_Analytical (const double &mjd2000,
  22                                                                         const int &planet,
  23                                                                         double *position,
  24                                                                         double *velocity)
  25 {
  26         const double RAD = M_PI/180.0;
  27         const double KM_AU = 149597870.66; // Kilometers in an astronomical Unit
  28         //const double MuSun = 1.32712440018e+11; //Gravitational constant of Sun);
  29         const double MuSun = 1.327124280000000e+011;  //Gravitational constant of Sun);
  30         double Kepl_Par[6];
  31         double XM;
  32
  33         double T =(mjd2000 + 36525.00)/36525.00;
  34
  35         switch (planet)
  36         {
  37                 case(1):// Mercury
  38                         Kepl_Par[0]=(0.38709860);
  39                         Kepl_Par[1]=(0.205614210 + 0.000020460*T - 0.000000030*T*T);
  40                         Kepl_Par[2]=(7.002880555555555560 + 1.86083333333333333e-3*T - 1.83333333333333333e-5*T*T);
  41                         Kepl_Par[3]=(4.71459444444444444e+1 + 1.185208333333333330*T + 1.73888888888888889e-4*T*T);
  42                         Kepl_Par[4]=(2.87537527777777778e+1 + 3.70280555555555556e-1*T +1.20833333333333333e-4*T*T);
  43                         XM   = 1.49472515288888889e+5 + 6.38888888888888889e-6*T;
  44                         Kepl_Par[5]=(1.02279380555555556e2 + XM*T);
  45                 break;
  46                 case(2):// Venus
  47                         Kepl_Par[0]=(0.72333160);
  48                         Kepl_Par[1]=(0.006820690 - 0.000047740*T + 0.0000000910*T*T);
  49                         Kepl_Par[2]=(3.393630555555555560 + 1.00583333333333333e-3*T - 9.72222222222222222e-7*T*T);
  50                         Kepl_Par[3]=(7.57796472222222222e+1 + 8.9985e-1*T + 4.1e-4*T*T);
  51                         Kepl_Par[4]=(5.43841861111111111e+1 + 5.08186111111111111e-1*T -1.38638888888888889e-3*T*T);
  52                         XM =5.8517803875e+4 + 1.28605555555555556e-3*T;
  53                         Kepl_Par[5]=(2.12603219444444444e2 + XM*T);
  54                 break;
  55                 case(3):// Earth
  56                         Kepl_Par[0]=(1.000000230);
  57                         Kepl_Par[1]=(0.016751040 - 0.000041800*T - 0.0000001260*T*T);
  58                         Kepl_Par[2]=(0.00);
  59                         Kepl_Par[3]=(0.00);
  60                         Kepl_Par[4]=(1.01220833333333333e+2 + 1.7191750*T + 4.52777777777777778e-4*T*T + 3.33333333333333333e-6*T*T*T);
  61                          XM   = 3.599904975e+4 - 1.50277777777777778e-4*T - 3.33333333333333333e-6*T*T;
  62                         Kepl_Par[5]=(3.58475844444444444e2 + XM*T);
  63                 break;
  64                 case(4):// Mars
  65                         Kepl_Par[0]=(1.5236883990);
  66                         Kepl_Par[1]=(0.093312900 + 0.0000920640*T - 0.0000000770*T*T);
  67                         Kepl_Par[2]=(1.850333333333333330 - 6.75e-4*T + 1.26111111111111111e-5*T*T);
  68                         Kepl_Par[3]=(4.87864416666666667e+1 + 7.70991666666666667e-1*T - 1.38888888888888889e-6*T*T - 5.33333333333333333e-6*T*T*T);
  69                         Kepl_Par[4]=(2.85431761111111111e+2 + 1.069766666666666670*T +  1.3125e-4*T*T + 4.13888888888888889e-6*T*T*T);
  70                         XM   = 1.91398585e+4 + 1.80805555555555556e-4*T + 1.19444444444444444e-6*T*T;
  71                         Kepl_Par[5]=(3.19529425e2 + XM*T);
  72                 break;
  73                 case(5):// Jupiter
  74                         Kepl_Par[0]=(5.2025610);
  75                         Kepl_Par[1]=(0.048334750 + 0.000164180*T  - 0.00000046760*T*T -0.00000000170*T*T*T);
  76                         Kepl_Par[2]=(1.308736111111111110 - 5.69611111111111111e-3*T +  3.88888888888888889e-6*T*T);
  77                         Kepl_Par[3]=(9.94433861111111111e+1 + 1.010530*T + 3.52222222222222222e-4*T*T - 8.51111111111111111e-6*T*T*T);
  78                         Kepl_Par[4]=(2.73277541666666667e+2 + 5.99431666666666667e-1*T + 7.0405e-4*T*T + 5.07777777777777778e-6*T*T*T);
  79                         XM   = 3.03469202388888889e+3 - 7.21588888888888889e-4*T + 1.78444444444444444e-6*T*T;
  80                         Kepl_Par[5]=(2.25328327777777778e2 + XM*T);
  81                 break;
  82                 case(6):// Saturn
  83                         Kepl_Par[0]=(9.5547470);
  84                         Kepl_Par[1]=(0.055892320 - 0.00034550*T - 0.0000007280*T*T + 0.000000000740*T*T*T);
  85                         Kepl_Par[2]=(2.492519444444444440 - 3.91888888888888889e-3*T - 1.54888888888888889e-5*T*T + 4.44444444444444444e-8*T*T*T);
  86                         Kepl_Par[3]=(1.12790388888888889e+2 + 8.73195138888888889e-1*T -1.52180555555555556e-4*T*T - 5.30555555555555556e-6*T*T*T);
  87                         Kepl_Par[4]=(3.38307772222222222e+2 + 1.085220694444444440*T + 9.78541666666666667e-4*T*T + 9.91666666666666667e-6*T*T*T);
  88                         XM   = 1.22155146777777778e+3 - 5.01819444444444444e-4*T - 5.19444444444444444e-6*T*T;
  89                         Kepl_Par[5]=(1.75466216666666667e2 + XM*T);
  90                 break;
  91                 case(7):// Uranus
  92                         Kepl_Par[0]=(19.218140);
  93                         Kepl_Par[1]=(0.04634440 - 0.000026580*T + 0.0000000770*T*T);
  94                         Kepl_Par[2]=(7.72463888888888889e-1 + 6.25277777777777778e-4*T + 3.95e-5*T*T);
  95                         Kepl_Par[3]=(7.34770972222222222e+1 + 4.98667777777777778e-1*T + 1.31166666666666667e-3*T*T);
  96                         Kepl_Par[4]=(9.80715527777777778e+1 + 9.85765e-1*T - 1.07447222222222222e-3*T*T - 6.05555555555555556e-7*T*T*T);
  97                         XM   = 4.28379113055555556e+2 + 7.88444444444444444e-5*T + 1.11111111111111111e-9*T*T;
  98                         Kepl_Par[5]=(7.26488194444444444e1 + XM*T);
  99                 break;
 100                 case(8)://Neptune
 101                         Kepl_Par[0]=(30.109570);
 102                         Kepl_Par[1]=(0.008997040 + 0.0000063300*T - 0.0000000020*T*T);
 103                         Kepl_Par[2]=(1.779241666666666670 - 9.54361111111111111e-3*T - 9.11111111111111111e-6*T*T);
 104                         Kepl_Par[3]=(1.30681358333333333e+2 + 1.0989350*T + 2.49866666666666667e-4*T*T - 4.71777777777777778e-6*T*T*T);
 105                         Kepl_Par[4]=(2.76045966666666667e+2 + 3.25639444444444444e-1*T + 1.4095e-4*T*T + 4.11333333333333333e-6*T*T*T);
 106                         XM   = 2.18461339722222222e+2 - 7.03333333333333333e-5*T;
 107                         Kepl_Par[5]=(3.77306694444444444e1 + XM*T);
 108                 break;
 109                 case(9):// Pluto
 110                         //Fifth order polynomial least square fit generated by Dario Izzo
 111                         //(ESA ACT). JPL405 ephemerides (Charon-Pluto barycenter) have been used to produce the coefficients.
 112                         //This approximation should not be used outside the range 2000-2100;
 113                         T =mjd2000/36525.00;
 114                         Kepl_Par[0]=(39.34041961252520 + 4.33305138120726*T - 22.93749932403733*T*T + 48.76336720791873*T*T*T - 45.52494862462379*T*T*T*T + 15.55134951783384*T*T*T*T*T);
 115                         Kepl_Par[1]=(0.24617365396517 + 0.09198001742190*T - 0.57262288991447*T*T + 1.39163022881098*T*T*T - 1.46948451587683*T*T*T*T + 0.56164158721620*T*T*T*T*T);
 116                         Kepl_Par[2]=(17.16690003784702 - 0.49770248790479*T + 2.73751901890829*T*T - 6.26973695197547*T*T*T + 6.36276927397430*T*T*T*T - 2.37006911673031*T*T*T*T*T);
 117                         Kepl_Par[3]=(110.222019291707 + 1.551579150048*T - 9.701771291171*T*T + 25.730756810615*T*T*T - 30.140401383522*T*T*T*T + 12.796598193159 * T*T*T*T*T);
 118                         Kepl_Par[4]=(113.368933916592 + 9.436835192183*T - 35.762300003726*T*T + 48.966118351549*T*T*T - 19.384576636609*T*T*T*T - 3.362714022614 * T*T*T*T*T);
 119                         Kepl_Par[5]=(15.17008631634665 + 137.023166578486*T + 28.362805871736*T*T - 29.677368415909*T*T*T - 3.585159909117*T*T*T*T + 13.406844652829 * T*T*T*T*T);
 120                 break;
 121
 122         }
 123
 124         // conversion of AU into KM
 125         Kepl_Par[0] *= KM_AU;
 126
 127         // conversion of DEG into RAD
 128         Kepl_Par[2] *= RAD;
 129         Kepl_Par[3] *= RAD;
 130         Kepl_Par[4] *= RAD;
 131         Kepl_Par[5] *= RAD;
 132         Kepl_Par[5] =  fmod(Kepl_Par[5], 2.0*M_PI);
 133
 134         // Conversion from Mean Anomaly to Eccentric Anomaly via Kepler's equation
 135         Kepl_Par[5] = Mean2Eccentric(Kepl_Par[5], Kepl_Par[1]);
 136
 137         // Position and Velocity evaluation according to j2000 system
 138         Conversion(Kepl_Par, position, velocity, MuSun);
 139 }
 140
 141 void Custom_Eph(const double &jd,
 142                                 const double &epoch,
 143                                 const double keplerian[],
 144                                 double *position,
 145                                 double *velocity)
 146 {
 147         const double RAD = M_PI/180.0;
 148         const double KM_AU = 149597870.66; // Kilometers in an astronomical Unit
 149         const double muSUN = 1.32712428e+11;    // Gravitational constant of Sun
 150         double a,e,i,W,w,M,jdepoch,DT,n,E;
 151         double V[6];
 152
 153         a=keplerian[0]*KM_AU; // in km
 154     e=keplerian[1];
 155     i=keplerian[2];
 156         W=keplerian[3];
 157     w=keplerian[4];
 158     M=keplerian[5];
 159         jdepoch = epoch +2400000.5;
 160     DT=(jd-jdepoch)*86400;
 161     n=sqrt(muSUN/pow(a,3));
 162
 163     M=M/180.0*M_PI;
 164     M+=n*DT;
 165     M=fmod(M,2*M_PI);
 166     E=Mean2Eccentric(M,e);
 167         V[0] = a; V[1] = e; V[2] = i*RAD;
 168         V[3] = W*RAD; V[4] = w*RAD; V[5] = E;
 169
 170         Conversion(V,position,velocity,muSUN);
 171 }