1
by acapnotic
Initial revision |
1 |
#!/usr/bin/env python
|
2 |
# moon.py, based on code by John Walker (http://www.fourmilab.ch/)
|
|
3 |
# ported to Python by Kevin Turner <acapnotic@twistedmatrix.com>
|
|
4 |
# on June 6, 2001 (JDN 2452066.52491), under a full moon.
|
|
5 |
#
|
|
6 |
# This program is in the public domain: "Do what thou wilt shall be
|
|
7 |
# the whole of the law".
|
|
8 |
||
9 |
"""Functions to find the phase of the moon.
|
|
10 |
||
11 |
Ported from \"A Moon for the Sun\" (aka moontool.c), a program by the
|
|
12 |
venerable John Walker. He used algoritms from \"Practical Astronomy
|
|
13 |
With Your Calculator\" by Peter Duffett-Smith, Second Edition.
|
|
14 |
||
15 |
For the full history of the code, as well as references to other
|
|
16 |
reading material and other entertainments, please refer to John
|
|
17 |
Walker's website,
|
|
18 |
http://www.fourmilab.ch/
|
|
19 |
(Look under the Science/Astronomy and Space heading.)
|
|
20 |
||
21 |
The functions of primary interest provided by this module are phase(),
|
|
22 |
which gives you a variety of data on the status of the moon for a
|
|
23 |
given date; and phase_hunt(), which given a date, finds the dates of
|
|
24 |
the nearest full moon, new moon, etc.
|
|
25 |
"""
|
|
26 |
||
3.1.1
by Kevin Turner
remove import *, comment out calculations that are not used for phase calculation. |
27 |
from math import sin, cos, floor, sqrt, pi, tan, atan # asin, atan2 |
1
by acapnotic
Initial revision |
28 |
import bisect |
29 |
try: |
|
30 |
import DateTime |
|
31 |
except ImportError: |
|
32 |
from mx import DateTime |
|
33 |
||
34 |
__TODO__ = [ |
|
35 |
'Add command-line interface.', |
|
36 |
'Make front-end modules for ASCII and various GUIs.', |
|
37 |
]
|
|
38 |
||
39 |
# Precision used when describing the moon's phase in textual format,
|
|
40 |
# in phase_string().
|
|
41 |
PRECISION = 0.05 |
|
42 |
NEW = 0 / 4.0 |
|
43 |
FIRST = 1 / 4.0 |
|
44 |
FULL = 2 / 4.0 |
|
45 |
LAST = 3 / 4.0 |
|
46 |
NEXTNEW = 4 / 4.0 |
|
47 |
||
48 |
||
49 |
class MoonPhase: |
|
50 |
"""I describe the phase of the moon.
|
|
51 |
||
52 |
I have the following properties:
|
|
53 |
date - a DateTime instance
|
|
54 |
phase - my phase, in the range 0.0 .. 1.0
|
|
55 |
phase_text - a string describing my phase
|
|
56 |
illuminated - the percentage of the face of the moon illuminated
|
|
57 |
angular_diameter - as seen from Earth, in degrees.
|
|
58 |
sun_angular_diameter - as seen from Earth, in degrees.
|
|
59 |
||
60 |
new_date - the date of the most recent new moon
|
|
61 |
q1_date - the date the moon reaches 1st quarter in this cycle
|
|
62 |
full_date - the date of the full moon in this cycle
|
|
63 |
q3_date - the date the moon reaches 3rd quarter in this cycle
|
|
64 |
nextnew_date - the date of the next new moon
|
|
65 |
"""
|
|
66 |
||
67 |
def __init__(self, date=DateTime.now()): |
|
68 |
"""MoonPhase constructor.
|
|
69 |
||
70 |
Give me a date, as either a Julian Day Number or a DateTime
|
|
71 |
object."""
|
|
72 |
||
73 |
if not isinstance(date, DateTime.DateTimeType): |
|
74 |
self.date = DateTime.DateTimeFromJDN(date) |
|
75 |
else: |
|
76 |
self.date = date |
|
77 |
||
78 |
self.__dict__.update(phase(self.date)) |
|
79 |
||
80 |
self.phase_text = phase_string(self.phase) |
|
81 |
||
82 |
def __getattr__(self, a): |
|
83 |
if a in ['new_date', 'q1_date', 'full_date', 'q3_date', |
|
84 |
'nextnew_date']: |
|
85 |
||
86 |
(self.new_date, self.q1_date, self.full_date, |
|
87 |
self.q3_date, self.nextnew_date) = phase_hunt(self.date) |
|
88 |
||
89 |
return getattr(self,a) |
|
6
by Kevin Turner
change exception and map usage to be python3 compatible |
90 |
raise AttributeError(a) |
1
by acapnotic
Initial revision |
91 |
|
92 |
def __repr__(self): |
|
5
by Kevin Turner
remove type imports. |
93 |
if type(self.date) is int: |
1
by acapnotic
Initial revision |
94 |
jdn = self.date |
95 |
else: |
|
96 |
jdn = self.date.jdn |
|
97 |
||
98 |
return "<%s(%d)>" % (self.__class__, jdn) |
|
99 |
||
100 |
def __str__(self): |
|
5
by Kevin Turner
remove type imports. |
101 |
if type(self.date) is int: |
1
by acapnotic
Initial revision |
102 |
d = DateTime.DateTimeFromJDN(self.date) |
103 |
else: |
|
104 |
d = self.date |
|
105 |
s = "%s for %s, %s (%%%.2f illuminated)" %\ |
|
106 |
(self.__class__, d.strftime(), self.phase_text, |
|
107 |
self.illuminated * 100) |
|
108 |
||
109 |
return s |
|
110 |
||
111 |
||
112 |
class AstronomicalConstants: |
|
113 |
||
114 |
# JDN stands for Julian Day Number
|
|
115 |
# Angles here are in degrees
|
|
116 |
||
117 |
# 1980 January 0.0 in JDN
|
|
118 |
# XXX: DateTime(1980).jdn yields 2444239.5 -- which one is right?
|
|
119 |
epoch = 2444238.5 |
|
120 |
||
121 |
# Ecliptic longitude of the Sun at epoch 1980.0
|
|
122 |
ecliptic_longitude_epoch = 278.833540 |
|
123 |
||
124 |
# Ecliptic longitude of the Sun at perigee
|
|
125 |
ecliptic_longitude_perigee = 282.596403 |
|
126 |
||
127 |
# Eccentricity of Earth's orbit
|
|
128 |
eccentricity = 0.016718 |
|
129 |
||
130 |
# Semi-major axis of Earth's orbit, in kilometers
|
|
131 |
sun_smaxis = 1.49585e8 |
|
132 |
||
133 |
# Sun's angular size, in degrees, at semi-major axis distance
|
|
134 |
sun_angular_size_smaxis = 0.533128 |
|
135 |
||
136 |
## Elements of the Moon's orbit, epoch 1980.0
|
|
137 |
||
138 |
# Moon's mean longitude at the epoch
|
|
139 |
moon_mean_longitude_epoch = 64.975464 |
|
140 |
# Mean longitude of the perigee at the epoch
|
|
141 |
moon_mean_perigee_epoch = 349.383063 |
|
142 |
||
143 |
# Mean longitude of the node at the epoch
|
|
144 |
node_mean_longitude_epoch = 151.950429 |
|
145 |
||
146 |
# Inclination of the Moon's orbit
|
|
147 |
moon_inclination = 5.145396 |
|
148 |
||
149 |
# Eccentricity of the Moon's orbit
|
|
150 |
moon_eccentricity = 0.054900 |
|
151 |
||
152 |
# Moon's angular size at distance a from Earth
|
|
153 |
moon_angular_size = 0.5181 |
|
154 |
||
155 |
# Semi-mojor axis of the Moon's orbit, in kilometers
|
|
156 |
moon_smaxis = 384401.0 |
|
157 |
# Parallax at a distance a from Earth
|
|
158 |
moon_parallax = 0.9507 |
|
159 |
||
160 |
# Synodic month (new Moon to new Moon), in days
|
|
161 |
synodic_month = 29.53058868 |
|
162 |
||
163 |
# Base date for E. W. Brown's numbered series of lunations (1923 January 16)
|
|
164 |
lunations_base = 2423436.0 |
|
165 |
||
166 |
## Properties of the Earth
|
|
167 |
||
168 |
earth_radius = 6378.16 |
|
169 |
||
170 |
c = AstronomicalConstants() |
|
171 |
||
172 |
# Little handy mathematical functions.
|
|
173 |
||
174 |
fixangle = lambda a: a - 360.0 * floor(a/360.0) |
|
175 |
torad = lambda d: d * pi / 180.0 |
|
176 |
todeg = lambda r: r * 180.0 / pi |
|
177 |
dsin = lambda d: sin(torad(d)) |
|
178 |
dcos = lambda d: cos(torad(d)) |
|
179 |
||
180 |
def phase_string(p): |
|
181 |
phase_strings = ( |
|
182 |
(NEW + PRECISION, "new"), |
|
183 |
(FIRST - PRECISION, "waxing crescent"), |
|
184 |
(FIRST + PRECISION, "first quarter"), |
|
185 |
(FULL - PRECISION, "waxing gibbous"), |
|
186 |
(FULL + PRECISION, "full"), |
|
187 |
(LAST - PRECISION, "waning gibbous"), |
|
188 |
(LAST + PRECISION, "last quarter"), |
|
189 |
(NEXTNEW - PRECISION, "waning crescent"), |
|
190 |
(NEXTNEW + PRECISION, "new")) |
|
191 |
||
6
by Kevin Turner
change exception and map usage to be python3 compatible |
192 |
i = bisect.bisect([a[0] for a in phase_strings], p) |
1
by acapnotic
Initial revision |
193 |
|
194 |
return phase_strings[i][1] |
|
195 |
||
196 |
||
197 |
def phase(phase_date=DateTime.now()): |
|
198 |
"""Calculate phase of moon as a fraction:
|
|
199 |
||
200 |
The argument is the time for which the phase is requested,
|
|
201 |
expressed in either a DateTime or by Julian Day Number.
|
|
202 |
||
203 |
Returns a dictionary containing the terminator phase angle as a
|
|
204 |
percentage of a full circle (i.e., 0 to 1), the illuminated
|
|
205 |
fraction of the Moon's disc, the Moon's age in days and fraction,
|
|
206 |
the distance of the Moon from the centre of the Earth, and the
|
|
207 |
angular diameter subtended by the Moon as seen by an observer at
|
|
208 |
the centre of the Earth."""
|
|
209 |
||
210 |
# Calculation of the Sun's position
|
|
211 |
||
212 |
# date within the epoch
|
|
213 |
if hasattr(phase_date, "jdn"): |
|
214 |
day = phase_date.jdn - c.epoch |
|
215 |
else: |
|
216 |
day = phase_date - c.epoch |
|
217 |
||
218 |
# Mean anomaly of the Sun
|
|
219 |
N = fixangle((360/365.2422) * day) |
|
220 |
# Convert from perigee coordinates to epoch 1980
|
|
221 |
M = fixangle(N + c.ecliptic_longitude_epoch - c.ecliptic_longitude_perigee) |
|
222 |
||
223 |
# Solve Kepler's equation
|
|
224 |
Ec = kepler(M, c.eccentricity) |
|
225 |
Ec = sqrt((1 + c.eccentricity) / (1 - c.eccentricity)) * tan(Ec/2.0) |
|
226 |
# True anomaly
|
|
227 |
Ec = 2 * todeg(atan(Ec)) |
|
228 |
# Suns's geometric ecliptic longuitude
|
|
229 |
lambda_sun = fixangle(Ec + c.ecliptic_longitude_perigee) |
|
230 |
||
231 |
# Orbital distance factor
|
|
232 |
F = ((1 + c.eccentricity * cos(torad(Ec))) / (1 - c.eccentricity**2)) |
|
233 |
||
234 |
# Distance to Sun in km
|
|
235 |
sun_dist = c.sun_smaxis / F |
|
236 |
sun_angular_diameter = F * c.sun_angular_size_smaxis |
|
237 |
||
238 |
########
|
|
239 |
#
|
|
240 |
# Calculation of the Moon's position
|
|
241 |
||
242 |
# Moon's mean longitude
|
|
243 |
moon_longitude = fixangle(13.1763966 * day + c.moon_mean_longitude_epoch) |
|
244 |
||
245 |
# Moon's mean anomaly
|
|
246 |
MM = fixangle(moon_longitude - 0.1114041 * day - c.moon_mean_perigee_epoch) |
|
247 |
||
248 |
# Moon's ascending node mean longitude
|
|
3.1.1
by Kevin Turner
remove import *, comment out calculations that are not used for phase calculation. |
249 |
# MN = fixangle(c.node_mean_longitude_epoch - 0.0529539 * day)
|
1
by acapnotic
Initial revision |
250 |
|
251 |
evection = 1.2739 * sin(torad(2*(moon_longitude - lambda_sun) - MM)) |
|
252 |
||
253 |
# Annual equation
|
|
254 |
annual_eq = 0.1858 * sin(torad(M)) |
|
255 |
||
256 |
# Correction term
|
|
257 |
A3 = 0.37 * sin(torad(M)) |
|
258 |
||
259 |
MmP = MM + evection - annual_eq - A3 |
|
260 |
||
261 |
# Correction for the equation of the centre
|
|
262 |
mEc = 6.2886 * sin(torad(MmP)) |
|
263 |
||
264 |
# Another correction term
|
|
265 |
A4 = 0.214 * sin(torad(2 * MmP)) |
|
266 |
||
267 |
# Corrected longitude
|
|
268 |
lP = moon_longitude + evection + mEc - annual_eq + A4 |
|
269 |
||
270 |
# Variation
|
|
271 |
variation = 0.6583 * sin(torad(2*(lP - lambda_sun))) |
|
272 |
||
273 |
# True longitude
|
|
274 |
lPP = lP + variation |
|
275 |
||
3.1.1
by Kevin Turner
remove import *, comment out calculations that are not used for phase calculation. |
276 |
#
|
277 |
# Calculation of the Moon's inclination
|
|
278 |
# unused for phase calculation.
|
|
279 |
||
1
by acapnotic
Initial revision |
280 |
# Corrected longitude of the node
|
3.1.1
by Kevin Turner
remove import *, comment out calculations that are not used for phase calculation. |
281 |
# NP = MN - 0.16 * sin(torad(M))
|
1
by acapnotic
Initial revision |
282 |
|
283 |
# Y inclination coordinate
|
|
3.1.1
by Kevin Turner
remove import *, comment out calculations that are not used for phase calculation. |
284 |
# y = sin(torad(lPP - NP)) * cos(torad(c.moon_inclination))
|
1
by acapnotic
Initial revision |
285 |
|
286 |
# X inclination coordinate
|
|
3.1.1
by Kevin Turner
remove import *, comment out calculations that are not used for phase calculation. |
287 |
# x = cos(torad(lPP - NP))
|
1
by acapnotic
Initial revision |
288 |
|
289 |
# Ecliptic longitude (unused?)
|
|
3.1.1
by Kevin Turner
remove import *, comment out calculations that are not used for phase calculation. |
290 |
# lambda_moon = todeg(atan2(y,x)) + NP
|
1
by acapnotic
Initial revision |
291 |
|
292 |
# Ecliptic latitude (unused?)
|
|
3.1.1
by Kevin Turner
remove import *, comment out calculations that are not used for phase calculation. |
293 |
# BetaM = todeg(asin(sin(torad(lPP - NP)) * sin(torad(c.moon_inclination))))
|
1
by acapnotic
Initial revision |
294 |
|
295 |
#######
|
|
296 |
#
|
|
297 |
# Calculation of the phase of the Moon
|
|
298 |
||
299 |
# Age of the Moon, in degrees
|
|
300 |
moon_age = lPP - lambda_sun |
|
301 |
||
302 |
# Phase of the Moon
|
|
303 |
moon_phase = (1 - cos(torad(moon_age))) / 2.0 |
|
304 |
||
305 |
# Calculate distance of Moon from the centre of the Earth
|
|
306 |
moon_dist = (c.moon_smaxis * (1 - c.moon_eccentricity**2))\ |
|
307 |
/ (1 + c.moon_eccentricity * cos(torad(MmP + mEc))) |
|
308 |
||
309 |
# Calculate Moon's angular diameter
|
|
310 |
moon_diam_frac = moon_dist / c.moon_smaxis |
|
311 |
moon_angular_diameter = c.moon_angular_size / moon_diam_frac |
|
312 |
||
313 |
# Calculate Moon's parallax (unused?)
|
|
3.1.1
by Kevin Turner
remove import *, comment out calculations that are not used for phase calculation. |
314 |
# moon_parallax = c.moon_parallax / moon_diam_frac
|
1
by acapnotic
Initial revision |
315 |
|
316 |
res = { |
|
317 |
'phase': fixangle(moon_age) / 360.0, |
|
318 |
'illuminated': moon_phase, |
|
319 |
'age': c.synodic_month * fixangle(moon_age) / 360.0 , |
|
320 |
'distance': moon_dist, |
|
321 |
'angular_diameter': moon_angular_diameter, |
|
322 |
'sun_distance': sun_dist, |
|
323 |
'sun_angular_diameter': sun_angular_diameter |
|
324 |
}
|
|
325 |
||
326 |
return res |
|
327 |
# phase()
|
|
328 |
||
329 |
||
330 |
def phase_hunt(sdate=DateTime.now()): |
|
331 |
"""Find time of phases of the moon which surround the current date.
|
|
332 |
||
333 |
Five phases are found, starting and ending with the new moons
|
|
334 |
which bound the current lunation.
|
|
335 |
"""
|
|
336 |
||
337 |
if not hasattr(sdate,'jdn'): |
|
338 |
sdate = DateTime.DateTimeFromJDN(sdate) |
|
339 |
||
340 |
adate = sdate + DateTime.RelativeDateTime(days=-45) |
|
341 |
||
342 |
k1 = floor((adate.year + ((adate.month - 1) * (1.0/12.0)) - 1900) * 12.3685) |
|
343 |
||
344 |
nt1 = meanphase(adate, k1) |
|
345 |
adate = nt1 |
|
346 |
||
347 |
sdate = sdate.jdn |
|
348 |
||
349 |
while 1: |
|
350 |
adate = adate + c.synodic_month |
|
351 |
k2 = k1 + 1 |
|
352 |
nt2 = meanphase(adate,k2) |
|
353 |
if nt1 <= sdate < nt2: |
|
354 |
break
|
|
355 |
nt1 = nt2 |
|
356 |
k1 = k2 |
|
357 |
||
6
by Kevin Turner
change exception and map usage to be python3 compatible |
358 |
phases = list(map(truephase, |
1
by acapnotic
Initial revision |
359 |
[k1, k1, k1, k1, k2], |
6
by Kevin Turner
change exception and map usage to be python3 compatible |
360 |
[0/4.0, 1/4.0, 2/4.0, 3/4.0, 0/4.0])) |
1
by acapnotic
Initial revision |
361 |
|
362 |
return phases |
|
363 |
# phase_hunt()
|
|
364 |
||
365 |
||
366 |
def meanphase(sdate, k): |
|
367 |
"""Calculates time of the mean new Moon for a given base date.
|
|
368 |
||
369 |
This argument K to this function is the precomputed synodic month
|
|
370 |
index, given by:
|
|
371 |
||
372 |
K = (year - 1900) * 12.3685
|
|
373 |
||
374 |
where year is expressed as a year and fractional year.
|
|
375 |
"""
|
|
376 |
||
377 |
# Time in Julian centuries from 1900 January 0.5
|
|
378 |
if not hasattr(sdate,'jdn'): |
|
379 |
delta_t = sdate - DateTime.DateTime(1900,1,1,12).jdn |
|
380 |
t = delta_t / 36525 |
|
381 |
else: |
|
382 |
delta_t = sdate - DateTime.DateTime(1900,1,1,12) |
|
383 |
t = delta_t.days / 36525 |
|
384 |
||
385 |
# square for frequent use
|
|
386 |
t2 = t * t |
|
387 |
# and cube
|
|
388 |
t3 = t2 * t |
|
389 |
||
390 |
nt1 = ( |
|
391 |
2415020.75933 + c.synodic_month * k + 0.0001178 * t2 - |
|
392 |
0.000000155 * t3 + 0.00033 * dsin(166.56 + 132.87 * t - |
|
393 |
0.009173 * t2) |
|
394 |
)
|
|
395 |
||
396 |
return nt1 |
|
397 |
# meanphase()
|
|
398 |
||
399 |
||
400 |
def truephase(k, tphase): |
|
401 |
"""Given a K value used to determine the mean phase of the new
|
|
402 |
moon, and a phase selector (0.0, 0.25, 0.5, 0.75), obtain the
|
|
403 |
true, corrected phase time."""
|
|
404 |
||
6
by Kevin Turner
change exception and map usage to be python3 compatible |
405 |
apcor = False |
1
by acapnotic
Initial revision |
406 |
|
407 |
# add phase to new moon time
|
|
408 |
k = k + tphase |
|
409 |
# Time in Julian centuries from 1900 January 0.5
|
|
410 |
t = k / 1236.85 |
|
411 |
||
412 |
t2 = t * t |
|
413 |
t3 = t2 * t |
|
414 |
||
415 |
# Mean time of phase
|
|
416 |
pt = ( |
|
417 |
2415020.75933 + c.synodic_month * k + 0.0001178 * t2 - |
|
418 |
0.000000155 * t3 + 0.00033 * dsin(166.56 + 132.87 * t - |
|
419 |
0.009173 * t2) |
|
420 |
)
|
|
421 |
||
422 |
# Sun's mean anomaly
|
|
423 |
m = 359.2242 + 29.10535608 * k - 0.0000333 * t2 - 0.00000347 * t3 |
|
424 |
||
425 |
# Moon's mean anomaly
|
|
426 |
mprime = 306.0253 + 385.81691806 * k + 0.0107306 * t2 + 0.00001236 * t3 |
|
427 |
||
428 |
# Moon's argument of latitude
|
|
429 |
f = 21.2964 + 390.67050646 * k - 0.0016528 * t2 - 0.00000239 * t3 |
|
430 |
||
431 |
if (tphase < 0.01) or (abs(tphase - 0.5) < 0.01): |
|
432 |
||
433 |
# Corrections for New and Full Moon
|
|
434 |
||
435 |
pt = pt + ( |
|
436 |
(0.1734 - 0.000393 * t) * dsin(m) |
|
437 |
+ 0.0021 * dsin(2 * m) |
|
438 |
- 0.4068 * dsin(mprime) |
|
439 |
+ 0.0161 * dsin(2 * mprime) |
|
440 |
- 0.0004 * dsin(3 * mprime) |
|
441 |
+ 0.0104 * dsin(2 * f) |
|
442 |
- 0.0051 * dsin(m + mprime) |
|
443 |
- 0.0074 * dsin(m - mprime) |
|
444 |
+ 0.0004 * dsin(2 * f + m) |
|
445 |
- 0.0004 * dsin(2 * f - m) |
|
446 |
- 0.0006 * dsin(2 * f + mprime) |
|
447 |
+ 0.0010 * dsin(2 * f - mprime) |
|
448 |
+ 0.0005 * dsin(m + 2 * mprime) |
|
449 |
)
|
|
450 |
||
6
by Kevin Turner
change exception and map usage to be python3 compatible |
451 |
apcor = True |
1
by acapnotic
Initial revision |
452 |
elif (abs(tphase - 0.25) < 0.01) or (abs(tphase - 0.75) < 0.01): |
453 |
||
454 |
pt = pt + ( |
|
455 |
(0.1721 - 0.0004 * t) * dsin(m) |
|
456 |
+ 0.0021 * dsin(2 * m) |
|
457 |
- 0.6280 * dsin(mprime) |
|
458 |
+ 0.0089 * dsin(2 * mprime) |
|
459 |
- 0.0004 * dsin(3 * mprime) |
|
460 |
+ 0.0079 * dsin(2 * f) |
|
461 |
- 0.0119 * dsin(m + mprime) |
|
462 |
- 0.0047 * dsin(m - mprime) |
|
463 |
+ 0.0003 * dsin(2 * f + m) |
|
464 |
- 0.0004 * dsin(2 * f - m) |
|
465 |
- 0.0006 * dsin(2 * f + mprime) |
|
466 |
+ 0.0021 * dsin(2 * f - mprime) |
|
467 |
+ 0.0003 * dsin(m + 2 * mprime) |
|
468 |
+ 0.0004 * dsin(m - 2 * mprime) |
|
469 |
- 0.0003 * dsin(2 * m + mprime) |
|
470 |
)
|
|
471 |
if (tphase < 0.5): |
|
472 |
# First quarter correction
|
|
473 |
pt = pt + 0.0028 - 0.0004 * dcos(m) + 0.0003 * dcos(mprime) |
|
474 |
else: |
|
475 |
# Last quarter correction
|
|
476 |
pt = pt + -0.0028 + 0.0004 * dcos(m) - 0.0003 * dcos(mprime) |
|
6
by Kevin Turner
change exception and map usage to be python3 compatible |
477 |
apcor = True |
1
by acapnotic
Initial revision |
478 |
|
479 |
if not apcor: |
|
6
by Kevin Turner
change exception and map usage to be python3 compatible |
480 |
raise ValueError( |
481 |
"TRUEPHASE called with invalid phase selector", |
|
482 |
tphase) |
|
1
by acapnotic
Initial revision |
483 |
|
484 |
return DateTime.DateTimeFromJDN(pt) |
|
485 |
# truephase()
|
|
486 |
||
487 |
||
488 |
def kepler(m, ecc): |
|
489 |
"""Solve the equation of Kepler."""
|
|
490 |
||
491 |
epsilon = 1e-6 |
|
492 |
||
493 |
m = torad(m) |
|
494 |
e = m |
|
495 |
while 1: |
|
496 |
delta = e - ecc * sin(e) - m |
|
497 |
e = e - delta / (1.0 - ecc * cos(e)) |
|
498 |
||
499 |
if abs(delta) <= epsilon: |
|
500 |
break
|
|
501 |
||
502 |
return e |
|
503 |
||
504 |
#
|
|
505 |
##
|
|
506 |
#
|
|
507 |
||
508 |
if __name__ == '__main__': |
|
509 |
m = MoonPhase() |
|
510 |
s = """The moon is %s, %.1f%% illuminated, %.1f days old.""" %\ |
|
511 |
(m.phase_text, m.illuminated * 100, m.age) |
|
6
by Kevin Turner
change exception and map usage to be python3 compatible |
512 |
print (s) |