Hello.
This script presents a report, over the difference between the true noon and mean noon, or just a dialog if you only enter one date. The accuracy is within +/- 20 seconds of the numbers represented in the Wikipedia article about the equation of time.
All faults are mine, but it is unfair to Nigel to say that I wrote this, because this script would never have seen a shed of light, without his contributions, which are sustainable.
You’ll need Satimage.osax in order to run this.
property parent : AppleScript
# http://en.wikipedia.org/wiki/Equation_of_time
script equTimeLowAcc
property parent : AppleScript
property EpochJD2000 : 2.451545E+6
property JulianDaysInCentury : 36525
property ScriptTitle : "Equation of Time:"
property JD0 : («data isot303030312D31312D3234543132» as date) - ((«data isot343731342D31312D3234543132» as date) - ((«data isot303030312D31312D3234543132» as date) - 366 * days))
-- Nigel Garvey http://macscripter.net/viewtopic.php?pid=167391#p167391 post # 6
property Debug : true
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} ¬
}
on run
set {dateRangeStart, dateRangeEnd} to getDateRange((current date), "Month", ScriptTitle)
if ((dateRangeEnd - dateRangeStart) / days) < 1 then
set singular to true
else
set singular to false
set oldMonth to dateRangeStart's month
_MAKEDOC("Equation of Time in the period from " & dateRangeStart's short date string & " to " & dateRangeEnd's short date string & return, 500, 800)
end if
set JDT_cur to anasdatetojdn(dateRangeStart) div 1 + 0.5 # 0h
set JDT_end to anasdatetojdn(dateRangeEnd) div 1 + 0.5
repeat while JDT_cur ≤ JDT_end
local eqTime, trueNoon, h, m, s, eqtimestr
set eqTime to computeEquationOfTime(JDT_cur)
set trueNoon to 12 - eqTime + JDTdiffGmt(JDT_cur)
set {h, m, s} to decDegToHexagesDeg(trueNoon)
set noonstr to formattime(h, m, s)
set {h, m, s} to decDegToHexagesDeg(eqTime)
if eqTime < 0 then
set eqtimestr to "-" & pad(-m) & " minutes and " & pad(s div 1) & " seconds."
else
set eqtimestr to pad(m) & " minutes and " & pad(s div 1) & " seconds."
end if
if singular then
tell application (path to frontmost application as text)
activate
display dialog "Yield: circa " & eqtimestr & return & return & "The sun passes through the Meredian at " & noonstr & " when the Sun is in the Zenith." with title "Equation of Time:" buttons {"Ok"} default button 1 cancel button 1 with icon note
end tell
else
set {yr, mnt, _day} to gregorianCalendarTriplet from JDT_cur
if mnt ≠oldMonth then
_PRINT(linefeed)
set oldMonth to mnt
end if
printEquationOfTime(eqtimestr, noonstr, yr, mnt, _day)
end if
set JDT_cur to JDT_cur + 1
end repeat
end run
# returns the time difference between phystime to gmt,
# and civil time to gmt.
on JDTdiffGmt(JDT) # **
local pgmt, cgmt, diff
set pgmt to physOffsetToGmt()
set gregDt to GregorianDate for JDT
set cgmt to civilTimeToGmtFromLocalTZInSeconds(gregDt) / hours
if cgmt < 0 then
set diff to pgmt - cgmt
else
set diff to cgmt - pgmt
end if
return diff
end JDTdiffGmt
on physOffsetToGmt() # **
global localLongitude, localPhysOffsetToGMT
local A
try
set A to localPhysOffsetToGMT
on error
set localLongitude to longitudeAsInTimePreferences()
set localPhysOffsetToGMT to (offsetOfhourAngleOfLocalToGMTMeredian by localLongitude)
end try
return localPhysOffsetToGMT
end physOffsetToGmt
# set decimalAngleForMyLongitude to longitudeAsInTimePreferences()
# Returns the hourAngle offset from the Greenwhich Meredian, according to your Time Preferences.
# see: http://macscripter.net/viewtopic.php?pid=166581#p166581
# post 44
to longitudeAsInTimePreferences()
local LocalRegion, decimalDegree, theMinutes, theSeconds
try
tell (do shell script "ls -al /etc/localtime") to set LocalRegion to word -2 & "." & word -1
on error
error "longitudeAsInTimePreferences: something is wrong, do you have a position specified in Time Preferences?" number 5000
end try
tell (do shell script "sed -ne '/" & LocalRegion & "/ s/.*\\([-+][[:digit:]]*\\).*/\\1/p' </usr/share/zoneinfo/zone.tab")
try
if length of it is 6 then
set {decimalDegree, theMinutes} to {(text 1 thru 4 of it as number), (text 5 thru -1 of it as number)}
if theMinutes ≠0 then
set decimalDegree to decimalDegree + theMinutes / 60
end if
else # length of it is 8 Thanks to kel1 fo showing me.
set {decimalDegree, theMinutes, theSeconds} to {(text 1 thru 4 of it as number), (text 5 thru 6 of it as number), (text 7 thru -1 of it as number)}
if theMinutes ≠0 then
set theMinutes to theMinutes / 60
end if
if theSeconds ≠0 then
set theSeconds to theSeconds / 3600
end if
set decimalDegree to decimalDegree + theMinutes + theSeconds
end if
on error
error "longitudeAsInTimePreferences: something is wrong, do you have a /usr/share/zoneinfo/zone.tab file on your system?" number 5000
end try
end tell
return decimalDegree
end longitudeAsInTimePreferences
# Returns the phsyical timezone for a longitude
to offsetOfhourAngleOfLocalToGMTMeredian by decimalAngle
# see: http://macscripter.net/viewtopic.php?pid=166581#p166581
# post 43
# we remove any superfluos degrees first.
set decimalAngle to decimalAngle mod 360
# we make angles ≥ 180 into negatives
if decimalAngle ≥ 180 then
set decimalAngle to decimalAngle - 360
end if
local hourAngle
if decimalAngle < 0 then
set decimalAngle to (decimalAngle * -1) - 0.1
# we have to skew the numbers by -0.1 to keep
# boundary conditions
set hourAngle to (24 - ((decimalAngle + 7.5) div 15)) mod 24
else
set hourAngle to (decimalAngle + 7.5) div 15
end if
return hourAngle
end offsetOfhourAngleOfLocalToGMTMeredian
on TZtoTZ(TZ1date, TZ1, TZ2)
return (do shell script ("eraTime=$(TZ=" & TZ1 & " date -jf '%FT%T' '" & (TZ1date as «class isot» as string) & "' '+%s') ; TZ=" & TZ2 & " date -r \"$eraTime\" '+%FT%T'") as «class isot») as date
end TZtoTZ
on civilTimeToGmtFromLocalTZInSeconds(anASDate)
return (anASDate - TZtoTZ(anASDate, (do shell script "readlink /etc/localtime |sed -n 's_/[^/]*/[^/]*/[^/]*/\\([^/]*\\)/\\([^/]*\\).*_\\1/\\2_p'"), "GMT"))
end civilTimeToGmtFromLocalTZInSeconds
on GregorianDate for theJDN
-- Nigel Garvey
return JD0 + theJDN * days
end GregorianDate
on anasdatetojdn(anASDate)
tell (anASDate)
# yyyy.mm.dd hh:mm:ss
return (my JDN(year of it, (month of it as integer), day of it, hours of it, minutes of it, seconds of it))
end tell
end anasdatetojdn
on printEquationOfTime(eqtimestr, noonstr, yr, mnt, _day)
global us_citizen
local A
try
set A to us_citizen
on error
set us_citizen to getUsCitizenShip()
end try
if us_citizen then
_PRINT(pad(mnt) & "/" & pad(_day) & "/" & yr & ": " & eqtimestr & " The sun passes through the Meredian at " & noonstr & linefeed)
else
_PRINT(pad(_day) & "/" & pad(mnt) & "/" & yr & ": " & eqtimestr & " The sun passes through the Meredian at " & noonstr & linefeed)
end if
end printEquationOfTime
on pad(N)
-- 22. September 2013, improved pad, for the cases where n < 0 and n > 100.
-- Thanks to Yvan Koenig.
-- Fixed a bug
if N < -100 then
return N as text
else if N < 0 then
return ("-" & text -2 thru -1 of ("0" & ((N as number) * -1) as text))
-- error "The pad handler assumes positive values" number 5000
else if N < 100 then
return text 2 thru -1 of ((100 + N) as text)
else
return N as text
end if
end pad
on _PRINT(txt)
# Nigel Garvey made the original version.
local A
tell application "TextEdit"
if not my Debug then
activate
else
do shell script "open -b \"com.apple.textedit\""
end if
if ((count its documents) is 0) then make new document
set A to (get path of its front document)
try
if A is "" then set A to missing value
on error
set A to missing value
end try
if A is not missing value then make new document
make new paragraph at end of text of front document with data (txt & linefeed) with properties {font:"Andale Mono", size:12}
end tell
return true
end _PRINT
on _MAKEDOC(heading, _width, _height)
tell application "TextEdit"
tell (make new document at front)
set text of it to heading
set font of every character to "Andale Mono"
set size of every character to 13
end tell
tell its front window
# set {item 3 of its bounds, item 4 of its bounds} to {800, 300}
set _bounds to bounds
set item 3 of _bounds to (item 1 of _bounds) + _width
set item 4 of _bounds to (item 2 of _bounds) + _height
set bounds to _bounds
end tell
end tell
end _MAKEDOC
on formatDate(yr, mnt, _day)
global us_citizen
local A
try
set A to us_citizen
on error
set us_citizen to getUsCitizenShip()
end try
if us_citizen then
return (pad(mnt) & "/" & pad(_day) & "/" & yr)
else
return (pad(_day) & "/" & pad(mnt) & "/" & yr)
end if
end formatDate
to getUsCitizenShip()
try
set localeString to (do shell script "defaults read .GlobalPreferences AppleLocale 2>/dev/null || echo en_US |tr '
' ' '")
on error E number N
return true
end try
if localeString is "en_US" then return true
return false
end getUsCitizenShip
on formattime(_hr, _min, _sec)
return (pad((_hr div 1)) & ":" & pad(_min) & ":" & pad((_sec div 1)))
end formattime
# 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 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
on gregorianCalendarTriplet from JDN
# This is an algorithm by Richards to convert a Julian Day Number, J, to a date in the Gregorian calendar
# (proleptic, when applicable). Richards does not state which dates the algorithm is valid for.
# it seems to work for dates after, BUT NOT INCLUDING 15th of oct 1582.
# I have extended it to return a hour, min, sec triplet for the fractional part of the JDN.
local y, m, N, R, p, V, s, w, b, C, f, g, h, D, MT, yr, JD
set y to 4716
set j to 1401
set m to 2
set N to 12
set R to 4
set p to 1461
set V to 3
set u to 5
set s to 153
set w to 2
set b to 274277
set C to -38
set JD to JDN div 1
set f to JD + j + (((4 * JD + b) div 146097) * 3) div 4 + C
set E to R * f + V
set g to (E mod p) div R
set h to u * g + w
set D to (h mod s) div u + 1
set MT to (h div s + m) mod N + 1
set yr to E div p - y + (N + m - MT) div N
return {yr, MT, D}
end gregorianCalendarTriplet
-- Ask the user for the range of dates to be covered.
-- from a start date, up to, and including enddate.
on getDateRange(initialDate, initialRange, ScriptTitle)
-- Nigel Garvey http://macscripter.net/viewtopic.php?id=41273
-- customized a little by McUsr (All faults are mine!),
-- so it is able to deliver just one date denoting a range of 1 day.
-- There are preset ranges of Month, Week, and Year, based on the initalDate.
set anIcon to note
set theMessage to "Enter the required date range:"
set today to initialDate
repeat
if initialRange = "Month" then
set enddate to today + 30 * days
set day of enddate to day of today
set i to 0
repeat while month of enddate > ((month of today) + 1)
set month of enddate to ((month of enddate) - 1)
set day of enddate to ((day of enddate) - 1)
set i to i + 1
end repeat
if i > 0 then
set enddate to today + (30 - i) * days
end if
else if initialRange = "Week" then
set enddate to today + 6 * days
else if initialRange = "Year" then
set enddate to today + 364 * days
# We'll add one day for a leap year
local yr
if today's month ≤ 3 then
set yr to today's year
if ((yr mod 4) = 0) and (((yr mod 100) ≠0) or ((yr mod 400) = 0)) then set enddate to endate + 1 * days
else
set yr to enddate's year
if ((yr mod 4) = 0) and (((yr mod 100) ≠0) or ((yr mod 400) = 0)) then set enddate to endate + 1 * days
end if
end if
set d1 to today's short date string
set d2 to enddate's short date string
tell application (path to frontmost application as text)
set dateRange to text returned of (display dialog theMessage default answer d1 & " - " & d2 with title ScriptTitle buttons {"Quit", "Ok"} cancel button 1 default button 2 with icon anIcon)
end tell
try
set {tids, AppleScript's text item delimiters, gotOne} to {AppleScript's text item delimiters, "-", false}
set {dateRangeStart, gotOne} to {date (text item 1 of dateRange), true}
set dateRangeEnd to date (text item 2 of dateRange)
set AppleScript's text item delimiters to tids
exit repeat
on error
if gotOne and length of text items of dateRange = 1 then
set dateRangeEnd to date (text item 1 of dateRange)
set AppleScript's text item delimiters to tids
exit repeat
else if gotOne then
set today to dateRangeStart
else
set today to initialDate
end if
set AppleScript's text item delimiters to tids
if theMessage does not contain "wasn't" then
set theMessage to "The input wasn't valid!" & return & return & theMessage
set anIcon to caution
end if
end try
end repeat
set dateRangeEnd's time to days - 1 -- Sets the last date's time to 23:59:59, the last second of the range.
return {dateRangeStart, dateRangeEnd}
end getDateRange
# returns the number of Julian centuries that has passed
# since the JDE
on JulianCenturies by JDE
# Meeus (11.1)
return ((JDE - (my EpochJD2000)) / (my JulianDaysInCentury))
end JulianCenturies
on mapWithinModulus(aDecNumber, aModulus)
set aDecNumber to aDecNumber mod aModulus
if aDecNumber < 0 then set aDecNumber to aDecNumber + aModulus
return aDecNumber
end mapWithinModulus
on equationOfEquinoxes(JD)
# Returns correction in seconds of time
-- eqEq = returned. (Not necessarily used.)
# requires JulianDateLib that can be found at:
-- Epsilon0 = Mean ObliguityOfEcliptic.
-- deltaPsi =
-- deltaEpsilon =
-- Epsilon = Epsilon0 + deltaEpsilon = True Obliguity of Ecliptic.
# set { eqeq,epsilon0, deltaEpsilon, deltaPsi, Epsilon } to equationOfEquinoxes(myjde)
local T_, D_, M_, Mprime, F_, Omega, argInRads, U_, epsilon0, deltaEpsilon, deltaPsi, epsilon
set T_ to JulianCenturies 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)
# deltaPsi is the nutation in longitude deltaEpsilon is the nutation in obliquity.
set {deltaPsi, 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 deltaPsi to deltaPsi + ((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 {(deltaPsi * (cos (epsilon * degrad)) / 15), epsilon, epsilon0, deltaEpsilon, deltaPsi}
# Returns correction in seconds of time
end equationOfEquinoxes
# Meeus chapter 27.
# The Equation of time is the difference between apparent (true)
# it is the difference between the hour-angle of the true sun
# and the mean (fictious) sun.
# Equation of time can be positive or negative
# if E > 0 then true sun crosses observers mereidan before
# mean sun.
# This equation is necessary if you want to compute the height of
# sun while it is at zenith.
on computeEquationOfTime(JDT)
local dummy, epsilon, epsilon0, deltaPsi, epsilon
# set JDT to anasdatetojdn((current date)) div 1 + 0.5 # 0h
set JDT to JDT div 1 + 0.5
# set JDT to 2.4489085E+6
# Todo: adjust for diff phys time to gmt - civil time to gmt.
set T_ to JulianCenturies by JDT
set L0 to mapWithinModulus(geometricMeanLongitudeOfSun(T_), 360)
set m to mapWithinModulus(meanAnamolyOfSun(T_), 360)
set ecc to eccOfEarthOrbit(T_)
set equOfC to sunsEquationOfCenter(T_, m, L0)
set trueLongitude to mapWithinModulus((L0 + equOfC), 360)
# why true longitude?
set trueAnamoly to mapWithinModulus((m + equOfC), 360)
set R to sunsRadiusVector(ecc, trueAnamoly)
# Meeus (24.5)
set Omega to mapWithinModulus((125.04 - (T_ * 1934.136)), 360)
set apparentLong to (trueLongitude - 0.00569 - (0.00478 * (sin (Omega * pi / 180))))
# corrected for nutation, low accuracy.
# The suns right ascension, and declination, in equatorial coordinates
set {dummy, epsilon0, deltaEpsilon, deltaPsi, epsilon} to equationOfEquinoxes(JDT)
-- Epsilon0 = Mean ObliguityOfEcliptic.
-- deltaPsi =
-- deltaEpsilon =
-- Epsilon = Epsilon0 + deltaEpsilon = True Obliguity of Ecliptic.
-- eqEq = returned. (Not necessarily used.)
-- set ** meanObliguityOfEcliptic to (meanObliquityOfEcliptic_Laskar of AstronomyLib for JDT)
set meanObliguityOfEcliptic to epsilon0 + 0.00256 * (cos (Omega * pi / 180))
set Ra to mapWithinModulus(sunsRightAscension(meanObliguityOfEcliptic, apparentLong) * 180 / pi, 360)
# reworking before working backwards:
# 0ld call:
-- set ** epsilon to AstronomyLib's trueEpsilon(JDT)
-- set ** deltaPsi to (AstronomyLib's nutationInEclipticLongitude(JDT)) / 3600
# n3w call:
return (equationOftime(L0, Ra, (deltaPsi / 3600), epsilon) * 4) / 60
end computeEquationOfTime
on equationOftime(L0, alfa, deltaPsi, epsilon)
# L0 : suns mean longitude
# alfa : apparent right ascension of sun.
# taking into account the aberration and nutation.
# deltaRho : nutation in longitude.
# epsilon : Obliquity of the ecliptic.
# Lo, alfa and deltaRho is in degrees.
return (L0 - 0.057183 - alfa + deltaPsi * (cos epsilon))
# 0.057183 is the sum of mean value of aberration
end equationOftime
# set Decl to
# set sunsRA to cos
on geometricMeanLongitudeOfSun(T)
# Meeus (24.2)
return (280.46645 + (T * (3.600076983E+4 + (T * 3.032E-4))))
end geometricMeanLongitudeOfSun
on meanAnamolyOfSun(T)
# Meeus ( 24.3)
return (357.52919 + (T * (3.5999053E+4 + (T * (-1.559E-4 + (T * -4.8E-7))))))
end meanAnamolyOfSun
on eccOfEarthOrbit(T)
# Meeus (24.4)
return (0.016708617 + (T * (-4.2037E-5 - (T * 1.236E-7))))
end eccOfEarthOrbit
on sunsEquationOfCenter(T, m, L0)
return (((1.9146 - (T * (0.004817 - (T * 1.4E-5)))) * (sin (m * pi / 180))) + ((0.01993 - (T * 1.01E-4)) * (sin (2 * (m * pi / 180)))) + ((2.9E-4 * (sin (3 * (m * pi / 180))))))
end sunsEquationOfCenter
on sunsRadiusVector(eccentricity, trueAnamoly)
# Meeus (24.5)
return ((1.00000010018 * (1 - (eccentricity ^ 2))) / (1 + (eccentricity * (cos (trueAnamoly * pi / 180)))))
end sunsRadiusVector
on sunsRightAscension(trueObliguity, Longitude)
# (24.6)
return atan2 {((cos (trueObliguity * pi / 180)) * (sin (Longitude * pi / 180))), (cos (Longitude * pi / 180))}
end sunsRightAscension
end script
tell equTimeLowAcc to run