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)18void 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