Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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)?

share|improve this question
2  
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 Sep 28 '12 at 19:24
    
@cormullion +1 for the reference. Great page that I also used for many of my implementations. –  Markus Roellig Sep 28 '12 at 22:42

2 Answers 2

up vote 10 down vote accepted

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}
share|improve this answer
1  
Great! A pity that I can't give you another +1 for each edit :) –  Markus Roellig Sep 29 '12 at 17:34
2  
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. –  J. M. Oct 1 '12 at 1:51
    
BTW: Still mining this for useful stuff. +2, if I could! –  raxacoricofallapatorius Nov 20 '12 at 16:16
    
@rax, I'm glad to hear the effort I put into this was worthwhile. You're welcome. :) –  J. M. Nov 20 '12 at 16:22
1  
@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}]). –  raxacoricofallapatorius Dec 10 '12 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.) °]]
share|improve this answer
1  
For getting local sidereal time, the system variable $GeoLocation is quite useful... –  J. M. Sep 29 '12 at 0:17

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.