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
|
include <math.h>
# AST_VORBIT -- Radial velocity component of the Earth-Moon barycenter
# relative to the Sun.
procedure ast_vorbit (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, manom, lperi, oblq, eccen, tanom, slong, r, d, l, b, vorb
double ast_julday()
begin
# T is the number of Julian centuries since J1900.
t = (ast_julday (epoch) - 2415020d0) / 36525.
# MANOM is the mean anomaly of the Earth's orbit (degrees)
# LPERI is the mean longitude of perihelion (degrees)
# OBLQ is the mean obliquity of the ecliptic (degrees)
# ECCEN is the eccentricity of the Earth's orbit (dimensionless)
manom = 358.47583d0 +
t * (35999.04975d0 - t * (0.000150d0 + t * 0.000003d0))
lperi = 101.22083d0 +
t * (1.7191733d0 + t * (0.000453d0 + t * 0.000003d0))
oblq = 23.452294d0 -
t * (0.0130125d0 + t * (0.00000164d0 - t * 0.000000503d0))
eccen = 0.01675104d0 - t * (0.00004180d0 + t * 0.000000126d0)
# Convert to principle angles
manom = mod (manom, 360.0D0)
lperi = mod (lperi, 360.0D0)
# Convert to radians
r = DEGTORAD (ra * 15)
d = DEGTORAD (dec)
manom = DEGTORAD (manom)
lperi = DEGTORAD (lperi)
oblq = DEGTORAD (oblq)
# TANOM is the true anomaly (approximate formula) (radians)
tanom = manom + (2 * eccen - 0.25 * eccen**3) * sin (manom) +
1.25 * eccen**2 * sin (2 * manom) +
13./12. * eccen**3 * sin (3 * manom)
# SLONG is the true longitude of the Sun seen from the Earth (radians)
slong = lperi + tanom + PI
# L and B are the longitude and latitude of the star in the orbital
# plane of the Earth (radians)
call ast_coord (double (0.), double (0.), double (-HALFPI),
HALFPI - oblq, r, d, l, b)
# VORB is the component of the Earth's orbital velocity perpendicular
# to the radius vector (km/s) where the Earth's semi-major axis is
# 149598500 km and the year is 365.2564 days.
vorb = ((TWOPI / 365.2564d0) *
149598500.d0 / sqrt (1. - eccen**2)) / 86400.d0
# V is the projection onto the line of sight to the observation of
# the velocity of the Earth-Moon barycenter with respect to the
# Sun (km/s).
v = vorb * cos (b) * (sin (slong - l) - eccen * sin (lperi - l))
end
|
