How do I compute sidereal time and JD?

Question

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.
@cormullion +1 for the reference. Great page that I also used for many of my implementations.

J. M. · Accepted Answer · 2012-10-01 01:45:55Z

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 :)
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.
@rax, I'm glad to hear the effort I put into this was worthwhile. You're welcome. :)
@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}]).

Markus Roellig · Answer 2 · 2012-09-30 14:30:51Z

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 Degree]

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 tusing the following routines:

 FromJulianDay::usage = "FromJulianDay[JulianDay] converts a Julian day 
                        to a date in \Gregorian calendar.";
 GST::usage = "GST calculates the sideral 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 Sideral 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, \[CapitalOmega], \[CapitalDelta]\[Psi], \[CapitalDelta]\
         \[Epsilon]},
      (*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;
   \[CapitalOmega] = 125.04452 - 1934.136261 T;
   (*These three are in Degree*)
   \[CapitalDelta]\[Psi] = -17.20 Sin[\[CapitalOmega] Degree] - 1.32 Sin[2 L  
    Degree] - 0.23 Sin[2 Lprime Degree] + 0.21 Sin[2 \[CapitalOmega] Degree];
   \[CapitalDelta]\[Epsilon] = 9.20 Cos[\[CapitalOmega] Degree] + 0.57 
     Cos[2 L Degree] + 0.10 Cos[2 Lprime Degree] - 0.09 
     Cos[2 \[CapitalOmega] Degree];
   (*these two are in arc seconds*)
   GMST[{y, m, d, h, min, sec}] + (\[CapitalDelta]\[Psi]/
   3600.) Cos[(em + \[CapitalDelta]\[Epsilon]/3600.) Degree]]

For getting local sidereal time, the system variable $GeoLocation is quite useful...

asked	7 months ago
viewed	399 times
active	7 months ago

How do I compute sidereal time and JD?

2 Answers

Your Answer

Not the answer you're looking for? Browse other questions tagged packages date-and-time astronomy or ask your own question.

Linked

How do I compute sidereal time and JD?

2 Answers

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged packages date-and-time astronomy or ask your own question.

Linked

Related