How do I determine astronomical transit times?

Question

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?

Answer 1

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
 ]