1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
|
include <math.h>
# AST_VBARY -- Radial velocity component of center of the Earth relative to
# to the barycenter of the Earth-Moon system.
procedure ast_vbary (ra, dec, epoch, v)
double ra # Right ascension of observation (hours)
double dec # Declination of observation (degrees)
double epoch # Julian epoch of observation
double v # Component of orbital velocity (km/s)
double t, oblq, omega, llong, lperi, inclin, em, anom, vmoon
double r, d, l, b, lm, bm, ast_julday()
begin
# T is the number of Julian centuries since J1900.
t = (ast_julday (epoch) - 2415020) / 36525.
# OBLQ is the mean obliquity of the ecliptic
# OMEGA is the longitude of the mean ascending node
# LLONG is the mean lunar longitude (should be 13.1763965268)
# LPERI is the mean lunar longitude of perigee
# INCLIN is the inclination of the lunar orbit to the ecliptic
# EM is the eccentricity of the lunar orbit (dimensionless)
# All quantities except the eccentricity are in degrees.
oblq = 23.452294d0 -
t * (0.0130125d0 + t * (0.00000164d0 - t * 0.000000503d0))
omega = 259.183275d0 -
t * (1934.142008d0 + t * (0.002078d0 + t * 0.000002d0))
llong = 270.434164d0 +
t * (481267.88315d0 + t * (-0.001133d0 + t * 0.0000019d0)) - omega
lperi = 334.329556d0 +
t * (4069.034029d0 - t * (0.010325d0 + t * 0.000012d0)) - omega
em = 0.054900489d0
inclin = 5.1453964d0
# Determine true longitude. Compute mean anomaly, convert to true
# anomaly (approximate formula), and convert back to longitude.
# The mean anomaly is only approximate because LPERI should
# be the true rather than the mean longitude of lunar perigee.
lperi = DEGTORAD (lperi)
llong = DEGTORAD (llong)
anom = llong - lperi
anom = anom + (2. * em - 0.25 * em**3) * sin (anom) + 1.25 * em**2 *
sin (2. * anom) + 13./12. * em**3 * sin (3. * anom)
llong = anom + lperi
# L and B are the ecliptic longitude and latitude of the observation.
# LM and BM are the lunar longitude and latitude of the observation
# in the lunar orbital plane relative to the ascending node.
r = DEGTORAD (ra * 15)
d = DEGTORAD (dec)
omega = DEGTORAD (omega)
oblq = DEGTORAD (oblq)
inclin = DEGTORAD (inclin)
call ast_coord (double (0.), double (0.), double (-HALFPI),
HALFPI - oblq, r, d, l, b)
call ast_coord (omega, double (0.), omega-HALFPI, HALFPI-inclin,
l, b, lm, bm)
# VMOON is the component of the lunar velocity perpendicular to the
# radius vector. V is the projection onto the line of sight to the
# observation of the velocity of the Earth's center with respect to
# the Earth-Moon barycenter. The 81.53 is the ratio of the Earth's
# mass to the Moon's mass.
vmoon = (TWOPI / 27.321661d0) *
384403.12040d0 / sqrt (1. - em**2) / 86400.
v = vmoon * cos (bm) * (sin (llong - lm) - em * sin (lperi - lm))
v = v / 81.53
end
|
