/* ROCKS.CPP This code generates positions for 25 "rocks" (faint inner satellites of gas giants that can be modelled as precessing ellipses). The 25 "rocks" are: 505: Amalthea 706: 1986U7 Cordelia 808: 1989 N1 Proteus 514: 1979J2 Thebe 707: 1986U8 Ophelia 807: 1989 N2 Larissa 515: 1979J1 Adrastea 708: 1986U9 Bianca 805: 1989 N3 Despina 516: 1979J3 Metis 709: 1986U3 Cressida 806: 1989 N4 Galatea 710: 1986U6 Desdemona 804: 1989 N5 Thalassa 615: 1980S28 Atlas 711: 1986U2 Juliet 803: 1989 N6 Naiad 616: 1980S27 Prometheus 712: 1986U1 Portia 617: 1980S26 Pandora 713: 1986U4 Rosalind 618: 1981S13 Pan 714: 1986U5 Belinda 715: 1986U1 Puck 718: 1986U10 The orbital elements were provided by Mark Showalter, along with the following comments: (1) The listing refers to JPL ephemeris IDs JUP120, SAT080, URA039, and NEP022. If you select the SPICE ephemerides of the same name in Horizons, then you really should get almost identical results, because the SPICE files were generated using these exact elements. (2) Be aware that Saturn's moons Prometheus and Pandora are exhibiting peculiar variations in longitude that are not currently understood, probably related to interactions with the other nearby moons and rings. So don't expect your results for these bodies to be particularly accurate. Regards, Mark Showalter */ #include #define ROCK struct rock ROCK { int jpl_id; double epoch_jd, a, h, k, mean_lon0, p, q, apsis_rate, mean_motion; double node_rate, laplacian_pole_ra, laplacian_pole_dec; }; #define PI 3.14159265358979323846264338327950288 #define N_ROCKS 25 static const ROCK rocks[N_ROCKS] = { { 505, // Amalthea 2443937.5, // element epoch Julian date 1.813655740000000E+05, // a = semi-major axis (km) 1.605053686034257E-03, // h = e sin(periapsis longitude) 1.534237220548455E-03, // k = e cos(periapsis longitude) 1.927834989971504E+02 * PI / 180., // l = mean longitude (deg) -3.135065124778306E-03, // p = tan(i/2) sin(node) -6.070750467755681E-04, // q = tan(i/2) cos(node) 2.908806783226840E-05 * PI / 180., // apsis rate (deg/sec) 8.363793061627021E-03 * PI / 180., // mean motion (deg/sec) -2.898948153780610E-05 * PI / 180., // node rate (deg/sec) 2.680583409224650E+02 * PI / 180., // Laplacian plane pole ra (deg) 6.449442793849740E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 514, // Thebe 2443937.81706, // element epoch Julian date 2.218883060000000E+05, // a = semi-major axis (km) -1.407779640533111E-02, // h = e sin(periapsis longitude) 1.050334296178836E-02, // k = e cos(periapsis longitude) 3.135839828250385E+02 * PI / 180., // l = mean longitude (deg) 4.010968787918696E-03, // p = tan(i/2) sin(node) 8.667940544785077E-03, // q = tan(i/2) cos(node) 1.433400739278653E-05 * PI / 180., // apsis rate (deg/sec) 6.177086323833784E-03 * PI / 180., // mean motion (deg/sec) -1.430904354449910E-05 * PI / 180., // node rate (deg/sec) 2.680583409224650E+02 * PI / 180., // Laplacian plane pole ra (deg) 6.449442793849740E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 515, /* Adrastea */ 2443935.5492707, // element epoch Julian date 1.289799080000000E+05, // a = semi-major axis (km) 0.000000000000000E+00, // h = e sin(periapsis longitude) 0.000000000000000E+00, // k = e cos(periapsis longitude) 2.959936528085309E+02 * PI / 180., // l = mean longitude (deg) 0.000000000000000E+00, // p = tan(i/2) sin(node) 0.000000000000000E+00, // q = tan(i/2) cos(node) 9.727101238555177E-05 * PI / 180., // apsis rate (deg/sec) 1.396989439740030E-02 * PI / 180., // mean motion (deg/sec) -9.660246625153272E-05 * PI / 180., // node rate (deg/sec) 2.680583409224650E+02 * PI / 180., // Laplacian plane pole ra (deg) 6.449442793849740E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 516, /* Metis */ 2450392.938090278, // element epoch Julian date 1.279788680000000E+05, // a = semi-major axis (km) 0.000000000000000E+00, // h = e sin(periapsis longitude) 0.000000000000000E+00, // k = e cos(periapsis longitude) 8.488537675846435E+01 * PI / 180., // l = mean longitude (deg) 0.000000000000000E+00, // p = tan(i/2) sin(node) 0.000000000000000E+00, // q = tan(i/2) cos(node) 1.000079083666271E-04 * PI / 180., // apsis rate (deg/sec) 1.413489180972219E-02 * PI / 180., // mean motion (deg/sec) -9.930964475521999E-05 * PI / 180., // node rate (deg/sec) 2.680583409224650E+02 * PI / 180., // Laplacian plane pole ra (deg) 6.449442793849740E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 615, /* Atlas */ 2444786.50, // element epoch Julian date 1.376664620000000E+05, // a = semi-major axis (km) 0.000000000000000E+00, // h = e sin(periapsis longitude) 0.000000000000000E+00, // k = e cos(periapsis longitude) 1.865410967364615E+02 * PI / 180., // l = mean longitude (deg) 0.000000000000000E+00, // p = tan(i/2) sin(node) 0.000000000000000E+00, // q = tan(i/2) cos(node) 3.334584985561025E-05 * PI / 180., // apsis rate (deg/sec) 6.924845572071244E-03 * PI / 180., // mean motion (deg/sec) -3.318681015315788E-05 * PI / 180., // node rate (deg/sec) 4.058861887893110E+01 * PI / 180., // Laplacian plane pole ra (deg) 8.352533375484340E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 616, /* Prometheus */ 2444839.66820, // element epoch Julian date 1.393774260000000E+05, // a = semi-major axis (km) -1.242198628344100E-03, // h = e sin(periapsis longitude) -1.923542406971122E-03, // k = e cos(periapsis longitude) 1.885381493090702E+02 * PI / 180., // l = mean longitude (deg) 0.000000000000000E+00, // p = tan(i/2) sin(node) 0.000000000000000E+00, // q = tan(i/2) cos(node) 3.191385931437856E-05 * PI / 180., // apsis rate (deg/sec) 6.797328358741667E-03 * PI / 180., // mean motion (deg/sec) -3.176541076443906E-05 * PI / 180., // node rate (deg/sec) 4.058861887893110E+01 * PI / 180., // Laplacian plane pole ra (deg) 8.352533375484340E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 617, /* Pandora */ 2444839.66820, // element epoch Julian date 1.417126100000000E+05, // a = semi-major axis (km) 4.059003765493188E-03, // h = e sin(periapsis longitude) 1.621094455120198E-03, // k = e cos(periapsis longitude) 8.214726745964258E+01 * PI / 180., // l = mean longitude (deg) 0.000000000000000E+00, // p = tan(i/2) sin(node) 0.000000000000000E+00, // q = tan(i/2) cos(node) 3.008540659628256E-05 * PI / 180., // apsis rate (deg/sec) 6.629501272378168E-03 * PI / 180., // mean motion (deg/sec) -2.995009714624053E-05 * PI / 180., // node rate (deg/sec) 4.058861887893110E+01 * PI / 180., // Laplacian plane pole ra (deg) 8.352533375484340E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 618, /* Pan */ 2444839.294313820, // element epoch Julian date 1.335828000000000E+05, // a = semi-major axis (km) 0.000000000000000E+00, // h = e sin(periapsis longitude) 0.000000000000000E+00, // k = e cos(periapsis longitude) 1.220720210760172E+02 * PI / 180., // l = mean longitude (deg) 0.000000000000000E+00, // p = tan(i/2) sin(node) 0.000000000000000E+00, // q = tan(i/2) cos(node) 0.000000000000000E+00 * PI / 180., // apsis rate (deg/sec) 7.245879629629629E-03 * PI / 180., // mean motion (deg/sec) 0.000000000000000E+00 * PI / 180., // node rate (deg/sec) 4.058861887893110E+01 * PI / 180., // Laplacian plane pole ra (deg) 8.352533375484340E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 706, /* Cordelia */ 2446450.0, // element epoch Julian date 4.975172200000000E+04, // a = semi-major axis (km) 2.196155257809854E-05, // h = e sin(periapsis longitude) -2.616109646458082E-04, // k = e cos(periapsis longitude) 7.000654224200319E+01 * PI / 180., // l = mean longitude (deg) 4.593675494284073E-04, // p = tan(i/2) sin(node) 5.801117535036105E-04, // q = tan(i/2) cos(node) 1.739356380560684E-05 * PI / 180., // apsis rate (deg/sec) 1.243655458532564E-02 * PI / 180., // mean motion (deg/sec) -1.736934732358518E-05 * PI / 180., // node rate (deg/sec) 7.731126999999999E+01 * PI / 180., // Laplacian plane pole ra (deg) 1.517520000000000E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 707, /* Ophelia */ 2446450.0, // element epoch Julian date 5.376339000000000E+04, // a = semi-major axis (km) -3.133200269914386E-04, // h = e sin(periapsis longitude) -9.916862471808582E-03, // k = e cos(periapsis longitude) 2.980683571691387E+02 * PI / 180., // l = mean longitude (deg) 2.485085980620679E-04, // p = tan(i/2) sin(node) -8.694254731571202E-04, // q = tan(i/2) cos(node) 1.325233191309423E-05 * PI / 180., // apsis rate (deg/sec) 1.106977237274373E-02 * PI / 180., // mean motion (deg/sec) -1.323657157523269E-05 * PI / 180., // node rate (deg/sec) 7.731126999999999E+01 * PI / 180., // Laplacian plane pole ra (deg) 1.517520000000000E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 708, /* Bianca */ 2446450.0, // element epoch Julian date 5.916555000000000E+04, // a = semi-major axis (km) 8.971748009891265E-04, // h = e sin(periapsis longitude) -1.827531842136253E-04, // k = e cos(periapsis longitude) 2.399991076795366E+02 * PI / 180., // l = mean longitude (deg) 1.682254477338740E-03, // p = tan(i/2) sin(node) -9.465466268162635E-05, // q = tan(i/2) cos(node) 9.471202896165893E-06 * PI / 180., // apsis rate (deg/sec) 9.587823618966998E-03 * PI / 180., // mean motion (deg/sec) -9.462032392641446E-06 * PI / 180., // node rate (deg/sec) 7.731126999999999E+01 * PI / 180., // Laplacian plane pole ra (deg) 1.517520000000000E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 709, /* Cressida */ 2446450.0, // element epoch Julian date 6.176673000000000E+04, // a = semi-major axis (km) 2.135521241345163E-04, // h = e sin(periapsis longitude) -2.900698938540351E-04, // k = e cos(periapsis longitude) 1.743441300043773E+01 * PI / 180., // l = mean longitude (deg) 4.887836102488585E-05, // p = tan(i/2) sin(node) -8.094698786371765E-06, // q = tan(i/2) cos(node) 8.146067553304360E-06 * PI / 180., // apsis rate (deg/sec) 8.988222389105249E-03 * PI / 180., // mean motion (deg/sec) -8.138705872927823E-06 * PI / 180., // node rate (deg/sec) 7.731126999999999E+01 * PI / 180., // Laplacian plane pole ra (deg) 1.517520000000000E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 710, /* Desdemona */ 2446450.0, // element epoch Julian date 6.265836400000000E+04, // a = semi-major axis (km) 9.723239512740034E-05, // h = e sin(periapsis longitude) -7.979138390076314E-05, // k = e cos(periapsis longitude) 3.140004112591695E+02 * PI / 180., // l = mean longitude (deg) -7.934922021259821E-04, // p = tan(i/2) sin(node) 5.783817071912779E-04, // q = tan(i/2) cos(node) 7.747195444389255E-06 * PI / 180., // apsis rate (deg/sec) 8.796939111736320E-03 * PI / 180., // mean motion (deg/sec) -7.740437839303249E-06 * PI / 180., // node rate (deg/sec) 7.731126999999999E+01 * PI / 180., // Laplacian plane pole ra (deg) 1.517520000000000E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 711, /* Juliet */ 2446450.0, // element epoch Julian date 6.435822200000000E+04, // a = semi-major axis (km) 5.937201276777622E-04, // h = e sin(periapsis longitude) 2.899050111831932E-04, // k = e cos(periapsis longitude) 3.086703577437537E+02 * PI / 180., // l = mean longitude (deg) -1.968311842960469E-04, // p = tan(i/2) sin(node) -5.362697096802111E-04, // q = tan(i/2) cos(node) 7.054135446068334E-06 * PI / 180., // apsis rate (deg/sec) 8.450533968211433E-03 * PI / 180., // mean motion (deg/sec) -7.048279828285220E-06 * PI / 180., // node rate (deg/sec) 7.731126999999999E+01 * PI / 180., // Laplacian plane pole ra (deg) 1.517520000000000E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 712, /* Portia */ 2446450.0, // element epoch Julian date 6.609726500000000E+04, // a = semi-major axis (km) 4.439609698121046E-05, // h = e sin(periapsis longitude) -2.828284655050884E-05, // k = e cos(periapsis longitude) 3.408116975843583E+02 * PI / 180., // l = mean longitude (deg) -5.078474991455212E-04, // p = tan(i/2) sin(node) -8.893683431833462E-05, // q = tan(i/2) cos(node) 6.425626220921286E-06 * PI / 180., // apsis rate (deg/sec) 8.119056497467633E-03 * PI / 180., // mean motion (deg/sec) -6.420569130989476E-06 * PI / 180., // node rate (deg/sec) 7.731126999999999E+01 * PI / 180., // Laplacian plane pole ra (deg) 1.517520000000000E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 713, /* Rosalind */ 2446450.0, // element epoch Julian date 6.992679500000000E+04, // a = semi-major axis (km) 5.119640813135829E-05, // h = e sin(periapsis longitude) -1.018959466854515E-04, // k = e cos(periapsis longitude) 2.895039357650758E+02 * PI / 180., // l = mean longitude (deg) 5.408825347704433E-04, // p = tan(i/2) sin(node) 2.371745328519729E-03, // q = tan(i/2) cos(node) 5.276490952266129E-06 * PI / 180., // apsis rate (deg/sec) 7.461000211559000E-03 * PI / 180., // mean motion (deg/sec) -5.272963274877423E-06 * PI / 180., // node rate (deg/sec) 7.731126999999999E+01 * PI / 180., // Laplacian plane pole ra (deg) 1.517520000000000E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 714, /* Belinda */ 2446450.0, // element epoch Julian date 7.525561300000000E+04, // a = semi-major axis (km) -4.550504222119267E-05, // h = e sin(periapsis longitude) 5.770959489262290E-05, // k = e cos(periapsis longitude) 3.189675651318586E+02 * PI / 180., // l = mean longitude (deg) -2.637801571748877E-04, // p = tan(i/2) sin(node) 4.337154961306541E-05, // q = tan(i/2) cos(node) 4.082890883918376E-06 * PI / 180., // apsis rate (deg/sec) 6.682410746920550E-03 * PI / 180., // mean motion (deg/sec) -4.080413586943079E-06 * PI / 180., // node rate (deg/sec) 7.731126999999999E+01 * PI / 180., // Laplacian plane pole ra (deg) 1.517520000000000E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 715, /* Puck */ 2446450.0, // element epoch Julian date 8.600444400000000E+04, // a = semi-major axis (km) 1.180321117033704E-04, // h = e sin(periapsis longitude) 8.610832533193203E-06, // k = e cos(periapsis longitude) 3.316236003018918E+02 * PI / 180., // l = mean longitude (deg) -2.784974210178902E-03, // p = tan(i/2) sin(node) -6.156536661882880E-05, // q = tan(i/2) cos(node) 2.565678895181057E-06 * PI / 180., // apsis rate (deg/sec) 5.469266062887764E-03 * PI / 180., // mean motion (deg/sec) -2.564609522196947E-06 * PI / 180., // node rate (deg/sec) 7.731126999999999E+01 * PI / 180., // Laplacian plane pole ra (deg) 1.517520000000000E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 718, /* 1986 U10 */ 2446450.0, // element epoch Julian date 7.641700000000000E+04, // a = semi-major axis (km) 7.681292951425416E-04, // h = e sin(periapsis longitude) 7.861361182572046E-04, // k = e cos(periapsis longitude) 2.735200375536485E+01 * PI / 180., // l = mean longitude (deg) 0.000000000000000E+00, // p = tan(i/2) sin(node) 0.000000000000000E+00, // q = tan(i/2) cos(node) 3.865740741000000E-06 * PI / 180., // apsis rate (deg/sec) 6.530555555555556E-03 * PI / 180., // mean motion (deg/sec) -3.868332784035184E-06 * PI / 180., // node rate (deg/sec) 7.731126999999999E+01 * PI / 180., // Laplacian plane pole ra (deg) 1.517520000000000E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 808, /* Proteus, a.k.a. 1989N1 */ 2447757.0, // element epoch Julian date 1.176471083951680E+05, // a = semi-major axis (km) 4.312251344540417E-04, // h = e sin(periapsis longitude) -7.458925443725646E-05, // k = e cos(periapsis longitude) 2.139989422701393E+02 * PI / 180., // l = mean longitude (deg) 1.694800530014298E-04, // p = tan(i/2) sin(node) -2.975357731965572E-04, // q = tan(i/2) cos(node) 9.133585598101285E-07 * PI / 180., // apsis rate (deg/sec) 3.712562763632274E-03 * PI / 180., // mean motion (deg/sec) -8.950702955262804E-07 * PI / 180., // node rate (deg/sec) 2.991765842223010E+02 * PI / 180., // Laplacian plane pole ra (deg) 4.240406946265600E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 807, /* Larissa, a.k.a. 1989N2 */ 2447757.0, // element epoch Julian date 7.354833124009850E+04, // a = semi-major axis (km) 6.875683073658962E-04, // h = e sin(periapsis longitude) -1.202881276932585E-03, // k = e cos(periapsis longitude) 1.851599794084169E+02 * PI / 180., // l = mean longitude (deg) 3.142062025309209E-04, // p = tan(i/2) sin(node) 1.724283304388263E-03, // q = tan(i/2) cos(node) 4.550153261802315E-06 * PI / 180., // apsis rate (deg/sec) 7.512192673335218E-03 * PI / 180., // mean motion (deg/sec) -4.539452511249425E-06 * PI / 180., // node rate (deg/sec) 2.992637541641930E+02 * PI / 180., // Laplacian plane pole ra (deg) 4.289952994522980E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 805, /* Despina, a.k.a. 1989N3 */ 2447757.0, // element epoch Julian date 5.252594522658500E+04, // a = semi-major axis (km) 1.123331610676554E-04, // h = e sin(periapsis longitude) -8.163653609617338E-05, // k = e cos(periapsis longitude) 8.560400232763797E+01 * PI / 180., // l = mean longitude (deg) 2.422272552561876E-04, // p = tan(i/2) sin(node) -5.175089322884327E-04, // q = tan(i/2) cos(node) 1.477732764822587E-05 * PI / 180., // apsis rate (deg/sec) 1.245062680798036E-02 * PI / 180., // mean motion (deg/sec) -1.475538975056945E-05 * PI / 180., // node rate (deg/sec) 2.992707052843080E+02 * PI / 180., // Laplacian plane pole ra (deg) 4.293869732679440E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 806, /* Galatea, a.k.a. 1989N4 */ 2447757.0, // element epoch Julian date 6.195267233239820E+04, // a = semi-major axis (km) -7.697492391542343E-05, // h = e sin(periapsis longitude) -9.217446948395429E-05, // k = e cos(periapsis longitude) 4.697627474745159E+01 * PI / 180., // l = mean longitude (deg) 4.388943254654188E-04, // p = tan(i/2) sin(node) -1.817235120946147E-04, // q = tan(i/2) cos(node) 8.286590313969883E-06 * PI / 180., // apsis rate (deg/sec) 9.718284359168671E-03 * PI / 180., // mean motion (deg/sec) -8.273570255118341E-06 * PI / 180., // node rate (deg/sec) 2.992687097406640E+02 * PI / 180., // Laplacian plane pole ra (deg) 4.292745820252440E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 804, /* Thalassa, a.k.a. 1989N5 */ 2447757.0, // element epoch Julian date 5.007455197063890E+04, // a = semi-major axis (km) 1.129957809286051E-04, // h = e sin(periapsis longitude) 1.073198871736319E-04, // k = e cos(periapsis longitude) 2.400691948035834E+02 * PI / 180., // l = mean longitude (deg) 1.792401126562768E-03, // p = tan(i/2) sin(node) -5.702862029857829E-06, // q = tan(i/2) cos(node) 1.747529721870989E-05 * PI / 180., // apsis rate (deg/sec) 1.337680047670611E-02 * PI / 180., // mean motion (deg/sec) -1.744878934774183E-05 * PI / 180., // node rate (deg/sec) 2.992710300655420E+02 * PI / 180., // Laplacian plane pole ra (deg) 4.294052613950830E+01 * PI / 180. }, // Laplacian plane pole dec (deg) { 803, /* Naiad, a.k.a. 1989N6 */ 2447757.0, // element epoch Julian date 4.822730488030720E+04, // a = semi-major axis (km) 3.265904297054616E-04, // h = e sin(periapsis longitude) 2.487114121917455E-05, // k = e cos(periapsis longitude) 6.059250202880661E+01* PI / 180., // l = mean longitude (deg) 3.122033175813978E-02, // p = tan(i/2) sin(node) 2.714646565670971E-02, // q = tan(i/2) cos(node) 1.966548948527533E-05 * PI / 180., // apsis rate (deg/sec) 1.415328843692009E-02 * PI / 180., // mean motion (deg/sec) -1.983882754543970E-05 * PI / 180., // node rate (deg/sec) 2.992712348230880E+02 * PI / 180., // Laplacian plane pole ra (deg) 4.294167905377210E+01 * PI / 180. } // Laplacian plane pole dec (deg) }; #define SEC_PER_DAY 86400. /* In Win32, if this function is in a DLL, we gotta */ /* mention that fact up front. The following #ifdef */ /* lets us do that without confusing OSes that are */ /* blessedly ignorant of the weirdnesses of Windoze DLLs: */ #ifdef WIN32 #define DLL_FUNC __declspec( dllexport) #else #define DLL_FUNC #endif /* Given a JDE and a JPL ID number (see list at the top of this file), */ /* evaluate_rock( ) will compute the J2000 equatorial Cartesian */ /* position for that "rock" and will return 0. Otherwise, it returns */ /* -1 as an error condition. No other errors are returned... though */ /* hypothetically, something indicating you're outside the valid time */ /* coverage for the orbit in question would be nice. */ int DLL_FUNC evaluate_rock( const double jde, const int jpl_id, double *output_vect) { int i; const ROCK *rptr = rocks; for( i = N_ROCKS; i; i--, rptr++) if( rptr->jpl_id == jpl_id) { double dt = jde - rptr->epoch_jd; double mean_lon = rptr->mean_lon0 + dt * SEC_PER_DAY * rptr->mean_motion; double avect[3], bvect[3], cvect[3]; double h, k, p, q, tsin, tcos, r, e, omega, true_lon; double a_fraction, b_fraction, c_fraction, dot_prod; /* avect is at right angles to Laplacian pole, */ /* but in plane of the J2000 equator: */ avect[0] = -sin( rptr->laplacian_pole_ra); avect[1] = cos( rptr->laplacian_pole_ra); avect[2] = 0.; /* bvect is at right angles to Laplacian pole */ /* _and_ to avect: */ tsin = sin( rptr->laplacian_pole_dec); tcos = cos( rptr->laplacian_pole_dec); bvect[0] = -avect[1] * tsin; bvect[1] = avect[0] * tsin; bvect[2] = tcos; /* cvect is the Laplacian pole vector: */ cvect[0] = avect[1] * tcos; cvect[1] = -avect[0] * tcos; cvect[2] = tsin; /* Rotate the (h, k) vector to account for */ /* a constant apsidal motion: */ tsin = sin( dt * rptr->apsis_rate * SEC_PER_DAY); tcos = cos( dt * rptr->apsis_rate * SEC_PER_DAY); h = rptr->k * tsin + rptr->h * tcos; k = rptr->k * tcos - rptr->h * tsin; /* I'm sure there's a better way to do this... */ /* all I do here is to compute the eccentricity */ /* and omega, a.k.a. longitude of perihelion, */ /* and do a first-order correction to get the */ /* 'actual' r and true longitude values. */ e = sqrt( h * h + k * k); omega = atan2( h, k); true_lon = mean_lon + 2. * e * sin( mean_lon - omega) + 1.25 * e * e * sin( 2. * (mean_lon - omega)); r = rptr->a * (1. - e * e) / (1 + e * cos( true_lon - omega)); /* Just as we rotated (h,k), we gotta rotate */ /* the (p,q) vector to account for precession */ /* in the Laplacian plane: */ tsin = sin( dt * rptr->node_rate * SEC_PER_DAY); tcos = cos( dt * rptr->node_rate * SEC_PER_DAY); p = rptr->q * tsin + rptr->p * tcos; q = rptr->q * tcos - rptr->p * tsin; /* Now we evaluate the position in components */ /* along avect, bvect, cvect. I derived the */ /* formulae from scratch... sorry I can't */ /* give references: */ tsin = sin( true_lon); tcos = cos( true_lon); dot_prod = 2. * (q * tsin - p * tcos) / (1. + p * p + q * q); a_fraction = tcos + p * dot_prod; b_fraction = tsin - q * dot_prod; c_fraction = dot_prod; /* Now that we've got components on each axis, */ /* the remainder is trivial: */ for( i = 0; i < 3; i++) output_vect[i] = r * (a_fraction * avect[i] + b_fraction * bvect[i] + c_fraction * cvect[i]); return( 0); } return( -1); } #ifdef TEST_PROGRAM /* A simple piece of test code that, given a JPL ID and a JD, prints out */ /* the result of evaluate_rock. Comparison to Horizons is straightforward, */ /* and indicates agreement to better than a meter. */ #include #include #include void main( int argc, char **argv) { double vect[3], r2 = 0.; int i, jpl_id = atoi( argv[2]); evaluate_rock( atof( argv[1]), jpl_id, vect); printf( "%lf %lf %lf\n", vect[0], vect[1], vect[2]); for( i = 0; i < 3; i++) r2 += vect[i] * vect[i]; printf( "\n%lf", sqrt( r2)); } #endif