Moon eclipses calculator

An AppleScript script predicting Moon eclipses on a wide range of years.
The script provided is based on Jean MEEUS works published in his book:
Astronomical Formulae for Calculators © 1979 Willmann-Bell Editions.

Results are given with the date and time of the eclipse event, its type and magnitude.
The average precision for eclipses time calculated with this program is 2 minutes, maximum imprecision could reach 5 minutes.
Results will look like to this:
{“2017 08 07”, “Partial E. @ UTC 18 h 19 min / M: 0,256”}
Where date is YYYY MM dd; the type of eclipse can be “Penumbral”, “Partial” or “Total” ; “M” magnitude is, here, relative to the only Earth’s deep shadow (for partial and total eclipse), not its whole extent including the penumbra cone.
The time represents the moment of eclipse at its maximum.

Keep in mind that this script will give only a fast and rather precise reminder for interesting lunar event, not a reference tool.
For very accurate and extended informations about eclipses, please look at the master work by Fred ESPENAK on NASA site:
http://eclipse.gsfc.nasa.gov/lunar.html

Script tested and running OK from Lion up to El Capitan.
The script is ready to run in Applescript Editor, no need of Satimage or such “scripting addition”. Nonetheless, if you have a trigonometric osax installed, you can use it by getting rid of the pipes around |sin| and |cos| in formulae and commenting out or deleting the two functions provided at the bottom of the script.
Feel free to point out and correct any error encountered in the script or it’s results.

Enjoy!


(*Written by Joseph25b / oct 2015. Moon's eclipses calculator based on Jean MEEUS algorithms*)

set allMoonEclipses to yearMoonEclipses(2013, 2018) --Here, set the range of needed years

on yearMoonEclipses(Year1, Year2)
	set Refdate to my makeASDateFrom(1900, 1, 1, 6, 13, 26) --(1st january 1900 06:13:26 AM):DO NOT modify the Refdate!!
	set rad to pi / 180
	set moonEclipsesList to {}
	(*After Jean Meeus's algorithms published in "Astronomical Formulae for Calculators" - 1979*)
	repeat with aYear from Year1 to Year2
		set JJD to {}
		set iQ0 to 0
		set floorDate to my JJulien(31.99999, 12, (aYear - 1))
		set ceilingDate to my JJulien(1, 1, aYear + 1)
		set calendarCible to my makeASDateFrom(aYear, 1, 1, 0, 0, 0) -- (1st january aYear 0:0:0:)
		set Ktest to ((calendarCible - Refdate) / 29.53058868) / 86400
		set entierK to Ktest div 1
		set resteK to Ktest mod 1
		if resteK <= 0.5 then
			set K1 to entierK + 0.5
		else
			set K1 to entierK + 1.5
		end if
		repeat with i from 0 to 12 -- 
			set typeEclipse to ""
			set K to K1 + i
			set T to K / 1236.85
			set T2 to T * T
			set T3 to T2 * T
			set JJm to (2.41502075933E+6 + 29.53058832 * K + 1.337E-4 * T2 - 1.5E-7 * T3 + 3.3E-4 * (|sin|(rad * (166.56 + 132.87 * T - 0.009173 * T2))))
			set M to (359.2242 + 29.10535608 * K - 3.33E-5 * T2 - 3.47E-6 * T3)
			set M to (M mod 360) * rad
			set N to (306.0253 + 385.81691806 * K + 0.0107306 * T2 + 1.236E-5 * T3)
			set N to (N mod 360) * rad
			set F to (21.2964 + 390.67050646 * K - 0.0016528 * T2 - 2.39E-6 * T3)
			set F to (F mod 360) * rad
			set testEclipse to (|sin|(F))
			if testEclipse > -0.3584 and testEclipse < 0.3584 then -- If it can ever be an eclipse, the Moon's latitude must be under 21° modulo 180°
				set corr to ((0.1734 - 3.93E-4 * T) * (|sin|(M))) + (0.0021 * (|sin|(2 * M))) - (0.4068 * (|sin|(N))) + (0.0161 * (|sin|(2 * N))) - (0.0051 * (|sin|(M + N))) - (0.0074 * (|sin|(M - N))) - (0.0104 * (|sin|(2 * F)))
				set JJ to (JJm + corr) --the date and time of the maximum of eclipse, if any...
				--log JJ
				set S to my calculationOfS(M, N)
				set C to my calculationOfC(M, N, F)
				set Mu to my calculationOfMu(M, N)
				set gama to (S * (|sin|(F)) + C * (|cos|(F)))
				if gama < 0 then set gama to (-1 * gama)
				set grandeurPenombre to (1.5572 + Mu - gama) / 0.545
				set Magn to ((1.0129 - Mu - gama) / 0.545) + 0.009
				set Magn to ((Magn * 1000) div 1) / 1000
				if grandeurPenombre > 0 then -- a negative value means that there is no eclipse
					set typeEclipse to 1 -- at least a penumbral eclipse
					if Magn >= 1 then --Total Moon's eclipse in Earth's shadow cone 
						set typeEclipse to 3
					else if Magn > 0 then -- only a part of the Moon shadowed by Earth's shadow cone
						set typeEclipse to 2
					end if
					if JJ >= ceilingDate then exit repeat -- avoid last item's belonging to next year
					if JJ <= floorDate then "nothing" -- avoid first item's belonging to previous year
					if JJ > floorDate then copy {JJ, typeEclipse, Magn} as list to end of JJD -- lets'go
				end if
			end if
		end repeat
		repeat with nbEclipse from 1 to (count items in JJD)
			set tempEclipseDate to my (JJUlToCalendar(JJD's item nbEclipse)) as list
			set moonEclipsesList to moonEclipsesList & (tempEclipseDate)
		end repeat
	end repeat
	return moonEclipsesList
end yearMoonEclipses

on calculationOfS(M, N) -- part of J Meeus algorithms
	set tempVal to 5.19595 - (0.0048 * (|cos|(M))) + (0.002 * (|cos|(2 * M))) - (0.3283 * (|cos|(N))) - (0.006 * (|cos|(M + N))) + (0.0041 * (|cos|(M - N)))
	return tempVal
end calculationOfS

on calculationOfC(M, N, F) -- part of J Meeus algorithms
	set tempVal to (0.207 * (|sin|(M))) + (0.0024 * (|sin|(2 * M))) - (0.039 * (|sin|(N))) + (0.0115 * (|sin|(2 * N))) - (0.0073 * (|sin|(M + N))) - (0.0067 * (|sin|(M - N))) + (0.0117 * (|sin|(2 * F)))
	return tempVal
end calculationOfC

on calculationOfMu(M, N) -- part of J Meeus algorithms
	set tempVal to 0.0059 + (0.0046 * (|cos|(M))) - (0.0182 * (|cos|(N))) + (4.0E-4 * (|cos|(2 * N))) - (5.0E-4 * (|cos|(M + N)))
	return tempVal
end calculationOfMu

on JJUlToCalendar(JJeclipse) --  after J Meeus's algorithms
	set tEclipses to {"Penumbral E. ", "Partial E. ", "Total E. "}
	set eclipsindex to item 2 of JJeclipse
	set typEclp to (tEclipses's item eclipsindex)
	set magnit to "M: " & JJeclipse's item 3
	set JJul to (item 1 of JJeclipse) + 0.5
	set Z to (JJul div 1)
	set HMS to (JJul mod 1) * 86400
	set Fh to HMS div 3600
	set Fh to nb00j(Fh)
	set Fmn to (HMS mod 3600) div 60
	set Fmn to nb00j(Fmn)
	if Z >= 2299161 then -- Gregorian calendar
		copy ((Z - 1.86721625E+6) / 3.652425E+4) div 1 to alpha
		copy Z + 1 + alpha - ((alpha / 4) div 1) to A
	else -- Julian calendar
		copy Z to A
	end if
	copy A + 1524 to B
	copy ((B - 122.1) / 365.25) div 1 to C
	copy (365.25 * C) div 1 to D
	copy ((B - D) / 30.6001) div 1 to E
	copy B - D - ((30.6001 * E) div 1) to DD
	set DD to nb00j(DD)
	copy E - 1 to MM
	if E > 13.5 then copy E - 13 to MM
	copy nb00j(MM) to MMft
	copy C - 4716 to YYYY
	if MM < 2.5 then copy C - 4715 to YYYY --Year
	set ref_J_Eclipse to (YYYY & " " & MMft & " " & DD) as string
	return {((ref_J_Eclipse) as list) & (typEclp & "@ UTC " & Fh & " h " & Fmn & " min / " & magnit)} as list
end JJUlToCalendar

on JJulien(DD, MM, YY) --  J Meeus's algorithm
	copy YY + (MM / 100) + (DD / 10000) to Gregor
	if MM < 3 then
		copy MM + 12 to MM
		copy YY - 1 to YY
	end if
	copy (365.25 * YY) div 1 + (30.6001 * (MM + 1)) div 1 + DD + (1.7209945E+6) to JourJulien
	if Gregor > 1582.1015 then -- Step from julian to gregorian calendar
		copy (YY div 100) to K
		set JourJulien to JourJulien + (2 - K) + (K div 4)
	end if
	return JourJulien
end JJulien

--*********** Annexes handlers *******************************************
(*NB: Rather than unpredictable result from AS when making a "date" from a string, Shane Stanley suggests this function below*)
on makeASDateFrom(theYear, theMonth, theDay, theHours, theMinutes, theSeconds)
	tell (current date) to set {theASDate, year, day, its month, day, time} to {it, theYear, 1, theMonth, theDay, theHours * hours + theMinutes * minutes + theSeconds}
	return theASDate
end makeASDateFrom

on nb00j(val) -- returning two digits dates as string
	set temp to (100 + val) --as string
	return items 2 thru -1 of (temp as string) as text
end nb00j

on Bissextile(YY)
	if YY mod 4 is not  0 then return 0
	if YY mod 100 = 0 and YY mod 400 is not  0 then return 0
	return 1
end Bissextile

on |sin|(X) -- x en radian
	set loopList to {3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27} --enough iterations
	set X to (X as real) mod (2 * pi)
	set snus to X as real
	set Z to 1
	set X2 to X * X
	repeat with j in loopList
		set Z to Z * (j - 1) * j --factorial loop
		set X to X * X2 -- increment x^3, x^5, x^7, etc.
		set pow to (j / 2) - 0.5
		set snus to snus + ((-1) ^ pow) * (X / Z)
	end repeat
	set ctrl to snus
	if ctrl < 0 then set ctrl to ctrl * (-1)
	if ctrl < 9.999999999999E-14 then set snus to 0.0
	return snus
end |sin|

on |cos|(X) -- x en radian
	set loopList to {2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26} --enough iterations
	set X to (X as real) mod (2 * pi)
	set cosnus to 1 as real
	set Z to 1
	set X2 to X * X
	set X to 1
	repeat with j in loopList
		set X to X * X2 -- increment x^2, x^4, x^6, etc.
		set Z to Z * (j - 1) * j --factorial loop
		set pow to (j / 2)
		set cosnus to cosnus + ((-1) ^ pow) * (X / Z)
	end repeat
	set ctrl to cosnus
	if ctrl < 0 then set ctrl to ctrl * (-1)
	if ctrl < 9.9999999999E-14 then set cosnus to 0.0
	return cosnus
end |cos|