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
77
78
79
80
81
82
|
include <math.h>
# AST_HJD -- Helocentric Julian Day from Epoch
procedure ast_hjd (ra, dec, epoch, lt, hjd)
double ra # Right ascension of observation (hours)
double dec # Declination of observation (degrees)
double epoch # Julian epoch of observation
double lt # Light travel time in seconds
double hjd # Helocentric Julian Day
double ast_julday()
begin
call ast_jd_to_hjd (ra, dec, ast_julday(epoch), lt, hjd)
end
# AST_JD_TO_HJD -- Helocentric Julian Day from UT Julian date
procedure ast_jd_to_hjd (ra, dec, jd, lt, hjd)
double ra # Right ascension of observation (hours)
double dec # Declination of observation (degrees)
double jd # Geocentric Julian date of observation
double lt # Light travel time in seconds
double hjd # Helocentric Julian Day
double t, manom, lperi, oblq, eccen, tanom, slong, r, d, l, b, rsun
begin
# JD is the geocentric Julian date.
# T is the number of Julian centuries since J1900.
t = (jd - 2415020d0) / 36525d0
# 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)
# R is the distance to the Sun.
rsun = (1. - eccen**2) / (1. + eccen * cos (tanom))
# LT is the light travel difference to the Sun.
lt = -0.005770d0 * rsun * cos (b) * cos (l - slong)
hjd = jd + lt
end
|
