Sidereal time

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

Hello.

The Library for sidereal time has been updated, and it is completely accurate, with regards to test towards online sidereal calculator found here, as well as towards numerous other sources.