Hello.
This isn’t eveybody’s cup of tea, but if you want it, then getting the correct sidereal time is hard to get by using the formulaes that first pop up in google.
If you are curious about what it is for, then stars position are given with their right ascension to the vernal equinox in an angle, your apparent local siderealtime - right ascension, is the current hour-angle of the star.
# Version 2 Completely working and accurate library for sidereal time.
# GMST: Greenwhich sidereal time. (set time to 0 h 0 min, 0 sec to get GMST 0h)
# GAST: Greenwhichapparent time.
# LMST: Local MeanSidereal time.
# LAST: Local Apparent Sidereal Time.
property EpochJD2000 : 2.451545E+6
property JulianDaysInCentury : 36525
property degrad : pi / 180
property nutTbl : {¬
{0, 0, 0, 0, 1, -171996, -174.2, 92025, 8.9}, ¬
{-2, 0, 0, 2, 2, -13187, -1.6, 5736, -3.1}, ¬
{0, 0, 0, 2, 2, -2274, -0.2, 977, -0.5}, ¬
{0, 0, 0, 0, 2, 2062, 0.2, -895, 0.5}, ¬
{0, 1, 0, 0, 0, 1426, -3.4, 54, -0.1}, ¬
{0, 0, 1, 0, 0, 712, 0.1, -7, 0}, ¬
{-2, 1, 0, 2, 2, -517, 1.2, 224, -0.6}, ¬
{0, 0, 0, 2, 1, -386, -0.4, 200, 0}, ¬
{0, 0, 1, 2, 2, -301, 0, 129, -0.1}, ¬
{-2, -1, 0, 2, 2, 217, -0.5, -95, 0.3}, ¬
{-2, 0, 1, 0, 0, -158, 0, 0, 0}, ¬
{-2, 0, 0, 2, 1, 129, 0.1, -70, 0}, ¬
{0, 0, -1, 2, 2, 123, 0, -53, 0}, ¬
{2, 0, 0, 0, 0, 63, 0, 0, 0}, ¬
{0, 0, 1, 0, 1, 63, 0.1, -33, 0}, ¬
{2, 0, -1, 2, 2, -59, 0, 26, 0}, ¬
{0, 0, -1, 0, 1, -58, -0.1, 32, 0}, ¬
{0, 0, 1, 2, 1, -51, 0, 27, 0}, ¬
{-2, 0, 2, 0, 0, 48, 0, 0, 0}, ¬
{0, 0, -2, 2, 1, 46, 0, -24, 0}, ¬
{2, 0, 0, 2, 2, -38, 0, 16, 0}, ¬
{0, 0, 2, 2, 2, -31, 0, 13, 0}, ¬
{0, 0, 2, 0, 0, 29, 0, 0, 0}, ¬
{-2, 0, 1, 2, 2, 29, 0, -12, 0}, ¬
{0, 0, 0, 2, 0, 26, 0, 0, 0}, ¬
{-2, 0, 0, 2, 0, -22, 0, 0, 0}, ¬
{0, 0, -1, 2, 1, 21, 0, -10, 0}, ¬
{0, 2, 0, 0, 0, 17, -0.1, 0, 0}, ¬
{2, 0, -1, 0, 1, 16, 0, -8, 0}, ¬
{-2, 2, 0, 2, 2, -16, 0.1, 7, 0}, ¬
{0, 1, 0, 0, 1, -15, 0, 9, 0}, ¬
{-2, 0, 1, 0, 1, -13, 0, 7, 0}, ¬
{0, -1, 0, 0, 1, -12, 0, 6, 0}, ¬
{0, 0, 2, -2, 0, 11, 0, 0, 0}, ¬
{2, 0, -1, 2, 1, -10, 0, 5, 0}, ¬
{2, 0, 1, 2, 2, -8, 0, 3, 0}, ¬
{0, 1, 0, 2, 2, 7, 0, -3, 0}, ¬
{-2, 1, 1, 0, 0, -7, 0, 0, 0}, ¬
{0, -1, 0, 2, 2, -7, 0, 3, 0}, ¬
{2, 0, 0, 2, 1, -7, 0, 3, 0}, ¬
{2, 0, 1, 0, 0, 6, 0, 0, 0}, ¬
{-2, 0, 2, 2, 2, 6, 0, -3, 0}, ¬
{-2, 0, 1, 2, 1, 6, 0, -3, 0}, ¬
{2, 0, -2, 0, 1, -6, 0, 3, 0}, ¬
{2, 0, 0, 0, 1, -6, 0, 3, 0}, ¬
{0, -1, 1, 0, 0, 5, 0, 0, 0}, ¬
{-2, -1, 0, 2, 1, -5, 0, 3, 0}, ¬
{-2, 0, 0, 0, 1, -5, 0, 3, 0}, ¬
{0, 0, 2, 2, 1, -5, 0, 3, 0}, ¬
{-2, 0, 2, 0, 1, 4, 0, 0, 0}, ¬
{-2, 1, 0, 2, 1, 4, 0, 0, 0}, ¬
{0, 0, 1, -2, 0, 4, 0, 0, 0}, ¬
{-1, 0, 1, 0, 0, -4, 0, 0, 0}, ¬
{-2, 1, 0, 0, 0, -4, 0, 0, 0}, ¬
{1, 0, 0, 0, 0, -4, 0, 0, 0}, ¬
{0, 0, 1, 2, 0, 3, 0, 0, 0}, ¬
{0, 0, -2, 2, 2, -3, 0, 0, 0}, ¬
{-1, -1, 1, 0, 0, -3, 0, 0, 0}, ¬
{0, 1, 1, 0, 0, -3, 0, 0, 0}, ¬
{0, -1, 1, 2, 2, -3, 0, 0, 0}, ¬
{2, -1, -1, 2, 2, -3, 0, 0, 0}, ¬
{0, 0, 3, 2, 2, -3, 0, 0, 0}, ¬
{2, -1, 0, 2, 2, -3, 0, 0, 0} ¬
}
# References Jean Meeus Astronomical Algorithms chp 11
(*
on run
# set theDateGMT to my date_lib's set_date_long(2013, 3, 19, 12, 0, 0)
set theDateGMT to (current date) # + (time to GMT)
# set theDateGMT to my date_lib's set_date_long(1994, 6, 16, 18, 0, 0)
# set GMST_NOW to greenwichMeanSidereal_ASD for theDateGMT
tell _makeDateTimeSextet(theDateGMT)
set JD to (my _JDN(item 1, item 2, item 3, item 4, item 5, item 6))
end tell
# set jd0 to JDN_0h(JD)
# Calculations
# GMST
# Greenwhich mean sidereal is tested for both 0h, and other points
# in time.
set GMST_NOW to greenwichMeanSideral(JD)
set gmst_decimal to _mapWithinModulus(_rth(GMST_NOW), 24)
set GMST_hourTriple to _decDegToHexagesDeg(gmst_decimal)
# GAST
# Greenwhich appareant sidereal is tested for both 0h and other points.
set gast_decimal to gmst_decimal + (equationOfEquinoxes(JD) / 3600)
# the equationOfEquinoxes returns result in seconds - converts to degrees
set gast_decimal2 to _mapWithinModulus((_rth(greenwichApparentSidereal(JD))), 24)
set GAST_hourTriple to _decDegToHexagesDeg(gast_decimal)
set GAST_hourTriple2 to _decDegToHexagesDeg(gast_decimal2)
# LMST:
# Local Mean Sidereal Time.
set mdl to (_makeDecimalDegree for 8 by 46 against 57)
-- local mean sideereal
set LMST_NOW to (localMeanSiderealTime for JD against mdl) - ((_politicalTimeToGmtFromLocalTZInSeconds(theDateGMT)) * 1.00273790934 / 3600 * pi / 12)
set lmst_decimal to _mapWithinModulus(_rth(LMST_NOW), 24)
set LMST_hourTriple to _decDegToHexagesDeg(lmst_decimal)
# LAST:
# Local Apparent sidereal time. This is the same as the hour angle to the.
set LAST_NOW to (localApparentSiderealTime for JD against mdl) - ((_politicalTimeToGmtFromLocalTZInSeconds(theDateGMT)) * 1.00273790934 / 3600 * pi / 12)
set last_decimal to _mapWithinModulus(_rth(LAST_NOW), 24)
set LAST_hourtriple to _decDegToHexagesDeg(lmst_decimal)
end run
*)
# Mean sidereal time (GMST) in radians, it works for both 0h, and different
# times than 0h, by subtracting time, and adding back sidereal seconds
on greenwichMeanSideral(aJulianDate)
# http://www.celestrak.com/columns/v02n02/
local T_, UT_0h, UT_excess, GMST_deg, GMST_sidSec, GMST_rad, t_G
# 2000 January 1, 12h UT, is the Julian date 2451545.0:
set UT_excess to (aJulianDate + 0.5) mod 1
# we remove any excess time from 0hUT
set UT_0h to aJulianDate - UT_excess
set T_ to (UT_0h - 2.451545E+6) / 36525
set GMST_deg to 2.411054841E+4 + T_ * (8.640184812866E+6 + T_ * (0.093104 - T_ * 6.2E-6))
# Meeus (11.2) first term converted to seconds ! Only valid for 0h UT
set GMST_sidSec to (GMST_deg + 8.64E+4 * 1.00273790934 * UT_excess) mod 8.64E+4
set GMST_rad to 2 * pi * (GMST_sidSec / 8.64E+4)
# Better off normalizing degrees! Since we may add stuff!
if GMST_rad < 0 then set GMST_rad to (2 * pi) + GMST_rad
return GMST_rad
end greenwichMeanSideral
# Apparent sidereal time, in radians (GAST)
# Same as right ascension (RA)
on greenwichApparentSidereal(aJulianDate)
return _mapWithinModulus(greenwichMeanSideral(aJulianDate) ¬
+ ((equationOfEquinoxes(aJulianDate) / 3600) * pi / 12), (2 * pi))
end greenwichApparentSidereal
# returns the local mean sidereal time in hours with a decimal fraction
on localMeanSiderealTime for aJulianDate against aDecimalLongitude
if ((aDecimalLongitude < -180.0) or (aDecimalLongitude > 180)) then error "Longitude outside the range +/- 0,180Ë™." number 6001
return _mapWithinModulus((greenwichMeanSideral(aJulianDate) + (aDecimalLongitude * 1.00273790935 * pi / 180)), (2 * pi))
end localMeanSiderealTime
# returns the apparent local sidereal time in hours with a decimal fraction
on localApparentSiderealTime for aJulianDate against aDecimalLongitude
if ((aDecimalLongitude < -180.0) or (aDecimalLongitude > 180)) then error "Longitude outside the range +/- 0,180Ë™." number 6001
return _mapWithinModulus((greenwichApparentSidereal(aJulianDate) + (aDecimalLongitude * 1.00273790935 * pi / 180)), (2 * pi))
end localApparentSiderealTime
on greenwichMeanSidereal_0hASD for anASDate
local JD, GMST_0H
# Find time to gmt for this date ASDates are assumed to be adjusted up front?
tell (my date_lib's makeDateTimeSextet(anASDate))
set JD to (my JDN(item 1, item 2, item 3, 0, 0, 0))
end tell
return (greenwichMeanSideral(JD))
end greenwichMeanSidereal_0hASD
on greenwichMeanSidereal_ASD for anASDate
local JD, GMST_0H
# Find time to gmt for this date ASDates are assumed to be adjusted up front?
tell (my date_lib's makeDateTimeSextet(anASDate))
set JD to (my JDN(item 1, item 2, item 3, item 4, item 5, item 6))
end tell
return (greenwichMeanSideral(JD))
end greenwichMeanSidereal_ASD
# What is a DateTimeSextet?
# A DateTimeSextet is a list that consists of { yr,mnth,dy,hr,min,sec }
# set msx to makeDateTimeSextet from (current date)
on _makeDateTimeSextet(anASDate)
tell (anASDate)
# yyyy.mm.dd hh:mm:ss
return ({year of it as integer, month of it as integer, day of it as number, hours of it, minutes of it, seconds of it})
end tell
end _makeDateTimeSextet
on _JDN(_year, _month, _day, _hour, _min, _sec)
# From Wikipedia
# The Julian Date (JD) of any instant is the Julian day number
# for the preceding noon plus the fraction of the day since
# that instant. Julian Dates are expressed as a Julian day
# number with a decimal fraction added. For example, The Julian
# Date for 00:30:00.0 1 January 2013 is 2456293.520833.
local a, y, m, _JDN, time_fraction
set a to ((14 - _month) div 12)
set y to _year + 4800 - a
set m to _month + (12 * a) - 3
# Gregorian Calendar starts October 15, 1582
if _year > 1582 then
set _JDN to _day + ((153 * m + 2) div 5) + 365 * y + (y div 4) - y div 100 + y div 400 - 32045
else if _year = 1582 and _month > 10 then
set _JDN to _day + ((153 * m + 2) div 5) + 365 * y + (y div 4) - y div 100 + y div 400 - 32045
else if _year = 1582 and _month = 10 and _day ≥ 15 then
set _JDN to _day + ((153 * m + 2) div 5) + 365 * y + (y div 4) - y div 100 + y div 400 - 32045
else
# Julian calendar before oct 15, 1582
set _JDN to _day + ((153 * m + 2) div 5) + 365 * y + (y div 4) - 32083
end if
set _hour to (_hour - 12) / 24
if _min ≠0 then
set _min to _min / 1440 # minutes per day
end if
if _sec ≠0 then
set _sec to _sec / 86400 # seconds per day
end if
set _JDN to _JDN + _hour + _min + _sec
return _JDN
end _JDN
# The range should usually be 24 or 360
on _mapWithinModulus(aDecNumber, aModulus)
set aDecNumber to aDecNumber mod aModulus
if aDecNumber < 0 then set aDecNumber to aDecNumber + aModulus
return aDecNumber
end _mapWithinModulus
# converts a decimal degree number, that is a degree, with
# minutes and seconds as a decimal fraction into a real
# hexagesimal number with separate terms for degrees,
# minutes and seconds
-- Converts from the decimalDegree coordinat system (58.4603)
-- to DMS ( 58 27 37 )
on _decDegToHexagesDeg(decimalDegrees)
# new version, faster with increased precision
local tempvar, hexagesDegs, hexagesMins, hexagesSecs
set hexagesDegs to decimalDegrees div 1
set tempvar to (decimalDegrees mod 1) * 3600
# have to check for sign, and negate if have degrees,
# and degrees are negative.
if tempvar < 0 then
if hexagesDegs ≠0 then
set tempvar to tempvar * -1
end if
end if
set hexagesMins to tempvar div 60
# this sign always correct
set hexagesSecs to tempvar mod 60
# have to adjus if neg, and has degrees
if hexagesSecs < 0 then
if hexagesMins ≠0 then
set hexagesSecs to hexagesSecs * -1
end if
end if
return {hexagesDegs, hexagesMins, hexagesSecs}
end _decDegToHexagesDeg
on _time_TinJulianCenturies by JDE
# Meeus (11.1)
return ((JDE - (my EpochJD2000)) / (my JulianDaysInCentury))
end _time_TinJulianCenturies
on _politicalTimeToGmtFromLocalTZInSeconds(anASDate)
return (anASDate - _TZtoTZ(anASDate, (do shell script "readlink /etc/localtime |sed -n 's_/[^/]*/[^/]*/[^/]*/\\([^/]*\\)/\\([^/]*\\).*_\\1/\\2_p'"), "GMT"))
end _politicalTimeToGmtFromLocalTZInSeconds
on _TZtoTZ(TZ1date, TZ1, TZ2)
# Nigel Garvey
return (do shell script ("eraTime=$(TZ=" & TZ1 & " date -jf '%Y-%m-%dT%H:%M:%S' '" & (TZ1date as «class isot» as string) & "' '+%s') ; TZ=" & TZ2 & " date -r \"$eraTime\" '+%Y-%m-%dT%H:%M:%S'") as «class isot») as date
end _TZtoTZ
# radians converted to hours
on _rth(r)
return (r * 180 / pi / 15)
end _rth
to _makeDecimalDegree for tDegrees by tMinutes against tSeconds
return (tDegrees + ((tMinutes + (tSeconds / 60)) / 60))
end _makeDecimalDegree
on equationOfEquinoxes(JD)
# requires JulianDateLib that can be found at:
local T_, D_, M_, Mprime, F_, Omega, deltaRho, deltaEpsilon, argInRads, U_, epsilon0, Epsilon
set T_ to _time_TinJulianCenturies by JD
# set T to (JDE - EpochJD2000) / JulianDaysInCentury
set T_squared to T_ * T_
set T_cubed to T_squared * T_
# set d to (meanElongationOfMoonFromSun for T)
set D_ to _mapWithinModulus((297.85036 + 4.4526711148E+5 * T_ - 0.0019142 * T_squared + T_cubed / 189474), 360)
# set m to (meanAnamolyOfSun for T)
set M_ to _mapWithinModulus((357.52772 + 3.599905034E+4 * T_ - 1.603E-4 * T_squared - T_cubed / 300000), 360)
# set M_mark to (meanAnaomlyOfMoon for T)
set Mprime to _mapWithinModulus((134.96298 + 4.77198867398E+5 * T_ + 0.0086972 * T_squared + T_cubed / 56250), 360)
# set f to (moonsArgOfLatitude for T)
set F_ to _mapWithinModulus((93.27191 + 4.83202017538E+5 * T_ - 0.0036825 * T_squared + T_cubed / 327270), 360)
# set Omega to (moonsLongitudeForAscendingNodeOfMoonsMeanOrbitOnEclipticFromMeanEquinoxOfDate for T)
set Omega to _mapWithinModulus((125.04452 - 1934.136261 * T_ + 0.0020708 * T_squared + T_cubed / 450000), 360)
# deltaRho is the nutation in longitude deltaEpsilon is the nutation in obliquity.
set {deltaRho, deltaEpsilon} to {0, 0}
repeat with i from 1 to (length of nutTbl)
set argInRads to ((item 1 of item i of nutTbl) * D_ + (item 2 of item i of nutTbl) * M_ ¬
+ (item 3 of item i of nutTbl) * Mprime + (item 4 of item i of nutTbl) * F_ ¬
+ (item 5 of item i of nutTbl) * Omega) * degrad
set deltaRho to deltaRho + ((item 6 of item i of nutTbl) + (item 7 of item i of nutTbl) * T_) ¬
* (sin argInRads) * 1.0E-4
set deltaEpsilon to deltaEpsilon + ((item 8 of item i of nutTbl) + (item 9 of item i of nutTbl) * T_) ¬
* (cos argInRads) * 1.0E-4
end repeat
# calculating Epsilon_0 still Meus page 135 (21.2)
# set epsilon0 to 23.439291111111 - 0.013004166667 * T - (1.638888888889E-7 * T_squared) + (5.036111111111E-7 * T_cubed)
# trying with new version of epsilon0
set U_ to (JD - 2451545) / 3652500
set epsilon0 to (23.439291111111 ¬
+ (U_ * (-1.300258333333 + (U_ * (-4.305555555556E-4 + (U_ * (0.555347222222 + ¬
(U_ * (-0.014272222222 + (U_ * (0.069352777778 + (U_ * (-0.010847222222 + ¬
(U_ * (0.001977777778 + (U_ * (0.007741666667 + (U_ * (0.001608333333 + ¬
(6.805555555556E-4 * U_))))))))))))))))))))
set Epsilon to epsilon0 + (deltaEpsilon / 3600)
return (deltaRho * (cos (Epsilon * degrad)) / 15)
end equationOfEquinoxes