I'm trying to write a code to determine sunrise and sunset for the current date and location. I use LocateMe for the latter, and I call it with do shell script "~...LocateMe"
, then use the -f
function to get the latitude, longitude, and altitude. As the altitude is in meters, I multiply it by 3.28084 to get feet.
The remainder of the code simply uses the sunrise equation to calculate sunrise for that time and location. I have Satimage.osax to handle the trig functions, but as it calculates in radians while the above equation for some reason uses degrees, I'm careful to multiply the degrees by π/180 and the inverse trig formulae by 180/π.
The code outputs a Julian date, as it should. But it outputs sunrise for today, for instance, as July 18.914151470176876 2017 - just before midnight tonight. A few minutes off for a miscalculation in coordinates is one thing, but it can't be that far off. So where did I go wrong?
set n to ((current date) - (date "Saturday, January 1, 2000 at 12:00:00 PM")) / days
set n to n div 1
set longitude to do shell script "~/Library/Application\\ Support/LocateMe/LocateMe -f \"{LON}\""
set latitude to do shell script "~/Library/Application\\ Support/LocateMe/LocateMe -f \"{LAT}\""
set altitude to do shell script "~/Library/Application\\ Support/LocateMe/LocateMe -f \"{ALT}\""
set altitude to altitude * 3.28084
set Jstar to (n - (longitude / 360))
set M to ((357.5291 + (0.98560028 * Jstar)) mod 360)
set x to 1.9148 * (sin (M * pi / 180))
set y to 0.02 * (sin (2 * M * pi / 180))
set z to 3.0E-4 * (sin (3 * M * pi / 180))
set C to (x + y + z)
set lambda to ((M + C + 180 * 102.9372) mod 360)
set q to 2.4515455E+6
set r to (0.0053 * (sin (M*pi/180)))
set s to (-0.0069 * (sin (2 * lambda*pi/180)))
set Jtransit to (q + Jstar + r + s)
set d to ((asin (sin (lambda * pi / 180) * (sin (23.44 * pi / 180)))) * 180 / pi)
set numerator to (sin ((-0.83 - (1.15 * (sqrt (altitude))) / 60) * pi / 180) - (sin (latitude * pi / 180) * (sin (d * pi / 180))))
set denominator to (cos (latitude * pi / 180) * (cos (d * pi / 180)))
set HRA to ((acos (numerator / denominator)) * 180 / pi)
set sunrise to (Jtransit - (HRA / 360))
set sunset to (Jtransit + (HRA / 360))
After noting a location in which I wasn't so careful in my π/180's, the formula now consistently logs sunrise and sunset at roughly 19 hours later than it should be. So, for instance, sunrise this morning was 6:30 AM or so. The formula outputs tonight at 1:45 AM instead. I have confirmed that all of the calculations are correct; the problem is that it outputs the wrong Julian day.
EDIT: I realized that n
in the equation is today's day. My formula originally included the leftover fraction of the day. I added line 2 to account for that; it's perhaps a little weird, but that was the best I could think of.
That said, the formula now gives 7:30 PM for sunrise - still way off, but closer than before.