# How do I compute sidereal time and JD?

There was once a package that provided a wide range of functions for computing astronomical information. This package now appears to be out of date; but as near as I can tell, many of the useful functions it provided have not been replicated in AstronomicalData. Are there, for example, functions I can use to compute sidereal time, or Julian Date (JD)?

• A good online source for astronomical calculations is Paul Schlyter's page. Code is in C, of course, but the descriptions are nicely written. Sep 28, 2012 at 19:24
• @cormullion +1 for the reference. Great page that I also used for many of my implementations. Sep 28, 2012 at 22:42
• N.B. SiderealTime[] and JulianDate[] are now built-in, but for some reason, I have been unable to match the outputs of those functions against Meeus's test cases. Aug 2, 2016 at 8:26
• @J.M.: Yeah, there's something — um — unique about the way MMA treats dates in general. Aug 2, 2016 at 11:35
• In the meantime, did you already get yourself a copy of Meeus like I told you before? :) Aug 2, 2016 at 11:38

Jean Meeus's Astronomical Algorithms (as well as the related book Astronomical Formulæ for Calculators) is what you should start looking at whenever you need to deal with algorithms for quantities of astronomical interest.

For instance, here is a translation of Meeus's method for the Julian Date:

Options[jd] = {"Calendar" -> "Gregorian"};

jd[{yr_Integer, mo_Integer, da_?NumericQ, rest__}, opts : OptionsPattern[]] :=
Module[{y = yr, m = mo, h}, If[m < 3, y--; m += 12];
h = Switch[OptionValue["Calendar"],
"Gregorian", (Quotient[#, 4] - # + 2) &[Quotient[y, 100]],
"Julian", 0,
_, Return[$Failed]]; Floor[365.25 y] + Floor[30.6001 (m + 1)] + da + FromDMS[{rest}]/24 + 1720994.5 + h ] jd[{yr_Integer, mo_Integer, da_?NumericQ}, opts : OptionsPattern[]] := jd[{yr, mo, da, 0, 0, 0}, opts] jd[opts : OptionsPattern[]] := jd[DateList[], opts]  Some examples: jd[{1957, 10, 4}] // InputForm 2.4361155*^6 jd[{1977, 4, 26}] // InputForm 2.4432595*^6  If you're sure you'll only be dealing with dates in the Gregorian calendar, one could use AbsoluteTime[] to compute JD, like so: jd[args___] := AbsoluteTime[args]/86400 + 2.4150205*^6  Here's a short routine for Greenwich mean sidereal time (in {hours, minutes, seconds} format), based on formulæ given here: GMST[{yr_Integer, mo_Integer, da_?NumericQ, rest__}, opts : OptionsPattern[]] := DMSList[Mod[6.697374558 + 0.06570982441908 (jd[{yr, mo, da}, opts] - 2.451545*^6) + 1.00273790935 FromDMS[{rest}] + 0.000026 ((jd[{yr, mo, da, rest}, opts] - 2.451545*^6)/36525)^2, 24]] GMST[{yr_Integer, mo_Integer, da_?NumericQ}, opts : OptionsPattern[]] := GMST[{yr, mo, da, 0, 0, 0}, opts] GMST[opts : OptionsPattern[]] := GMST[DateList[], opts]  Some examples: GMST[{1987, 4, 10, 19, 21, 0}] {8, 34, 57.089598878778816} GMST[{2001, 10, 3, 6, 30, 0}] {7, 18, 8.328661911859285}  One could use the output of FindGeoLocation[] to adjust the value of GMST[] to the corresponding local mean sidereal time, of course. I'll leave that as an exercise (the method needed is given in the USNO link I gave). Getting the apparent sidereal time from the mean sidereal time requires a rather complicated correction term involving a series of trigonometric functions (the so-called equation of equinoxes). Here's a bunch of routines for computing Greenwich apparent sidereal time: (* IAU 2000B coefficients (77 terms); from http://dx.doi.org/10.1023/A:1021762727016 *)$nuϕcoeffs = {{0, 0, 0, 0, 1}, {0, 0, 2, -2, 2}, {0, 0, 2, 0, 2}, {0, 0, 0, 0, 2},
{0, 1, 0, 0, 0}, {0, 1, 2, -2, 2}, {1, 0, 0, 0, 0}, {0, 0, 2, 0, 1},
{1, 0, 2, 0, 2}, {0, -1, 2, -2, 2}, {0, 0, 2, -2, 1}, {-1, 0, 2, 0, 2},
{-1, 0, 0, 2, 0}, {1, 0, 0, 0, 1}, {-1, 0, 0, 0, 1}, {-1, 0, 2, 2, 2},
{1, 0, 2, 0, 1}, {-2, 0, 2, 0, 1}, {0, 0, 0, 2, 0}, {0, 0, 2, 2, 2},
{0, -2, 2, -2, 2}, {-2, 0, 0, 2, 0}, {2, 0, 2, 0, 2}, {1, 0, 2, -2, 2},
{-1, 0, 2, 0, 1}, {2, 0, 0, 0, 0}, {0, 0, 2, 0, 0}, {0, 1, 0, 0, 1},
{-1, 0, 0, 2, 1}, {0, 2, 2, -2, 2}, {0, 0, -2, 2, 0}, {1, 0, 0, -2, 1},
{0, -1, 0, 0, 1}, {-1, 0, 2, 2, 1}, {0, 2, 0, 0, 0}, {1, 0, 2, 2, 2},
{-2, 0, 2, 0, 0}, {0, 1, 2, 0, 2}, {0, 0, 2, 2, 1}, {0, -1, 2, 0, 2},
{0, 0, 0, 2, 1}, {1, 0, 2, -2, 1}, {2, 0, 2, -2, 2}, {-2, 0, 0, 2, 1},
{2, 0, 2, 0, 1}, {0, -1, 2, -2, 1}, {0, 0, 0, -2, 1}, {-1, -1, 0, 2, 0},
{2, 0, 0, -2, 1}, {1, 0, 0, 2, 0}, {0, 1, 2, -2, 1}, {1, -1, 0, 0, 0},
{-2, 0, 2, 0, 2}, {3, 0, 2, 0, 2}, {0, -1, 0, 2, 0}, {1, -1, 2, 0, 2},
{0, 0, 0, 1, 0}, {-1, -1, 2, 2, 2}, {-1, 0, 2, 0, 0}, {0, -1, 2, 2, 2},
{-2, 0, 0, 0, 1}, {1, 1, 2, 0, 2}, {2, 0, 0, 0, 1}, {-1, 1, 0, 1, 0},
{1, 1, 0, 0, 0}, {1, 0, 2, 0, 0}, {-1, 0, 2, -2, 1}, {1, 0, 0, 0, 2},
{-1, 0, 0, 1, 0}, {0, 0, 2, 1, 2}, {-1, 0, 2, 4, 2}, {-1, 1, 0, 1, 1},
{0, -2, 2, -2, 1}, {1, 0, 2, 2, 1}, {-2, 0, 2, 2, 2}, {-1, 0, 0, 0, 2},
{1, 1, 2, -2, 2}};

$Δψcoeffs = {{-172064161, -174666, 33386}, {-13170906, -1675, 13696}, {-2276413, -234, 2796}, {2074554, 207, -698}, {1475877, -3633, 11817}, {-516821, 1226, -524}, {711159, 73, -872}, {-387298, -367, 380}, {-301461, -36, 816}, {215829, -494, 111}, {128227, 137, 181}, {123457, 11, 19}, {156994, 10, -168}, {63110, 63, 27}, {-57976, -63, -189}, {-59641, -11, 149}, {-51613, -42, 129}, {45893, 50, 31}, {63384, 11, -150}, {-38571, -1, 158}, {32481, 0, 0}, {-47722, 0, -18}, {-31046, -1, 131}, {28593, 0, -1}, {20441, 21, 10}, {29243, 0, -74}, {25887, 0, -66}, {-14053, -25, 79}, {15164, 10, 11}, {-15794, 72, -16}, {21783, 0, 13}, {-12873, -10, -37}, {-12654, 11, 63}, {-10204, 0, 25}, {16707, -85, -10}, {-7691, 0, 44}, {-11024, 0, -14}, {7566, -21, -11}, {-6637, -11, 25}, {-7141, 21, 8}, {-6302, -11, 2}, {5800, 10, 2}, {6443, 0, -7}, {-5774, -11, -15}, {-5350, 0, 21}, {-4752, -11, -3}, {-4940, -11, -21}, {7350, 0, -8}, {4065, 0, 6}, {6579, 0, -24}, {3579, 0, 5}, {4725, 0, -6}, {-3075, 0, -2}, {-2904, 0, 15}, {4348, 0, -10}, {-2878, 0, 8}, {-4230, 0, 5}, {-2819, 0, 7}, {-4056, 0, 5}, {-2647, 0, 11}, {-2294, 0, -10}, {2481, 0, -7}, {2179, 0, -2}, {3276, 0, 1}, {-3389, 0, 5}, {3339, 0, -13}, {-1987, 0, -6}, {-1981, 0, 0}, {4026, 0, -353}, {1660, 0, -5}, {-1521, 0, 9}, {1314, 0, 0}, {-1283, 0, 0}, {-1331, 0, 8}, {1383, 0, -2}, {1405, 0, 4}, {1290, 0, 0}}; (* Coefficients of "complementary terms" for equation of equinoxes (12 terms); from IERS 2010 conventions *)$eeAcoeffs =
{2640.96, 63.52, 11.75, 11.21, -4.55, 2.02, 1.98, -1.72, -1.41, -1.26, -.63, -.63};

$eeϕcoeffs = {{0, 0, 0, 0, 1}, {0, 0, 0, 0, 2}, {0, 0, 2, -2, 3}, {0, 0, 2, -2, 1}, {0, 0, 2, -2, 2}, {0, 0, 2, 0, 3}, {0, 0, 2, 0, 1}, {0, 0, 0, 0, 3}, {0, 1, 0, 0, 1}, {0, 1, 0, 0, -1}, {1, 0, 0, 0, 1}, {1, 0, 0, 0, -1}}; (* mean obliquity of the ecliptic ε; formula 5.12 from USNO Circular 179 *) MeanEclipticObliquity[args___] := Module[{T}, T = (jd[args] - 2451545)/36525; (84381.406 + T (-46.836769 + T (-0.0001831 + T (0.0020034 + T (-5.76*^-7 - 4.34*^-8 T)))))/3600 ] EquationOfEquinoxes[args___] := Module[{T, ℓ☽, ℓ☉, ℱ, \[ScriptCapitalD], Ω, ϕterms}, T = (jd[args] - 2451545)/36525; (* ℓ☽ - mean lunar anomaly; ℓ☉ - mean solar anomaly; ℱ - mean lunar argument of latitude; \[ScriptCapitalD] - mean lunar elongation; Ω - mean longitude of lunar ascending node; *) {ℓ☽, ℓ☉, ℱ, \[ScriptCapitalD], Ω} = {485868.249036 + T (1.7179159232178*^9 + T (31.8792 + T (0.051635 - 0.0002447 T))), 1.28710479305*^6 + T (1.295965810481*^8 + T (-0.5532 + T (0.000136 - 0.00001149 T))), 335779.526232 + T (1.7395272628478*^9 + T (-12.7512 + T (-0.001037 + 4.17*^-6 T))), 1.07226070369*^6 + T (1.602961601209*^9 + T (-6.3706 + T (0.006593 - 0.00003169 T))), 450160.398036 + T (-6.9628905431*^6 + T (7.4722 + T (0.007702 - 0.00005939 T)))}; ϕterms =$nuϕcoeffs.{ℓ☽, ℓ☉, ℱ, \[ScriptCapitalD], Ω} Degree/3600;

(Total[((#1 + #2 T) Sin[#4] + #3 Cos[#4]) & @@@ MapThread[Append, {$Δψcoeffs, ϕterms}], Method -> "CompensatedSummation"]/3.6*^10) Cos[MeanEclipticObliquity[args] Degree] +$eeAcoeffs.Sin[$eeϕcoeffs.{ℓ☽, ℓ☉, ℱ, \[ScriptCapitalD], Ω} Degree/3600]/3.6*^9 ] GAST[{yr_Integer, mo_Integer, da_?NumericQ, rest__}, opts : OptionsPattern[]] := DMSList[Mod[6.697374558 + 0.06570982441908 (jd[{yr, mo, da}, opts] - 2.451545*^6) + 1.00273790935 FromDMS[{rest}] + 0.000026 ((jd[{yr, mo, da, rest}, opts] - 2.451545*^6)/36525)^2 + EquationOfEquinoxes[{yr, mo, da, rest}]/15, 24]] GAST[{yr_Integer, mo_Integer, da_?NumericQ}, opts : OptionsPattern[]] := GAST[{yr, mo, da, 0, 0, 0}, opts] GAST[opts : OptionsPattern[]] := GAST[DateList[], opts]  An example: GAST[{2012, 10, 1, 0, 0, 0}] // InputForm {0, 40, 31.48976449269867}  • Great! A pity that I can't give you another +1 for each edit :) Sep 29, 2012 at 17:34 • N.B. Formatting the coefficients for the involved series took a lot more time than writing out the programs for evaluating the series; I sure hope the effort expended was worth it... if one wants to use the more accurate IAU 2000A series instead (1365 terms!), writing the corresponding routines for that is left as an exercise. Oct 1, 2012 at 1:51 • BTW: Still mining this for useful stuff. +2, if I could! Nov 20, 2012 at 16:16 • @rax, I'm glad to hear the effort I put into this was worthwhile. You're welcome. :) Nov 20, 2012 at 16:22 • @J.M.: You may want to add If[y < 0, y++];(* Correct negative years to astronomical reckoning *) to the top of your jd function. Mathematica (9.0) does not use astronomical year reckoning, though 0 is still a valid synonym for the year before 1 CE (try DateList[{0}]). Dec 10, 2012 at 18:46 Search for brfASTRO.m which is a fantastic astronomy package that Peter Breitfeld wrote. He offers wonderful material on his homepage. I used many of his routines (thank you Peter!) Here is my (long but untested) version which I probably also copied partly and forgot from where. Please contact me if you feel that proper credit is due! This is a small part of a large package I compiled and collected over time. I can give more infos on what else is in there if needed. Please go ahead and test! Needs["Calendar"] JulianDate::usage = "JulianDate[year, month, day}] calculates the Julian date (JD) for a given date."; FromJulianDay::usage = "FromJulianDay[JulianDay] converts a Julian day to a date in Gregorian calendar."; ModifiedJulianDate::usage = "The ModifiedJulianDate corresponds to Julian Date - 2400000.5, i.e. MJD=0.0 corresponds to 1858 November 17 at \!$$\*SuperscriptBox[\\"0\", \"h\"]$$ UT. Contrary to the JD, the Modified Julian Day begins at Greenwich mean midnight."; LocalSiderealTime::timetx = "Unknown Sidereal time type."; LocalSiderealTime::usage = "LocalSiderealTime[date, East Longitude] calculates the Local Sidereal Time (mean or apparent) for a given date (or Julian day JD) and a given east longitude [in decimal degree]. The output is given in decimal degree. A high precision algorithm is used."; JulianDate[{y0_, m0_, d0_}] := DaysBetween[{-4712, 1, 1, 12, 0, 0}, CalendarChange[{y0, m0, d0}, Gregorian, Julian], Calendar -> Julian] + (0/24. + 0/(24 60) + 0/(24 60 60) - 0.5) JulianDate[{y0_, m0_, d0_, h0_, min0_, s0_}] := DaysBetween[{-4712, 1, 1, 12, 0, 0}, CalendarChange[{y0, m0, d0, h0, min0, s0}, Gregorian, Julian], Calendar -> Julian] + (h0/24. + min0/(24 60) + s0/(24 60 60) - 0.5) ModifiedJulianDate[{y0_, m0_, d0_, h0_, min0_, s0_}] := JulianDate[{y0, m0, d0, h0, min0, s0}] - 2400000.5 ModifiedJulianDate[{y0_, m0_, d0_}] := JulianDate[{y0, m0, d0, 12, 0, 0}] - 2400000.5 SetAttributes[LocalSiderealTime, Listable]; Options[LocalSiderealTime] = {SiderealTimeType -> "Mean"}; LocalSiderealTime[{y_, m_, d_, h_: 0, min_: 0, sec_: 0}, LongitudeEast_, opts : OptionsPattern[]] := LocalSiderealTime[JulianDate[{y, m, d, h, min, sec}], LongitudeEast, opts] LocalSiderealTime[JD_?NumberQ, LongitudeEast_, opts : OptionsPattern[]] := Module[{TJD, T, DayFraction, GMST0UT, nut, Obl, EquationOfEquinox, type, lst}, type = OptionValue[SiderealTimeType]; (* % convert JD to integer day+fraction of day *) TJD = Floor[JD - 0.5] + 0.5; DayFraction = JD - TJD; T = (TJD - 2451545.0)/36525.0; GMST0UT = 24110.54841 + 8640184.812866 T + 0.093104 T^2 - 6.2 10^-6 T^3; (* % convert to fraction of day in range[0 1) *) GMST0UT = GMST0UT/86400.0; GMST0UT = GMST0UT - Floor[GMST0UT]; (* lst = GMST0UT+1.0027379093*DayFraction+LongitudeEast/(2. Pi); Longitude should be given in Degree *) lst = GMST0UT + 1.0027379093*DayFraction + LongitudeEast/(360); lst = lst - Floor[lst]; Which[ type == "Mean", lst, type == "Apparent", nut = Nutation[JD]; Obl = Obliquity[JD]; EquationOfEquinox = (3600/Degree) nut[[1]] Cos[Obl]/15; lst = lst + EquationOfEquinox/86400, True, Message[LST::timetx] ]; lst = lst*360] SetAttributes[Nutation, Listable]; Options[Nutation] = {FullPrecision -> True}; Nutation[{y_, m_, d_, h_: 0, min_: 0, sec_: 0}, opts : OptionsPattern[]] := Nutation[JulianDate[{y, m, d, h, min, sec}], opts] Nutation[JD_?NumberQ, opts : OptionsPattern[]] := Module[{T, d, M, Mt, F, Om, DLon, DObl, Nut, fullPrecision, rad = 1/Degree}, fullPrecision = OptionValue[FullPrecision]; T = (JD - 2451545.0)/36525.0; (* Mean elongation of the Moon from the Sun *) d = (297.85036 + 445267.111480*T - 0.0019142 T^2 + T^3/189474)/rad; (* Mean anomaly of the Sun (Earth) *) M = (357.52772 + 35999.050340*T - 0.0001603 T^2 - T^3/300000)/rad; (* Mean anomaly of the Moon *) Mt = (134.96298 + 477198.867398*T + 0.0086972 T^2 + T^3/56250)/rad; (* Moon's argument of latitude *) F = (93.27191 + 483202.017538*T - 0.0036825 T^2 + T^3/327270)/rad; (* Longitude of ascending node of the Moon's mean orbit *) Om = (125.04452 - 1934.136261*T + 0.0020708 T^2 + T^3/450000)/rad; (* nutation in longitude[0."0001] *) DLon = (-171996 - 174.2 T) Sin[Om] + (-13187 + 1.6 T) Sin[ 2. (-d + F + Om)] + (-2274 - 0.2 T) Sin[ 2. (F + Om)] + (2062 + 0.2 T) Sin[2. Om] + (1426 - 3.4 T) Sin[ M] + (712 + 0.1 T) Sin[Mt] + (-517 + 1.2 T) Sin[ 2. (-d + F + Om) + M] + (-386 - 0.4 T) Sin[ 2. F + Om] + (-301) Sin[ Mt + 2. (F + Om)] + (217 - 0.5 T) Sin[-M + 2. (-d + F + Om)] + (-158) Sin[-2. d + Mt] + (129 + 0.1 T) Sin[ 2. (-d + F) + Om] + (123) Sin[-Mt + 2. (F + Om)] + (63) Sin[ 2. d] + (63 + 0.1 T) Sin[ Mt + Om] + (-59) Sin[-Mt + 2. (d + F + Om)] + (-58 - 0.1 T) Sin[-Mt + Om] + (-51) Sin[Mt + 2. F + Om] + (48) Sin[ 2. (-d + Mt)] + (46) Sin[Om + 2. (-Mt + F)] + (-38) Sin[ 2. (d + F + Om)] + (-31) Sin[2. (Mt + F + Om)] + (29) Sin[ 2. Mt] + (29) Sin[Mt + 2. (-d + F + Om)] + (26) Sin[ 2. F] + (-22) Sin[ 2. (-d + F)] + (21) Sin[-Mt + Om + 2. F] + (17 - 0.1 T) Sin[ 2. M] + (16) Sin[2. d - Mt + Om] + (-16 + 0.1 T) Sin[ 2. (-d + M + F + Om)] + (-15) Sin[ M + Om] + (-13) Sin[-2. d + Mt + Om] + (-12) Sin[-M + Om] + (11) Sin[2. (Mt - F)] + (-10) Sin[ 2. (d + F) - Mt + Om] + (-8) Sin[ 2. (d + F + Om) + Mt] + (7) Sin[ 2. (F + Om) + M] + (-7) Sin[-2. d + M + Mt] + (-7) Sin[-M + 2. (F + Om)] + (-7) Sin[2. (d + F) + Om] + (6) Sin[ 2. d + Mt] + (6) Sin[2. (-d + Mt + F + Om)] + (6) Sin[ 2. (-d + F) + Mt + Om] + (-6) Sin[2. (d - Mt) + Om] + (-6) Sin[ 2. d + Om] + (5) Sin[-M + Mt] + (-5) Sin[ 2. (F - d) + Om - M] + (-5) Sin[Om - 2. d] + (-5) Sin[ 2. (Mt + F) + Om] + (4) Sin[2. (Mt - d) + Om] + (4) Sin[ 2. (F - d) + M + Om] + (4) Sin[Mt - 2. F] + (-4) Sin[ Mt - d] + (-4) Sin[M - 2. d] + (-4) Sin[d] + (3) Sin[ Mt + 2. F] + (-3) Sin[2. (F + Om - Mt)] + (-3) Sin[ Mt - d - M] + (-3) Sin[M + Mt] + (-3) Sin[ Mt - M + 2. (F - Om)] + (-3) Sin[ 2. (d + F + Om) - M - Mt] + (-3) Sin[ 3. Mt + 2. (F + Om)] + (-3) Sin[2. (d + F + Om) - M]; (* % nutation in obliquity[0."0001] *) DObl = (92025 + 8.9 T) Cos[Om] + (5736 - 3.1 T) Cos[ 2. (-d + F + Om)] + (977 - 0.5 T) Cos[ 2. (F + Om)] + (-895 + 0.5 T) Cos[2. Om] + (54 - 0.1 T) Cos[ M] + (-7) Cos[Mt] + (224 - 0.6 T) Cos[ 2. (-d + F + Om) + M] + (200) Cos[ 2. F + Om] + (129 - 0.1 T) Cos[ Mt + 2. (F + Om)] + (-95 + 0.3 T) Cos[-M + 2. (-d + F + Om)] + (-70) Cos[ 2. (-d + F) + Om] + (-53) Cos[-Mt + 2. (F + Om)] + (-33) Cos[ Mt + Om] + (26) Cos[-Mt + 2. (d + F + Om)] + (32) Cos[-Mt + Om] + (27) Cos[Mt + 2. F + Om] + (-24) Cos[ Om + 2. (-Mt + F)] + (16) Cos[2. (d + F + Om)] + (13) Cos[ 2. (Mt + F + Om)] + (-12) Cos[ Mt + 2. (-d + F + Om)] + (-10) Cos[-Mt + Om + 2 F] + (-8) Cos[ 2. d - Mt + Om] + (7) Cos[2. (-d + M + F + Om)] + (9) Cos[ M + Om] + (7) Cos[-2. d + Mt + Om] + (6) Cos[-M + Om] + (5) Cos[2. (d + F) - Mt + Om] + (3) Cos[ 2. (d + F + Om) + Mt] + (-3) Cos[ 2. (F + Om) + M] + (3) Cos[-M + 2. (F + Om)] + (3) Cos[ 2. (d + F) + Om] + (-3) Cos[2. (-d + Mt + F + Om)] + (-3) Cos[ 2. (-d + F) + Mt + Om] + (3) Cos[2. (d - Mt) + Om] + (3) Cos[ 2. d + Om] + (3) Cos[2. (F - d) + Om - M] + (3) Cos[ Om - 2. d] + (3) Cos[2. (Mt + F) + Om]; (* convert to radians *) DLon = DLon*0.0001/(3600.*rad); DObl = DObl*0.0001/(3600.*rad); Nut = {DLon, DObl} ]; Options[Obliquity] = {FullPrecision -> True}; SetAttributes[Obliquity, Listable]; Obliquity[{y_, m_, d_, h_: 0, min_: 0, sec_: 0}, opts : OptionsPattern[]] := Obliquity[JulianDate[{y, m, d, h, min, sec}], opts] Obliquity[JD_?NumberQ, opts : OptionsPattern[]] := Module[{T, Obl, U, fullPrecision}, fullPrecision = OptionValue[FullPrecision]; T = (JD - 2451545.0)/36525.0; U = T/100; If[fullPrecision, Obl = 23.44484666666667 - (4680.93 U - 1.55 U^2 + 1999.25 U^3 - 51.38 U^4 - 249.67 U^5 - 39.05 U^6 + 7.12 U^7 + 27.87 U^8 + 5.79 U^9 + 2.45 U^10)/3600., Obl = 23.439291 - 0.0130042 T - 0.00000016 T^2 + 0.000000504 T^3;]; Obl °]  As a small test, the same date as above: JulianDate[{1957, 10, 4}] >> 2.43612*10^6  Edit: Using the routine FromJulianDate, written by Peter Breitfeld, GST, GMST and GAST can be computed using the following routines:  FromJulianDay::usage = "FromJulianDay[JulianDay] converts a Julian day to a date in \Gregorian calendar."; GST::usage = "GST calculates the sidereal time at the meridian of Greenwich, at\ \!$$\*SuperscriptBox[\"0\", \"h\"]$$ Universal Time of a given date."; GMST::usage = "GMST[date_List] calculates the Greenwich Mean Sidereal Time (GMST)\ for a given date. The result is given in Degree."; GAST::usage = "GAST[date_List] calculates the Greenwich apparent sidereal time (GAST) for a given date.The result is given in Degree. The low precision IAU 1980 nutation algorithm has been used."; FromJulianDay[JD_] := Module[{y, m, d, II, FF, AA, BB, CC, DD, EE, GG, hour, min, sec}, (* Inspired by: *)(* brfAstro.m-- Paket zur Implementierung von Funktionen der Astrometrie Autor: Peter Breitfeld *) II = IntegerPart[JD + 1/2]; FF = FractionalPart[JD + 1/2]; If[II > 2299160, AA = IntegerPart[(II - 1867216.25)/36524.25]; BB = II + 1 + AA - IntegerPart[AA/4], BB = II]; CC = BB + 1524; DD = IntegerPart[(CC - 122.1)/365.25]; EE = IntegerPart[365.25 DD]; GG = IntegerPart[(CC - EE)/30.6001]; d = CC - EE + FF - IntegerPart[30.6001 GG]; If[GG < 13.5, m = GG - 1, m = GG - 13]; If[m > 2.5, y = DD - 4716, y = DD - 4715]; hour = IntegerPart[24 FractionalPart[d]]; min = IntegerPart[60 FractionalPart[24 FractionalPart[d]]]; sec = 60 FractionalPart[60 FractionalPart[24 FractionalPart[d]]]; {y, m, IntegerPart[d], hour, min, sec}]; GST[{y0_, m0_, d0_, h0_, min0_, s0_}] := (JulianDate[{y0, m0, d0, h0, min0, s0}] - 2451545.0)/36525 GST[JD_?NumberQ] := GST[FromJulianDay[JD]] GST[{y0_, m0_, d0_}] := (JulianDate[{y0, m0, d0, 0, 0, 0}] - 2451545.0)/36525 GMST[{y_, m_, d_}] := GMST[{y, m, d, 12, 0, 0}] GMST[{y_, m_, d_, h_, min_, sec_}] := Module[{TUT1, gmst}, TUT1 = (JulianDate[{y, m, d, h, min, sec}] - 2451545.)/36525; gmst = 67310.54841 (* sec *)+ (876600 3600 (* hours *)+ 8640184.812866 (* secs *)) TUT1 + 0.093104 TUT1^2 - 6.2 10^-6 TUT1^3; gmst = If[gmst < 0, Mod[gmst, -86400.], Mod[gmst, 86400.]]; (* convert gmst in sec -> degree 1s = 1/240 Deg *) gmst = gmst/240.; If[gmst < 0, 360 + gmst, 360 - gmst]] GAST[{y_, m_, d_}] := GAST[{y, m, d, 12, 0, 0}] GAST[{y_, m_, d_, h_, min_, sec_}] := Module[{em, T, L, Lprime, Ω, Δψ, Δε}, (* This function calculates the apparent Greenwich sidereal time using the first few terms of the IAU 1980 nutation algorithm.*) (* em = mean obliquity of the ecliptic *) T = (JulianDate[{y, m, d, h, min, sec}] - 2451545.)/36525; em = FromDMS[{23, 26, 21.448}] - FromDMS[{0, 0, 46.8150}] T - FromDMS[{0, 0, 0.00059}] T^2 + FromDMS[{0, 0, 0.001813}] T^3; L = 280.4665 + 36000.7698 T; Lprime = 218.3165 + 481267.8813 T; Ω = 125.04452 - 1934.136261 T; (* These three are in Degree *) Δψ = -17.20 Sin[Ω °] - 1.32 Sin[2 L Degree] - 0.23 Sin[2 Lprime °] + 0.21 Sin[2 Ω °]; Δε = 9.20 Cos[Ω °] + 0.57 Cos[2 L °] + 0.10 Cos[2 Lprime °] - 0.09 Cos[2 Ω °]; (* these two are in arc seconds *) GMST[{y, m, d, h, min, sec}] + (Δψ/3600.) Cos[(em + Δε/3600.) °]]  • For getting local sidereal time, the system variable $GeoLocation is quite useful... Sep 29, 2012 at 0:17

I think CalendarConvert, which is new in version 10, might help:

CalendarConvert[DateObject[{2014, 10, 16, 14, 10, 4}, TimeZone -> 0], "Julian"]
`