1 /*****************************************************************************\
2  * Pluto.cpp
3  *
4  * Pluto is a class that can calculate the orbit of Pluto.
5  *
6  * author: mark huss (mark@mhuss.com)
7  * Based on Bill Gray's open-source code at projectpluto.com
8  *
9 \*****************************************************************************/
10 
11 #include <math.h>
12 #include "AstroOps.h"
13 #include "Pluto.h"
14 #include "PlutoTerms.h"  // data extracted from vsop.bin file
15 
16 // t is in julian centuries from J2000.0
17 
calcAllLocs(double & lat,double & lon,double & rad,const double t)18 void Pluto::calcAllLocs( double& lat, double& lon, double& rad, const double t) {
19 
20   // jupiter's mean longitude
21   double mlJup =  Astro::toRadians(34.35 + 3034.9057 * t);
22 
23   // saturn's mean longitude
24   double mlSat =  Astro::toRadians(50.08 + 1222.1138 * t);
25 
26   // pluto's mean longitude
27   double mlPl = Astro::toRadians(238.96 +  144.9600 * t);
28 
29   // use local vars for retvals, hoping to encourage FP reg optimizations
30   double lon_ = 238.956785 + 144.96 * t;
31   double lat_ = -3.908202;
32   double rad_ = 407.247248;  // temporarily in tenths of AUs; fixed at the end
33 
34   double arg;
35   for( int i=0; i<7; i++ ) {
36     if( i == 6)
37       arg = mlJup - mlPl;
38     else
39       arg = (double)(i + 1) * mlPl;
40 
41     double cosArg = cos(arg) * 1.e-6;
42     double sinArg = sin(arg) * 1.e-6;
43     long* plc = plutoLongCoeff[i];
44 
45     lon_ += (double)(plc[0]) * sinArg + (double)(plc[1]) * cosArg;
46     lat_ += (double)(plc[2]) * sinArg + (double)(plc[3]) * cosArg;
47     rad_   += (double)(plc[4]) * sinArg + (double)(plc[5]) * cosArg;
48   }
49 
50   PlutoCoeffs* pc = plutoCoeff;
51   for( int i=0; i<N_COEFFS; i++ ) {
52     if( pc->lon_a || pc->lon_b ||
53         pc->lat_a || pc->lat_b ||
54         pc->rad_a || pc->rad_b)
55     {
56       if( 0 == pc->j )
57         arg = 0.;
58       else
59         arg = ((1 == pc->j) ? mlJup : mlJup * (double)pc->j);
60 
61       if( pc->s < 0)
62         arg -= (-1 == pc->s) ? mlSat : mlSat + mlSat;
63 
64       if( pc->s > 0)
65         arg += (1 == pc->s) ? mlSat : mlSat + mlSat;
66 
67       if( pc->p)
68         arg += mlPl * (double)pc->p;
69 
70       double cosArg = cos(arg) * 1.e-6;
71       double sinArg = sin(arg) * 1.e-6;
72       lon_ += sinArg * (double)(pc->lon_a) + cosArg * (double)(pc->lon_b);
73       lat_ += sinArg * (double)(pc->lat_a) + cosArg * (double)(pc->lat_b);
74       rad_   += sinArg * (double)(pc->rad_a) + cosArg * (double)(pc->rad_b);
75     }
76     pc++;
77   }
78   lon = Astro::toRadians(lon_);
79   lat = Astro::toRadians(lat_);
80   rad = rad_ / 10.;                  // convert back to AUs
81 }
82 
83