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

« back to all changes in this revision

Viewing changes to moonpos.cc

  • 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:
1
 
/****************************************************************************
2
 
    Xplanet 0.94 - render an image of a planet into an X window
3
 
    Copyright (C) 2002 Hari Nair <hari@alumni.caltech.edu>
4
 
 
5
 
    This program is free software; you can redistribute it and/or modify
6
 
    it under the terms of the GNU General Public License as published by
7
 
    the Free Software Foundation; either version 2 of the License, or
8
 
    (at your option) any later version.
9
 
 
10
 
    This program is distributed in the hope that it will be useful,
11
 
    but WITHOUT ANY WARRANTY; without even the implied warranty of
12
 
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13
 
    GNU General Public License for more details.
14
 
 
15
 
    You should have received a copy of the GNU General Public License
16
 
    along with this program; if not, write to the Free Software
17
 
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
18
 
****************************************************************************/
19
 
 
20
 
/* Based on Chapter 30 of Astronomical Formulae for Calculators by Meeus */
21
 
 
22
 
#include <cmath>
23
 
using namespace std;
24
 
 
25
 
#include "util.h"
26
 
 
27
 
static double S(const double T)
28
 
{
29
 
    return sin(T*deg_to_rad);
30
 
}
31
 
 
32
 
static double C(const double T)
33
 
{
34
 
    return cos(T*deg_to_rad);
35
 
}
36
 
 
37
 
void 
38
 
moonpos(double &moonlon, double &moonlat, double &moondist, const double T)
39
 
{
40
 
    double M  = poly(358.475833, 35999.04975, -1.500E-4, -3.30E-6, T);
41
 
    double Lp = poly(270.434164, 481267.8831, -1.133E-3,  1.90E-6, T);
42
 
    double Mp = poly(296.104608, 477198.8491,  9.192E-3,  1.44E-5, T);
43
 
    double D  = poly(350.737486, 445267.1142, -1.436E-3,  1.90E-6, T);
44
 
    double F  = poly( 11.250889, 483202.0251, -3.211E-3, -3.00E-7, T);
45
 
    double omega = poly(259.183275,-1934.142, 2.078E-3, 2.2E-6, T);
46
 
    Lp = Lp + (  (2.330E-4) * S(51.2 + 20.2*T) 
47
 
                 + (1.964E-3) * S(omega)
48
 
                 + (3.964E-3) * S(poly (346.56, 132.87, -9.1731E-3, 0, T)));
49
 
    M = M - (1.778E-3) * S(51.2 + 20.2*T);
50
 
    Mp = Mp + (  (8.170E-4) * S(51.2 + 20.2*T) 
51
 
                 + (2.541E-3) * S(omega)
52
 
                 + (3.964E-3) * S(poly (346.56, 132.87, -9.1731E-3, 0, T)));
53
 
    D = D + (  (2.011E-3) * S(51.2 + 20.2*T) 
54
 
               + (1.964E-3) * S(omega)
55
 
               + (3.964E-3) * S(poly (346.56, 132.87, -9.1731E-3, 0, T)));
56
 
    F = F - (  (2.4691E-2) * S(omega)
57
 
               - (4.3280E-3) * S(omega + 275.05 - 2.3*T)
58
 
               + (3.9640E-3) * S(poly (346.56, 132.87, -9.1731E-3, 0, T)));
59
 
    double E = 1 - T * (2.495E-3 + 7.52E-6 * T);
60
 
    double E2 = E*E;
61
 
    double lambda = (Lp 
62
 
                     + 6.28875 * S(Mp)
63
 
                     + 1.274018 * S(2*D - Mp) + .658309 * S(2*D)
64
 
                     + .213616 * S(2*Mp) - E * .185596 * S(M)
65
 
                     - .114336*S(2*F) + .58793E-1*S(2*(D - Mp))
66
 
                     + E*.57212E-1*S(2*D - M - Mp)
67
 
                     + .5332E-1*S(2*D + Mp)
68
 
                     + E*.45874E-1*S(2*D - M)
69
 
                     + E*.41024E-1*S(Mp - M) - .34718E-1*S(D)
70
 
                     - E*.30465E-1*S(M + Mp) + .15326E-1*S(2*(D - F))
71
 
                     - .12528E-1*S(2*F - Mp) - .1098E-1*S(2*F - Mp)
72
 
                     + .10674E-1*S(4*D - Mp) + .10034E-1*S(3*Mp)
73
 
                     + .8548E-2*S(4*D - 2*Mp)
74
 
                     - E*.791E-2*S(M - Mp + 2*D)
75
 
                     - E*.6783E-2*S(2*D + M) + .5162E-2*S(Mp - D)
76
 
                     + E*.5E-2*S(M + D)
77
 
                     + E*.4049E-2*S(Mp - M + 2*D)
78
 
                     + .3996E-2*S(2*(Mp + D)) + .3862E-2*S(4*D)
79
 
                     + .3665E-2*S(2*D - 3*Mp)
80
 
                     + E*.2695E-2*S(2*Mp - M)
81
 
                     + .2602E-2*S(Mp - 2*F - 2*D)
82
 
                     + E*.2396E-2*S(2*D - M - 2*Mp) - .2349E-2*S(Mp+D)
83
 
                     + E2 * .2249E-2*S(2*(D - M))
84
 
                     - E*.2125E-2*S(2*Mp + M) - E2 * .2079E-2*S(2*M)
85
 
                     + E2 * .2059E-2*S(2*D - Mp - 2*M)
86
 
                     - .1773E-2*S(Mp + 2*D - 2*F)
87
 
                     - .1595E-2*S(2*(F + D))
88
 
                     + E*.122E-2*S(4*D - M - Mp)
89
 
                     - .111E-2*S(2*(Mp + F)) + .892E-3*S(Mp - 3*D)
90
 
                     - E*.811E-3*S(M + Mp + 2*D)
91
 
                     + E*.761E-3*S(4*D - M - 2*Mp)
92
 
                     + E2 * .717E-3*S(Mp - 2*M)
93
 
                     + E2 * .704E-3*S(Mp - 2*M - 2*D)
94
 
                     + E*.693E-3*S(M - 2*Mp + 2*D)
95
 
                     + E*.598E-3*S(2*D - M - 2*F) + .55E-3*S(Mp + 4*D)
96
 
                     + .538E-3*S(4*Mp) + E*.521E-3*S(4*D - M)
97
 
                     + .486E-3*S(2*Mp - D));
98
 
    double B = (5.128189*S(F) + .280606*S(Mp + F)
99
 
                + .277693*S(Mp - F) + .173238*S(2*D - F)
100
 
                + .55413E-1*S(2*D + F - Mp)
101
 
                + .46272E-1*S(2*D - F - Mp)
102
 
                + .32573E-1*S(2*D + F) + .17198E-1*S(2*Mp + F)
103
 
                + .9267E-2*S(2*D + Mp - F) + .8823E-2*S(2*Mp - F)
104
 
                + E*.8247E-2*S(2*D - M - F)
105
 
                + .4323E-2*S(2*D - F - 2*Mp)
106
 
                + .42E-2*S(2*D + F + Mp)
107
 
                + E*.3372E-2*S(F - M - 2*D)
108
 
                + E*.2472E-2*S(2*D + F - M - Mp)
109
 
                + E*.2222E-2*S(2*D + F - M)
110
 
                + E*.2072E-2*S(2*D - F - M - Mp)
111
 
                + E*.1877E-2*S(F - M + Mp)
112
 
                + .1828E-2*S(4*D - F - Mp)
113
 
                - E*.1803E-2*S(F + M) - .175E-2*S(3*F)
114
 
                + E*.157E-2*S(Mp - M - F) - .1487E-2*S(F + D)
115
 
                - E*.1481E-2*S(F + M + Mp)
116
 
                + E*.1417E-2*S(F - M - Mp)
117
 
                + E*.135E-2*S(F - M) + .133E-2*S(F - D)
118
 
                + .1106E-2*S(F + 3*Mp) + .102E-2*S(4*D - F)
119
 
                + .833E-3*S(F + 4*D - Mp) + .781E-3*S(Mp - 3*F)
120
 
                + .67E-3*S(F + 4*D - 2*Mp) + .606E-3*S(2*D - 3*F)
121
 
                + .597E-3*S(2*D + 2*Mp - F)
122
 
                + E *.492E-3*S(2*D + Mp - M - F)
123
 
                + .45E-3*S(2*Mp - F - 2*D) + .439E-3*S(3*Mp - F)
124
 
                + .423E-3*S(F + 2*D + 2*Mp)
125
 
                + .422E-3*S(2*D - F - 3*Mp)
126
 
                - E*.367E-3*S(M + F + 2*D - Mp)
127
 
                - E*.353E-3*S(M + F + 2*D) + .331E-3*S(F + 4*D)
128
 
                + E*.317E-3*S(2*D + F - M + Mp)
129
 
                + E2 * .306E-3*S(2*D - 2*M - F)
130
 
                - .283E-3*S(Mp + 3*F));
131
 
    double omega1 = .4664E-3*C(omega);
132
 
    double omega2 = .754E-4*C(omega + 275.05 - 2.3*T);
133
 
    double beta = B*(1 - omega1 - omega2);
134
 
    double parallax = (.9507234 + .51818E-1*C(Mp)
135
 
                       + .9531E-2*C(2*D - Mp) + .7843E-2*C(2*D)
136
 
                       + .2824E-2*C(2*Mp) + .857E-3*C(2*D + Mp)
137
 
                       + E*.533E-3*C(2*D - M)
138
 
                       + E*.401E-3*C(2*D - M - Mp)
139
 
                       + E*.32E-3*C(Mp - M) - .271E-3*C(D)
140
 
                       - .264E-3*C(M + Mp) - .198E-3*C(2*F - Mp)
141
 
                       + .173E-3*C(3*Mp) + .167E-3*C(4*D - Mp)
142
 
                       - E*.111E-3*C(M) + .103E-3*C(4*D - 2*Mp)
143
 
                       - .84E-4*C(2*(Mp - D)) - E*.83E-4*C(2*D + M)
144
 
                       + .79E-4*C(2*(D + Mp)) + .72E-4*C(4*D)
145
 
                       + E*.64E-4*C(2*D - M + Mp)
146
 
                       - E*.63E-4*C(2*D + M - Mp)
147
 
                       + E*.41E-4*C(M + D) + E*.35E-4*C(2*Mp - M)
148
 
                       - .33E-4*C(3*Mp - 2*D) - .3E-4*C(Mp + D)
149
 
                       - .29E-4*C(2*(F - D)) - E*.29E-4*C(2*Mp + M)
150
 
                       + E2 * .26E-4*C(2*(D - M))
151
 
                       - .23E-4* C(2*F - 2*D + Mp)
152
 
                       + E*.19E-4*C(4*D - M - Mp));
153
 
 
154
 
    moonlon = lambda * deg_to_rad;
155
 
    moonlat = beta * deg_to_rad;
156
 
    moondist = 1/S(parallax);
157
 
}