# How do I determine astronomical transit times?

Mathematica, through the "legacy" Scientific Astronomer package, used to have the ability to easily determine accurate transit times of astronomical objects simply, with Culmination[object, neardate]; but now lacks this feature. Is there a way to use AstronomicalData to do this?

• Just in case, would you still be cool with something that doesn't use AstronomicalData[]? Oct 22, 2012 at 23:00
• Yes, though I'm hoping to avoid (a) writing too much code and (b) concerning myself too much with reconciling different sources and benchmarks. I worry that if I try to kludge together too many methods and datasources I may end up with inconsistencies. Oct 22, 2012 at 23:09
• Looking at the algorithm description by Meeus, it would seem that the formula for transit times involves $\Delta T$, which, if memory serves, has to be obtained from an external source like this. Oct 23, 2012 at 0:34
• @J.M.: I don't own Meeus. It keeps coming up though. Oct 23, 2012 at 3:11
• Oct 23, 2012 at 3:31

This is not a full answer, but more a response to J.M.'s comment and provides a routine to calculate $\Delta T$ which was sitting on my hard disk. This is intended as a starting point for further calculations.

deltaT::usage =
"deltaT[date] calculate the arithmetic difference, in seconds, \
between the Terrestrial Dynamical Time (TD) and the  Universal \
Time(UT).\n ΔT = TD - UT is given in seconds.";
deltaT::yearx = "Invalid year.";

Options[deltaT] = {Method -> Interpolation};
deltaT[{y0_, m0_: 1, d0_: 1, h0_: 0, min0_: 0, s0_: 0},
opts : OptionsPattern[]] :=
Module[
{y = y0 + (m0 - 0.5)/12, u, ΔT, t, tab, meth},
(*from:http://eclipse.gsfc.nasa.gov/SEcat5/deltat.html*)
tab = {{-500, 17190}, {-400, 15530}, {-300, 14080}, {-200,
12790}, {-100, 11640}, {0, 10580}, {100, 9600}, {200,
8640}, {300, 7680}, {400, 6700}, {500, 5710}, {600, 4740}, {700,
3810}, {800, 2960}, {900, 2200}, {1000, 1570}, {1100,
1090}, {1200, 740}, {1300, 490}, {1400, 320}, {1500, 200}, {1600,
120}, {1700, 9}, {1750, 13}, {1800, 14}, {1850,
7}, {1900, -3}, {1950, 29}, {1955., 31.1}, {1960.,
33.2}, {1965., 35.7}, {1970., 40.2}, {1975.,
45.5}, {1980., 50.5}, {1985., 54.3}, {1990.,
56.9}, {1995., 60.8}, {2000., 63.8}, {2005., 64.7}};
meth = OptionValue[Method];
If[! MemberQ[{Interpolation, Polynomials}, meth],
meth = Interpolation];
If[meth === Interpolation,
If[-500 <= y0 <= 2005, ΔT =
Interpolation[tab, y],(*extrapolation beyond observed values*)
t = (y - 1820)/100;
ΔT = -20 + 32*t^2],
Which[y0 < -500, u = (y - 1820)/100;
ΔT = -20 + 32*u^2, -500 <= y0 <=
500,(*error not larger than 4 seconds*)u = y/100;
ΔT =
10583.6 - 1014.41*u + 33.78311*u^2 - 5.952053*u^3 -
0.1798452*u^4 + 0.022174192*u^5 + 0.0090316521*u^6,
500 < y0 <= 1600, u = (y - 1000)/100;
ΔT =
1574.2 - 556.01*u + 71.23472*u^2 + 0.319781*u^3 -
0.8503463*u^4 - 0.005050998*u^5 + 0.0083572073*u^6,
1600 < y0 <= 1700, t = y - 1600;
ΔT = 120 - 0.9808*t - 0.01532*t^2 + t^3/7129,
1700 < y0 <= 1800, t = y - 1700;
ΔT =
8.83 + 0.1603*t - 0.0059285*t^2 + 0.00013336*t^3 - t^4/1174000,
1800 < y0 <= 1860,
t = y - 1800; ΔT =
13.72 - 0.332447*t + 0.0068612*t^2 + 0.0041116*t^3 -
0.00037436*t^4 + 0.0000121272*t^5 - 0.0000001699*t^6 +
0.000000000875*t^7, 1860 < y0 <= 1900, t = y - 1860;
ΔT =
7.62 + 0.5737*t - 0.251754*t^2 + 0.01680668*t^3 -
0.0004473624*t^4 + t^5/233174, 1900 < y0 <= 1920, t = y - 1900;
ΔT = -2.79 + 1.494119*t - 0.0598939*t^2 +
0.0061966*t^3 - 0.000197*t^4, 1920 < y0 <= 1941, t = y - 1920;
ΔT =
21.20 + 0.84493*t - 0.076100*t^2 + 0.0020936*t^3,
1941 < y0 <= 1961, t = y - 1950;
ΔT = 29.07 + 0.407*t - t^2/233 + t^3/2547,
1961 < y0 <= 1986, t = y - 2000;
ΔT =
63.86 + 0.3345*t - 0.060374*t^2 + 0.0017275*t^3 +
0.000651814*t^4 + 0.00002373599*t^5, 1986 < y0 <= 2005,
t = y - 2000;
ΔT =
63.86 + 0.3345*t - 0.060374*t^2 + 0.0017275*t^3 +
0.000651814*t^4 + 0.00002373599*t^5, 2005 < y0 <= 2050,
t = y - 2000;
ΔT =
62.92 + 0.32217*t +
0.005589*
t^2,
(*This expression is derived from estimated values of ΔT in the years 2010 and 2050. The value for 2010 (66.9 seconds) is based on a linearly extrapolation from 2005 using 0.39 seconds/year (average from 1995 to 2005).The value for 2050 (93 seconds) is linearly extrapolated from 2010 using 0.66 seconds year (average rate from 1901 to 2000).*)
2050 < y0 <=
2150, ΔT = -20 + 32*((y - 1820)/100)^2 -
0.5628*(2150 -
y),
(*The last term is introduced to eliminate the discontinuity at 2050.*)
2150 < y0, u = (y - 1820)/100;
ΔT = -20 + 32*u^2, True, Message[deltaT::yearx]];
(*All values of ΔT based on Morrison and Stephenson[2004] assume a value for the Moon's secular acceleration of-26 arcsec/cy^2. However the ELP-2000  82 lunar ephemeris employed in the Canon uses a slightly different value of-25.858 arcsec/cy^2. Thus a small correction "c" must be added to the values derived from the polynomial expressions for ΔT before they can be used in the Canon c=-0.000012932*(y-1955)^2. Since the  values of ΔT for the interval 1955 to 2005 were derived independent of any lunar ephemeris, no correction is needed for this period.*)
(*References Morrison, L.and Stephenson,F.R., "Historical Values of the Earth's Clock Error ΔT and the Calculation of Eclipses",J.Hist.Astron.,Vol .35 Part 3, August 2004,No .120,pp 327-336 (2004).Stephenson F.R and Houlden M.A., Atlas of Historical Eclipse Maps,Cambridge Univ.Press.,Cambridge, 1986.*)
];
ΔT
]

• +1, and holy hell, that's complicated... :o Oct 23, 2012 at 8:44
• About 1999 I bought Scientific Astronomer add-on, with its GetLocation function which could calculate the location of a satellite from its Two-Line Elements (TLE) at specified times. This calculation was done in local memory based on specified TLEs. The current implicit lookup of a satellite's location using forms such as: "EntityValue[Entity["Satellite", "48938"], Dated["Position", DateObject[{2022, 7, 4, 12, 0, 2}]] ]" or ["PositionLine"], appears to get the info from Wolfram Servers, much slower than I would like. Is there an implementation of the SGP4 propagator algorithm available? Jul 22, 2022 at 3:41