← Back to branch summary

~ubuntu-branches/ubuntu/warty/xplanet/warty

« back to all changes in this revision

Viewing changes to src/libplanet/libmoons/neptune.cpp

  • Committer: Bazaar Package Importer
  • Author(s): LaMont Jones
  • Date: 2004-08-24 07:14:00 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20040824071400-2dr4qnjbjmm8z3ia
Tags: 1.0.6-1ubuntu1
Build-depend: libtiff4-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/libplanet/libmoons/neptune.cpp
 
1
#include <cmath>
 
2
#include <cstdio>
 
3
#include <cstdlib>
 
4
using namespace std;
 
5
 
 
6
/*
 
7
  Ephemerides for Triton and Nereid are described in Jacobson, 
 
8
  Astron. Astrophys. 231, 241-250 (1990)
 
9
*/
 
10
 
 
11
#include "body.h"
 
12
#include "xpUtil.h"
 
13
 
 
14
void
 
15
nepsat(const double jd, const body b, double &X, double &Y, double &Z)
 
16
{
 
17
    double td;       // Julian days from reference date
 
18
    double ty;       // Julian years from reference date
 
19
    double tc;       // Julian centuries from reference date
 
20
 
 
21
    double a;        // semimajor axis
 
22
    double L;        // mean longitude
 
23
    double e;        // eccentricity
 
24
    double w;        // longitude of periapse
 
25
    double i;        // inclination of orbit
 
26
    double o;        // longitude of ascending node
 
27
 
 
28
    double ma;       // mean anomaly
 
29
 
 
30
    double N;        // node of the orbital reference plane on the
 
31
                     // Earth equator B1950
 
32
    double J;        // inclination of orbital reference plane with
 
33
                     // respect to the Earth equator B1950
 
34
 
 
35
    switch (b)
 
36
    {
 
37
    case TRITON:
 
38
        td = jd - 2433282.5;
 
39
        ty = td/365.25;
 
40
        tc = ty/100;
 
41
 
 
42
        a = 354611.773;
 
43
        L = (49.85334766 + 61.25726751 * td) * deg_to_rad;
 
44
        e = 0.0004102259410;
 
45
        i = 157.6852321 * deg_to_rad;
 
46
        o = (151.7973992 + 0.5430763965 * ty) * deg_to_rad;
 
47
 
 
48
        w = (236.7318362 + 0.5295275852 * ty) * deg_to_rad;
 
49
 
 
50
        ma = L - w;
 
51
 
 
52
        w += o;
 
53
 
 
54
        // inclination and node of the invariable plane on the Earth
 
55
        // equator of 1950
 
56
        J = (90 - 42.51071244) * deg_to_rad;
 
57
        N = (90 + 298.3065940) * deg_to_rad;
 
58
    
 
59
        break;
 
60
    case NEREID:
 
61
        td = jd - 2433680.5;
 
62
        tc = td/36525;
 
63
 
 
64
        a = 5511233.255;
 
65
        L = (251.14984688 + 0.9996465329 * td) * deg_to_rad;
 
66
        e = 0.750876291;
 
67
        i = 6.748231850 * deg_to_rad;
 
68
        o = (315.9958928 - 3.650272562 * tc) * deg_to_rad;
 
69
 
 
70
        w = (251.7242240 + 0.8696048083 * tc) * deg_to_rad;
 
71
 
 
72
        ma = L - w;
 
73
 
 
74
        w -= o;
 
75
 
 
76
        // inclination and node of Neptune's orbit on the Earth
 
77
        // equator of 1950
 
78
        J = 22.313 * deg_to_rad;
 
79
        N = 3.522 * deg_to_rad;
 
80
        break;
 
81
    default:
 
82
        xpExit("Unknown Neptune satellite\n", __FILE__, __LINE__);
 
83
    }
 
84
 
 
85
    double E = kepler(e, ma);
 
86
 
 
87
    // convert semi major axis from km to AU
 
88
    a /= AU_to_km;
 
89
 
 
90
    // rectangular coordinates on the orbit plane, x-axis is toward
 
91
    // pericenter
 
92
    X = a * (cos(E) - e);
 
93
    Y = a * sqrt(1 - e*e) * sin(E);
 
94
    Z = 0;
 
95
 
 
96
    // rotate towards ascending node of the orbit
 
97
    rotateZ(X, Y, Z, -w);
 
98
 
 
99
    // rotate towards orbital reference plane
 
100
    rotateX(X, Y, Z, -i);
 
101
 
 
102
    // rotate towards ascending node of the orbital reference plane on
 
103
    // the Earth equator B1950
 
104
    rotateZ(X, Y, Z, -o);
 
105
 
 
106
    // rotate towards Earth equator B1950
 
107
    rotateX(X, Y, Z, -J);
 
108
 
 
109
    // rotate to vernal equinox
 
110
    rotateZ(X, Y, Z, -N);
 
111
 
 
112
    // precess to J2000
 
113
    precessB1950J2000(X, Y, Z);
 
114
}