Accurate Table of Moonphases

Hello.

The script below produces accurate tables of Moon Phases, besides tell you the current phase of the Moon.

Enjoy! :slight_smile:

Edit

You need Satimage.osax in order to run this!

-- © 2013 McUsr, sustainable share of code kindly provided by Nigel Garvey.
-- It is in public domain Licenced under GLPGL 3.0, which means 
-- that if you use this then you must provide the source code, 
-- and licence your product too under GLPGL 3.0.

script MoonPhases
	property scriptTitle : "Moon Phases"
	
	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 movingStates : {¬
		"Waning Crescent", ¬
		"Waxing Crescent", ¬
		"Waxing Gibbous", ¬
		"Waning Gibbous"}
	
	property exactStates : {¬
		"New Moon", ¬
		"First Quarter", ¬
		"Full Moon", ¬
		"Third Quarter"}
	
	property aNewMoon : 1
	property tFirstQuarter : 2
	property aFullMoon : 3
	property tThirdQuarter : 4
	
	property monthNames : {¬
		"January", ¬
		"February", ¬
		"Mars", ¬
		"April", ¬
		"May", ¬
		"June", ¬
		"July", ¬
		"August", ¬
		"September", ¬
		"October", ¬
		"November", ¬
		"December"}
	
	
	on createMoonPhasesTable()
		local dyr, startY, endY, MoonPhases
		
		set {startY, endY} to getRangeOfYears((current date)'s year, "Decade", scriptTitle)
		if endY - startY > 0 then
			_MAKEDOC("Lunar phases from " & startY & " to " & endY & linefeed & linefeed & linefeed, 500, 800)
		else
			_MAKEDOC("Lunar phases  for " & startY & linefeed & linefeed & linefeed, 500, 800)
		end if
		
		repeat with dyr from startY to endY
			
			repeat with mnth from 1 to 12
				_PRINT(item mnth of monthNames & " " & dyr & linefeed)
				set MoonPhases to moonphasesForMonth(dyr, mnth)
				repeat with aPhase in MoonPhases
					local aJde, yr_, mnt_, dy_, hr_, min_, sec_
					
					set {yr_, mnt_, dy_} to gregorianCalendarTriplet from item 1 of aPhase
					set dateStr to formatDate(yr_, mnt_, dy_)
					set {hr_, min_, sec_} to timeTriplet from item 1 of aPhase
					set timeStr to formatClockTime(hr_, min_, sec_)
					_PRINT(dateStr & " : " & timeStr & tab & item (item 2 of aPhase) of exactStates & return)
				end repeat
				_PRINT(linefeed)
			end repeat
			_PRINT(linefeed)
		end repeat
	end createMoonPhasesTable
	
	
	on run
		global am_pm_clock, us_citizen -- "System variables" for configuring the locale. (Date/clock).
		-- resets whatever values, so that any script copying will work with the receivers locale.
		set am_pm_clock to undef()
		set us_citizen to undef()
		
		
		moonPhaseOfToday()
		createMoonPhasesTable()
	end run
	
	on undef()
		return
	end undef
	
	
	on moonPhaseOfToday()
		script o
			property L : {}
		end script
		local YR, mnth, dy, curJDe
		tell (current date)
			set curJDe to (JulianDayNumber of me for it) + 0.5 div 1
			set {YR, mnth, dy} to {its year, its month as integer, its day}
		end tell
		
		set o's L to moonphasesForMonth(YR, mnth)
		local phase
		set phase to ""
		repeat with i from 1 to (count o's L)
			if curJDe < (item 1 of item i of o's L) div 1 then
				set phase to item (item 2 of item i of o's L) of movingStates
				exit repeat
			else if curJDe = (item 1 of item i of o's L) div 1 then
				set phase to item (item 2 of item i of o's L) of exactStates
				exit repeat
			end if
		end repeat
		if phase = "" then
			set phase to item (((item 2 of item (count o's L) of o's L) + 1) mod 4) of movingStates
		end if
		local msg
		#		set msg to formatDate(YR, mnth, dy) & return & return & "The current Phase is " & phase & "."
		set msg to formatDate(YR, mnth, dy) & return & return & phase & "."
		
		tell application (path to frontmost application as text)
			display dialog msg with title "Todays Moon Phase:" buttons {"Quit", "Make Table"} cancel button 1 default button 2 with icon note
		end tell
	end moonPhaseOfToday
	
	
	to sort_items from alist
		-- Kai
		script o
			property L : alist
		end script
		repeat with i from 2 to count o's L
			set V to o's L's item i
			repeat with i from (i - 1) to 1 by -1
				tell o's L's item i to if item 1 of V < item 1 of it then
					set o's L's item (i + 1) to it
				else
					set o's L's item (i + 1) to V
					exit repeat
				end if
			end repeat
			if i is 1 and item 1 of V < item 1 of o's L's beginning then set o's L's item 1 to V
		end repeat
	end sort_items
	
	
	on moonphasesForMonth(YR, mnth)
		script o
			property n : {}
		end script
		
		local tyr, tmnth
		# checks for previous last quarter		
		if mnth = 1 then
			set {tyr, tmnth} to {YR - 1, 11.5}
		else
			set {tyr, tmnth} to {YR, mnth - 1.5}
		end if
		
		local test, JD_adjusted, dyr, mt, dy, L
		
		set {test, L} to {newMoon(tyr, tmnth), {}}
		set JD_adjusted to adjustJDNforLocalTZ(test)
		set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
		if dyr = YR and mt = mnth then set {end of L, end of o's n} to {{JD_adjusted, aNewMoon}, JD_adjusted}
		
		set test to firstQuarter(tyr, tmnth)
		set JD_adjusted to adjustJDNforLocalTZ(test)
		set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
		if dyr = YR and mt = mnth then set {end of L, end of o's n} to {{JD_adjusted, tFirstQuarter}, JD_adjusted}
		
		
		set test to fullMoon(tyr, tmnth)
		set JD_adjusted to adjustJDNforLocalTZ(test)
		set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
		if dyr = YR and mt = mnth then set {end of L, end of o's n} to {{JD_adjusted, aFullMoon}, JD_adjusted}
		
		set test to lastQuarter(tyr, tmnth)
		set JD_adjusted to adjustJDNforLocalTZ(test)
		set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
		if dyr = YR and mt = mnth then set {end of L, end of o's n} to {{JD_adjusted, tThirdQuarter}, JD_adjusted}
		
		# and if we have an extra new moon!
		set test to newMoon(YR, mnth)
		set JD_adjusted to adjustJDNforLocalTZ(test)
		set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
		if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set {end of L, end of o's n} to {{JD_adjusted, aNewMoon}, JD_adjusted}
		
		if length of L < 5 then
			# There can't be more than five moonphases of a month, so the very moment we have five
			# - We bail out!
			set test to firstQuarter(YR, mnth)
			set JD_adjusted to adjustJDNforLocalTZ(test)
			set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
			if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set {end of L, end of o's n} to {{JD_adjusted, tFirstQuarter}, JD_adjusted}
			
			if length of L < 5 then
				set test to fullMoon(YR, mnth)
				set JD_adjusted to adjustJDNforLocalTZ(test)
				set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
				if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set {end of L, end of o's n} to {{JD_adjusted, aFullMoon}, JD_adjusted}
				
				if length of L < 5 then
					set test to lastQuarter(YR, mnth)
					set JD_adjusted to adjustJDNforLocalTZ(test)
					set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
					if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set {end of L, end of o's n} to {{JD_adjusted, tThirdQuarter}, JD_adjusted}
					
					if length of L < 5 then
						
						# Before someone starts to think that I should have coalesced
						# the four nearly identical handlers and iterated over the handler here
						# instead of having it spread out boringly one by oneinto one, I'd like
						#  to retort that making code shorter, by adding if tests, and loops only 
						# makes it easier to maintain, not faster. It will also add to overall 
						# complexity here.
						
						# the phases that should be there, but not necessarily are!
						set JD_adjusted to findNewMoonForMonth(YR, mnth)
						if JD_adjusted ≠ 0 then
							set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
							if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set {end of L, end of o's n} to {{JD_adjusted, aNewMoon}, JD_adjusted}
						end if
						if length of L < 5 then
							set JD_adjusted to findFirstQuarterForMonth(YR, mnth)
							if JD_adjusted ≠ 0 then
								set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
								if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set {end of L, end of o's n} to {{JD_adjusted, tFirstQuarter}, JD_adjusted}
							end if
							if length of L < 5 then
								set JD_adjusted to findFullMoonForMonth(YR, mnth)
								if JD_adjusted ≠ 0 then
									set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
									if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set {end of L, end of o's n} to {{JD_adjusted, aFullMoon}, JD_adjusted}
								end if
								if length of L < 5 then
									set JD_adjusted to findlastQuarterForMonth(YR, mnth)
									if JD_adjusted ≠ 0 then
										set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
										if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set {end of L, end of o's n} to {{JD_adjusted, tThirdQuarter}, JD_adjusted}
									end if
									
									# additional phases that eventually may be there
									if length of L < 5 then
										if mnth = 12 then
											set {tyr, tmnth} to {YR + 1, 11}
										else
											set {tyr, tmnth} to {YR, mnth + 1}
										end if
										
										set test to newMoon(tyr, tmnth)
										set JD_adjusted to adjustJDNforLocalTZ(test)
										set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
										if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set end of L to {JD_adjusted, aNewMoon}
										
										if length of L < 5 then
											set test to firstQuarter(tyr, tmnth)
											set JD_adjusted to adjustJDNforLocalTZ(test)
											set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
											if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set end of L to {JD_adjusted, tFirstQuarter}
											
											if length of L < 5 then
												set test to fullMoon(tyr, tmnth)
												set JD_adjusted to adjustJDNforLocalTZ(test)
												set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
												if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set end of L to {JD_adjusted, aFullMoon}
												
												if length of L < 5 then
													set test to lastQuarter(tyr, tmnth)
													set JD_adjusted to adjustJDNforLocalTZ(test)
													set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
													if dyr = YR and mt = mnth and JD_adjusted is not in o's n then set end of L to {JD_adjusted, tThirdQuarter}
												end if
											end if
										end if
									end if
								end if
							end if
						end if
					end if
				end if
			end if
		end if
		
		sort_items from L
		-- Kai Edwards!
		return L
	end moonphasesForMonth
	
	
	
	on findNewMoonForMonth(YR, mnth)
		local JD, origMonth, dyr, mt, dy, forwardDir
		set {origMonth, forwardDir} to {mnth, false}
		set mnth to mnth - 0.75 # probability!!
		repeat while true
			set JD to newMoon(YR, mnth)
			set JD_adjusted to adjustJDNforLocalTZ(JD)
			
			set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
			if mt = origMonth then
				exit repeat
			else if (mnth div 1 + 1) = 1 and mt = 12 then
				# the month we got returned was lesser than the
				# the one we looked for.
				set mnth to 1.25
			else if (mnth div 1 + 1) = 12 and mt = 1 then
				# the month we got returned was greater then the
				# one we looked for.
				set mnth to 10.75
			else if mt = 12 and origMonth = 1 then
				# the month we got returned was greater than the
				# one we looked for.
				set forwardDir to true
				set mnth to mnth + 0.125
			else if mt = 1 and origMonth = 12 then
				# the month we got returned was greater 
				set mnth to mnth - 0.25
			else if origMonth < 12 and mt < origMonth then
				# the month we got returned was lesser than 
				# we looked for
				set forwardDir to true
				set mnth to mnth + 0.125
			else
				# if we have been forward, and then thrown backwards
				# then the cycle simply doesn't exist for this month!
				if forwardDir then return 0
				# the month was greater than the one we looked for.
				set mnth to mnth - 0.25
			end if
		end repeat
		return JD_adjusted
	end findNewMoonForMonth
	
	
	on newMoon(YR, mnt)
		local k, t, JD, M_, MPrime, F_, Omega, E_
		set k to approxK(YR + (mnt / 12)) div 1 # new moon 
		set t to centurieJuliae(k)
		set JD to meanPhaseOfMoon(t, k)
		set JD to computeAdditionalCorrectionsAllPhases(JD, t, k)
		set M_ to sunsMeanAnamoly(t, k)
		set MPrime to moonsMeanAnamoly(t, k)
		set F_ to moonsArgOfLatitude(t, k)
		set Omega to longOfAscendingNodeOfLunarOrbit(t, k)
		set E_ to eccentricityOfEarthsOrbit(t)
		
		return JD + (-0.4072 * (sin MPrime) + ¬
			0.17241 * E_ * (sin M_) + ¬
			0.01608 * (sin (2 * MPrime)) + ¬
			0.01039 * (sin (2 * F_)) + ¬
			0.00739 * E_ * (sin (MPrime - M_)) + ¬
			-0.00514 * E_ * (sin (MPrime + M_)) + ¬
			0.00208 * E_ * E_ * (sin (2 * M_)) + ¬
			-0.00111 * (sin (MPrime - 2 * F_)) + ¬
			-5.7E-4 * (sin (MPrime + 2 * F_)) + ¬
			5.6E-4 * E_ * (sin (2 * MPrime + M_)) + ¬
			-4.2E-4 * (sin (3 * MPrime)) + ¬
			4.2E-4 * E_ * (sin (M_ + 2 * F_)) + ¬
			3.8E-4 * E_ * (sin (M_ - 2 * F_)) + ¬
			-2.4E-4 * E_ * (sin (2 * MPrime - M_)) + ¬
			-1.7E-4 * (sin Omega) + ¬
			-7.0E-5 * (sin (MPrime + 2 * M_)) + ¬
			4.0E-5 * (sin (2 * MPrime - 2 * F_)) + ¬
			4.0E-5 * (sin (3 * M_)) + ¬
			3.0E-5 * (sin (MPrime + M_ - 2 * F_)) + ¬
			3.0E-5 * (sin (2 * MPrime + 2 * F_)) + ¬
			-3.0E-5 * (sin (MPrime + M_ + 2 * F_)) + ¬
			3.0E-5 * (sin (MPrime - M_ + 2 * F_)) + ¬
			-2.0E-5 * (sin (MPrime - M_ - 2 * F_)) + ¬
			-2.0E-5 * (sin (3 * MPrime + M_)) + ¬
			2.0E-5 * (sin (4 * MPrime)))
		
	end newMoon
	
	# Before someone starts to think that I should have coalesced
	# the four nearly identical handlers into one I'd like to retort
	# That making code shorter, by adding if tests, only makes it 
	# easier to maintain, not faster.
	
	on findFirstQuarterForMonth(YR, mnth)
		local JD, origMonth, dyr, mt, dy, forwardDir
		set {origMonth, forwardDir} to {mnth, false}
		set mnth to mnth - 0.75
		repeat while true
			set JD to firstQuarter(YR, mnth)
			set JD_adjusted to adjustJDNforLocalTZ(JD)
			set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
			if mt = origMonth then
				exit repeat
			else if (mnth div 1 + 1) = 1 and mt = 12 then
				# the month we got returned was lesser than the
				# the one we looked for.
				set mnth to 1.25
			else if (mnth div 1 + 1) = 12 and mt = 1 then
				# the month we got returned was greater then the
				# one we looked for.
				set mnth to 10.75
			else if mt = 12 and origMonth = 1 then
				# the month we got returned was greater than the
				# one we looked for.
				set forwardDir to true
				set mnth to mnth + 0.125
			else if mt = 1 and origMonth = 12 then
				# the month we got returned was greater 
				set mnth to mnth - 0.25
			else if origMonth < 12 and mt < origMonth then
				# the month we got returned was lesser than 
				# we looked for
				set forwardDir to true
				set mnth to mnth + 0.125
			else
				# if we have been forward, and then thrown backwards
				# then the cycle simply doesn't exist for this month!
				if forwardDir then return 0
				# the month was greater than the one we looked for.
				set mnth to mnth - 0.25
			end if
		end repeat
		return JD_adjusted
	end findFirstQuarterForMonth
	
	on firstQuarter(YR, mnt)
		local JD, W_
		set {JD, W_} to commonQuarters(YR, mnt, 0.25)
		return JD + W_
	end firstQuarter
	
	on commonQuarters(YR, mnt, factor)
		local k, t, JD, M_, MPrime, F_, Omega, E_, W_
		set k to approxK(YR + (mnt / 12)) div 1 + factor # Last Quater
		set t to centurieJuliae(k)
		set JD to meanPhaseOfMoon(t, k)
		set JD to computeAdditionalCorrectionsAllPhases(JD, t, k)
		set M_ to sunsMeanAnamoly(t, k)
		set MPrime to moonsMeanAnamoly(t, k)
		set F_ to moonsArgOfLatitude(t, k)
		set Omega to longOfAscendingNodeOfLunarOrbit(t, k)
		set E_ to eccentricityOfEarthsOrbit(t)
		set W_ to 0.00306 - 3.8E-4 * E_ * (cos M_) + 2.6E-4 * (cos MPrime) - 2.0E-5 * (cos (MPrime - M_)) + 2.0E-5 * (cos (MPrime + M_)) + 2.0E-5 * (cos (2 * F_))
		
		return {JD - 0.62801 * (sin MPrime) + ¬
			0.17172 * E_ * (sin M_) + ¬
			-0.01183 * E_ * (sin (MPrime + M_)) + ¬
			0.00862 * (sin (2 * MPrime)) + ¬
			0.00804 * (sin (2 * F_)) + ¬
			0.00454 * E_ * (sin (MPrime - M_)) + ¬
			0.00204 * E_ * E_ * (sin (2 * M_)) + ¬
			-0.0018 * (sin (MPrime - (2 * F_))) + ¬
			-7.0E-4 * (sin (MPrime + (2 * F_))) + ¬
			-4.0E-4 * (sin (3 * MPrime)) + ¬
			-3.4E-4 * E_ * (sin (2 * MPrime - M_)) + ¬
			3.2E-4 * E_ * (sin (M_ + (2 * F_))) + ¬
			3.2E-4 * E_ * (sin (M_ - (2 * F_))) + ¬
			-2.8E-4 * E_ * E_ * (sin (MPrime + (2 * M_))) + ¬
			2.7E-4 * E_ * (sin ((2 * MPrime) + M_)) + ¬
			-1.7E-4 * (sin Omega) + ¬
			-5.0E-5 * (sin (MPrime - M_ - (2 * F_))) + ¬
			4.0E-5 * (sin ((2 * MPrime) + (2 * F_))) + ¬
			-4.0E-5 * (sin (MPrime + M_ + (2 * F_))) + ¬
			4.0E-5 * (sin (MPrime - (2 * M_))) + ¬
			3.0E-5 * (sin (MPrime + M_ - (2 * F_))) + ¬
			3.0E-5 * (sin (3 * M_)) + ¬
			2.0E-5 * (sin ((2 * MPrime) - (2 * F_))) + ¬
			2.0E-5 * (sin (MPrime - M_ + (2 * F_))) + ¬
			-2.0E-5 * (sin ((3 * MPrime) + M_)), W_}
	end commonQuarters
	
	on findFullMoonForMonth(YR, mnth)
		local JD, origMonth, dyr, mt, dy, forwardDir
		set {origMonth, forwardDir} to {mnth, false}
		set mnth to mnth - 0.75
		repeat while true
			set JD to fullMoon(YR, mnth)
			set JD_adjusted to adjustJDNforLocalTZ(JD)
			set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
			if mt = origMonth then
				exit repeat
			else if (mnth div 1 + 1) = 1 and mt = 12 then
				# the month we got returned was lesser than the
				# the one we looked for.
				set mnth to 1.25
			else if (mnth div 1 + 1) = 12 and mt = 1 then
				# the month we got returned was greater then the
				# one we looked for.
				set mnth to 10.75
			else if mt = 12 and origMonth = 1 then
				# the month we got returned was greater than the
				# one we looked for.
				set forwardDir to true
				set mnth to mnth + 0.125
			else if mt = 1 and origMonth = 12 then
				# the month we got returned was greater 
				set mnth to mnth - 0.25
			else if origMonth < 12 and mt < origMonth then
				# the month we got returned was lesser than 
				# we looked for
				set forwardDir to true
				set mnth to mnth + 0.125
			else
				# if we have been forward, and then thrown backwards
				# then the cycle simply doesn't exist for this month!
				if forwardDir then return 0
				# the month was greater than the one we looked for.
				set mnth to mnth - 0.25
			end if
		end repeat
		return JD_adjusted
	end findFullMoonForMonth
	
	on fullMoon(YR, mnt)
		local k, t, JD, M_, MPrime, F_, Omega, E_
		set k to approxK(YR + (mnt / 12)) div 1 + 0.5 # full  Moon 
		set t to centurieJuliae(k)
		set JD to meanPhaseOfMoon(t, k)
		set JD to computeAdditionalCorrectionsAllPhases(JD, t, k)
		set M_ to sunsMeanAnamoly(t, k)
		set MPrime to moonsMeanAnamoly(t, k)
		set F_ to moonsArgOfLatitude(t, k)
		set Omega to longOfAscendingNodeOfLunarOrbit(t, k)
		set E_ to eccentricityOfEarthsOrbit(t)
		
		return JD - 0.40614 * (sin MPrime) + ¬
			0.17302 * E_ * (sin M_) + ¬
			0.01614 * (sin (2 * MPrime)) + ¬
			0.01043 * (sin (2 * F_)) + ¬
			0.00734 * E_ * (sin (MPrime - M_)) + ¬
			-0.00514 * E_ * (sin (MPrime + M_)) + ¬
			0.00209 * E_ * E_ * (sin (2 * M_)) + ¬
			-0.00111 * (sin (MPrime - 2 * F_)) + ¬
			-5.7E-4 * (sin (MPrime + 2 * F_)) + ¬
			5.6E-4 * E_ * (sin (2 * MPrime + M_)) + ¬
			-4.2E-4 * (sin (3 * MPrime)) + ¬
			4.2E-4 * E_ * (sin (M_ + 2 * F_)) + ¬
			3.8E-4 * E_ * (sin (M_ - 2 * F_)) + ¬
			-2.4E-4 * E_ * (sin (2 * MPrime - M_)) + ¬
			-1.7E-4 * (sin Omega) + ¬
			-7.0E-5 * (sin (MPrime + 2 * M_)) + ¬
			4.0E-5 * (sin (2 * MPrime - 2 * F_)) + ¬
			4.0E-5 * (sin (3 * M_)) + ¬
			3.0E-5 * (sin (MPrime + M_ - 2 * F_)) + ¬
			3.0E-5 * (sin (2 * MPrime + 2 * F_)) + ¬
			-3.0E-5 * (sin (MPrime + M_ + 2 * F_)) + ¬
			3.0E-5 * (sin (MPrime - M_ + 2 * F_)) + ¬
			-2.0E-5 * (sin (MPrime - M_ - 2 * F_)) + ¬
			-2.0E-5 * (sin (3 * MPrime + M_)) + ¬
			2.0E-5 * (sin (4 * MPrime))
	end fullMoon
	
	on findlastQuarterForMonth(YR, mnth)
		local JD, origMonth, dyr, mt, dy, forwardDir
		set {origMonth, forwardDir} to {mnth, false}
		set mnth to mnth - 0.75
		repeat while true
			set JD to lastQuarter(YR, mnth)
			set JD_adjusted to adjustJDNforLocalTZ(JD)
			set {dyr, mt, dy} to gregorianCalendarTriplet from JD_adjusted
			if mt = origMonth then
				exit repeat
			else if (mnth div 1 + 1) = 1 and mt = 12 then
				# the month we got returned was lesser than the
				# the one we looked for.
				set mnth to 1.25
			else if (mnth div 1 + 1) = 12 and mt = 1 then
				# the month we got returned was greater then the
				# one we looked for.
				set mnth to 10.75
			else if mt = 12 and origMonth = 1 then
				# the month we got returned was greater than the
				# one we looked for.
				set forwardDir to true
				set mnth to mnth + 0.125
			else if mt = 1 and origMonth = 12 then
				# the month we got returned was greater 
				set mnth to mnth - 0.25
			else if origMonth < 12 and mt < origMonth then
				# the month we got returned was lesser than 
				# we looked for
				set forwardDir to true
				set mnth to mnth + 0.125
			else
				# if we have been forward, and then thrown backwards
				if forwardDir then return 0
				set mnth to mnth - 0.25
			end if
			
		end repeat
		return JD_adjusted
	end findlastQuarterForMonth
	
	on lastQuarter(YR, mnt)
		local JD, W_
		set {JD, W_} to commonQuarters(YR, mnt, 0.75)
		return JD - W_
	end lastQuarter
	
	on computeAdditionalCorrectionsAllPhases(JDE, t, k)
		return (JDE + 3.25E-4 * (sin (map360(299.77 + 0.107408 * k - 0.009173 * t * t) * degrad)) + ¬
			1.65E-4 * (sin (map360(251.88 + 0.016321 * k) * degrad)) + ¬
			1.64E-4 * (sin (map360(251.83 + 26.651886 * k) * degrad)) + ¬
			1.26E-4 * (sin (map360(349.42 + 36.412478 * k) * degrad)) + ¬
			1.1E-4 * (sin (map360(84.66 + 18.206239 * k) * degrad)) + ¬
			6.2E-5 * (sin (map360(141.74 + 53.303771 * k) * degrad)) + ¬
			6.0E-5 * (sin (map360(207.14 + 2.453732 * k) * degrad)) + ¬
			5.6E-5 * (sin (map360(154.84 + 7.30686 * k) * degrad)) + ¬
			4.7E-5 * (sin (map360(34.52 + 27.261239 * k) * degrad)) + ¬
			4.2E-5 * (sin (map360(207.19 + 0.121824 * k) * degrad)) + ¬
			4.0E-5 * (sin (map360(291.34 + 1.844379 * k) * degrad)) + ¬
			3.7E-5 * (sin (map360(161.72 + 24.198154 * k) * degrad)) + ¬
			3.5E-5 * (sin (map360(239.56 + 25.513099 * k) * degrad)) + ¬
			2.3E-5 * (sin (map360(331.55 + 3.592518 * k) * degrad)))
	end computeAdditionalCorrectionsAllPhases
	
	
	on eccentricityOfEarthsOrbit(t)
		# (Meeus 45.6)
		return (1 + (t * (-0.002516 + (t * -7.4E-6))))
	end eccentricityOfEarthsOrbit
	
	on longOfAscendingNodeOfLunarOrbit(k, t)
		# Meeus ( 47.7) 
		return (map360(124.7746 - (k * 1.5637558) + (t * t * (0.0020691 + (t * 2.15E-6)))) * degrad)
	end longOfAscendingNodeOfLunarOrbit
	
	on moonsArgOfLatitude(t, k)
		# meeus ( 47.6 )
		return (map360((160.7108 + (k * 390.67050274)) + (t * t * (-0.0016341 + (t * (-2.27E-6 + (t * 1.1E-8)))))) * degrad)
	end moonsArgOfLatitude
	
	on moonsMeanAnamoly(t, k)
		# Meeus (47.5)
		return (map360(201.5643 + (k * 385.81693528) + (t * t * (0.0107438 + (t * 1.239E-5 - (t * 5.8E-8))))) * degrad)
	end moonsMeanAnamoly
	
	on sunsMeanAnamoly(t, k)
		# Meeus (47.4)
		return (map360(2.5534 + (k * 29.10535669) + (t * t * (-2.18E-5 + (t * -1.1E-7)))) * degrad)
	end sunsMeanAnamoly
	
	on meanPhaseOfMoon(t, k)
		# Meeus (47.1)
		return 2.45155009765E+6 + (29.530588853 * k) + (t * (t * (1.337E-4 + (t * (-1.5E-7 + (t * 7.3E-10))))))
	end meanPhaseOfMoon
	
	on centurieJuliae(k)
		# Meeus (47.3)
		return (k / 1236.85)
	end centurieJuliae
	
	on approxK(decyear)
		# Meeus (47.2)
		set a to ((decyear - 2000) * 12.3685)
		return round ((decyear - 2000) * 12.3685) rounding to nearest
	end approxK
	
	# ----
	
	on map360(aDecNumber)
		set aDecNumber to aDecNumber mod 360
		if aDecNumber < 0 then set aDecNumber to aDecNumber + 360
		return aDecNumber
	end map360
	
	# end script
	
	on GregorianDate for theJDN
		-- Added 2013.10.7
		-- Nigel Garvey http://macscripter.net/viewtopic.php?pid=167391#p167391 post # 6
		return JD0 + theJDN * days
	end GregorianDate
	
	on TZtoTZ(TZ1date, TZ1, TZ2)
		# Nigel Garvey
		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 JulianDayNumber for theDate
		-- Added 2013.10.7
		-- Nigel Garvey http://macscripter.net/viewtopic.php?pid=167391#p167391 post # 6
		return (theDate - JD0) div days
	end JulianDayNumber
	
	
	on timeTriplet from JDN
		# This returns as standard time triplet, where the first point of time of a day is 00:00
		# and the last point of a time is 23:59:59.9999999 and so on..
		local timeFraction, _hours, _minutes, _seconds
		set timeFraction to (JDN mod 1) * 24
		set _hours to (timeFraction div 1 as integer)
		set timeFraction to (timeFraction - _hours)
		set _minutes to timeFraction * 60 div 1
		set _seconds to (timeFraction - (_minutes / 60)) * 3600
		return {_hours, _minutes, _seconds}
	end timeTriplet
	
	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
	
	on civilTimeToGmtFromLocalTZInSeconds(anASDate)
		# Nigel Garvey and Adam Bell made this possible!
		return (anASDate - TZtoTZ(anASDate, (do shell script "readlink /etc/localtime |sed -n 's_/[^/]*/[^/]*/[^/]*/\\([^/]*\\)/\\([^/]*\\).*_\\1/\\2_p'"), "GMT"))
	end civilTimeToGmtFromLocalTZInSeconds
	
	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 adjustJDNforLocalTZ(JDE)
		# it adds the time difference to GMT, and adds the 12 hour, to askew the day to the civil
		# time system, this may also increment the date, which is perfectly fine, since a Julian day
		# starts 12 hours before a civil day.
		return (JDE + ((civilTimeToGmtFromLocalTZInSeconds(GregorianDate for (JDE + 0.5))) / 86400) + 0.5)
	end adjustJDNforLocalTZ
	
	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
	
	on setClockFormat(m)
		local clock_format
		considering case
			tell (current date)'s time string
				if it contains "AM" or it contains "PM" then
					set clock_format to m's big_AMPM
				else if it contains "Am" or it contains "Pm" then
					set clock_format to m's ususal_AmPm
				else if it contains "am" or it contains "pm" then
					set clock_format to m's small_ampm
				else
					set clock_format to m's hour24
				end if
			end tell
		end considering
		return clock_format
	end setClockFormat
	
	on formatClockTime(_hr, _min, _sec)
		script o
			property big_AMPM : 3
			property ususal_AmPm : 2
			property small_ampm : 1
			property hour24 : 0
		end script
		global am_pm_clock
		local a
		try
			set a to am_pm_clock
		on error
			set am_pm_clock to setClockFormat(o)
		end try
		
		if am_pm_clock = o's hour24 then
			return (pad((_hr div 1)) & ":" & pad(_min) & ":" & pad((_sec div 1)))
		else if am_pm_clock = o's big_AMPM then
			if (_hr div 12) > 0 then
				return (padHourWithBlank((_hr div 12)) & ":" & pad(_min) & ":" & pad((_sec div 1)) & " PM")
			else
				return (padHourWithBlank((_hr div 1)) & ":" & pad(_min) & ":" & pad((_sec div 1)) & " AM")
			end if
		else if am_pm_clock = o's ususal_AmPm then
			if (_hr div 12) > 0 then
				return (padHourWithBlank((_hr div 12)) & ":" & pad(_min) & ":" & pad((_sec div 1)) & " Pm")
			else
				return (padHourWithBlank((_hr div 1)) & ":" & pad(_min) & ":" & pad((_sec div 1)) & " Am")
			end if
		else #  am_pm_clock = small_ampm  
			if (_hr div 12) > 0 then
				return (padHourWithBlank((_hr div 12)) & ":" & pad(_min) & ":" & pad((_sec div 1)) & " pm")
			else
				return (padHourWithBlank((_hr div 1)) & ":" & pad(_min) & ":" & pad((_sec div 1)) & " am")
			end if
		end if
	end formatClockTime
	
	on padHourWithBlank(n)
		-- "Special Handler" for padding the 12 hour clock with blanks.
		return (text -2 thru -1 of (" " & n as text))
	end padHourWithBlank
	
	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 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
	
	-- Ask the user for the range of dates to be covered.
	-- from a start date, up to, and including enddate.
	on getRangeOfYears(initialYear, initialRange, scriptTitle)
		-- Added 2013.10.7, used in Equation Of Time, high accuracy..
		-- 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.
		local anicon, theMessage, StartYear, EndYear
		set anicon to note
		set theMessage to "Enter a range of years or a single year."
		set StartYear to initialYear
		repeat
			if initialRange = "Century" then
				set EndYear to StartYear + 100
			else if initialRange = "SemiCentury" then
				set EndYear to StartYear + 50
			else if initialRange = "Decade" then
				set EndYear to StartYear + 10
			else if initialRange = "Year" then
				set EndYear to StartYear + 5
				# We'll add one day for a leap year
			end if
			local d1, d2, yearRange
			set d1 to StartYear as text
			set d2 to EndYear as text
			
			
			tell application (path to frontmost application as text)
				set yearRange 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
			local tids, gotOne, yearRangeStart, yearRangeEnd
			try
				set {tids, AppleScript's text item delimiters, gotOne} to {AppleScript's text item delimiters, "-", false}
				set {yearRangeStart, gotOne} to {text item 1 of yearRange as number, true}
				
				set yearRangeEnd to text item 2 of yearRange as number
				set AppleScript's text item delimiters to tids
				exit repeat
			on error
				if gotOne and length of text items of yearRange = 1 then
					set yearRangeEnd to text item 1 of yearRange as number
					set AppleScript's text item delimiters to tids
					exit repeat
				else if gotOne then
					set today to yearRangeStart
				else
					set today to initialYear
				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
		return {yearRangeStart, yearRangeEnd}
	end getRangeOfYears
end script

tell MoonPhases to run


Hello.

I added a property for setting a 12 hour am/pm clock, for those of you that are accustomed to that particular time format. I currently don’t know how to easily deduct that info, so you’ll have to set the property to true in order to turn the feature on!

I also neglected to say in the original post that the times of the lunar phases are returned according to your Mac’s date and time settings, in the preferences.

Hello.

I added functionality to automatically detect clock format: 24 hour, AM/PM, Am/Pm, and am/pm (for prosperity).

Hello.

I changed the detection of clock format that is am/pm or a 24 hour clock to be a private subsystem, not relying on global properties, as I deemed that the soundest way to implement it. Not that it is much better if the user just shares a copy of his script, as the previous clock format will still be in effect until the script is edited, but the use of this design principle makes it easier to implement the functionality (less editing), and steal. :slight_smile:

The script is still in post # 1, and you still need Satimage.osax.

Hello.

I just added forced initialization of the variables that controls the format of clock, and date according to locale.

I do that with an undef() handler that just has a return statement in it, effectively undefining the variable that is assigned the result.

It is a funny quirk here: It works if you do it in a single statement, but at least Script Debugger barks, when I try to do it simultaneously to several variables at once with a list statement. :slight_smile:

Hello.

Here is another script, that shows moonphases for a month at a time. Which should compute a little bit faster than the script above, you really need to have Satimage.osax installed, or provide your handlers for trigs and so on.

# © 2013 McUsr. You may not post this anywhere else. http://macscripter.net/viewtopic.php?pid=180859#p180859
# Revised after advice by Nigel Garvey,
# New Moon January the 6  year 2000
# http://kalender-365.de/manekalender.php?yy=2000
# http://www.calendar-365.com/moon/moon-phases.html
# http://en.wikipedia.org/wiki/Lunar_phase
# ----
# http://en.wikipedia.org/wiki/Full_moon_cycle
# http://en.wikipedia.org/wiki/New_moon
# http://en.wikipedia.org/wiki/Full_moon

property fullMoonsPerFullmoonCycle : 14
property daysPerFullmoonCycle : 411.78443
# what kind of moon cycle to follow, a synodic month contains the
# phases cf. link.
property synodicMonth : 29.530588853
property EpochJD2000 : 2.451545E+6

global w, exc
set TGMT to ((time to GMT) / hours) / 24
# log "TGMT : " & TGMT

# ENTER YEAR AND MONTH FOR THE MOON HERE
set sYear to 2015
set sMonth to 5
set m to fullMoonNumberForCalendricMonth(sYear, sMonth, 1)

set JD to JDN(sYear, sMonth, 1, 0, 0, 0)
# log "JD = " & JD

set JDE to calculateJDE(JD, sYear, sMonth)
# log "JDE = " & JDE

set t to (T_inJulianCenturies for JDE)
# log "T = " & T

set MAS to (MeanAnomalyForSunInRad for t)
# log "MAS = " & MAS
set MAM to (MeanAnomalyForMoonInRad for t)
# log "MAM = " & MAM
set exc to (MoonsAmplitude for t)
# log "E = " & exc
# Correction of quarters
set w to 0.0036 - (3.8E-4 * (MoonsAmplitude for t) * (cos MAS)) + (2.6E-4 * (cos MAM)) - (2.0E-5 * (cos (MAM - MAS))) + (2.0E-5 * (cos (MAS + MAM))) + (2.0E-5 * (cos (2 * (MoonsArgOfLatitude for t))))
# log " W = " & w
# NEW MOON
#log "moon number = " & m
repeat while true
	set dnr to meanSyzygyForMoonNr(m, 5.597661) # New
	#	log "mean syzygy = " & dnr
	#  Calculations in order to find the true syzygy
	
	set tdnr to (trueSyzygyForNewMoonNr(dnr, MAM, MAS) + (EpochJD2000) - 1)
	set tdnr to tdnr - ((deltaT for (decimalYear(sYear, sMonth))) / 86400) + 0.5 - TGMT
	#	log "true syzygy = " & tdnr
	#	log "tdnr - Jd = " & (tdnr - JD)
	if tdnr - JD > 30 then
		set m to m - 1
	else if JD - tdnr > 30 then
		set m to m + 1
	else
		exit repeat
	end if
end repeat
set {yr, MT, d} to (gregorianCalendarTriplet from tdnr)
set MTEXT to "Calendar date for New Moon = Yr " & yr & " mt = " & MT & " d = " & d & return

# FIRST QUARTER

set dnr to meanSyzygyForMoonNr(m, (12.980308 + (0.5 / exc))) # FirstQuarter
# log "mean syzygy = " & dnr
#  Calculations in order to find the true syzygy

set tdnr to (trueSyzygyForFirstQuarterMoonNr(dnr, MAM, MAS) + EpochJD2000 - 1)
set tdnr to tdnr - ((deltaT for (decimalYear(sYear, sMonth))) / 86400) # + 0.5 - TGMT
# log "true syzygy = " & tdnr
set {yr, MT, d} to (gregorianCalendarTriplet from tdnr)
set MTEXT to MTEXT & "Calendar date for First Quarter = Yr " & yr & " mt = " & MT & " d = " & d & return

# FULL MOON

set dnr to meanSyzygyForMoonNr(m, 20.362955) # Full
# log "mean syzygy = " & dnr
#  Calculations in order to find the true syzygy

set tdnr to (trueSyzygyForFullMoonNr(dnr, MAM, MAS) + EpochJD2000 - 1)
set tdnr to tdnr - ((deltaT for (decimalYear(sYear, sMonth))) / 86400) # + 0.5 - TGMT
# log "true syzygy = " & tdnr


set {yr, MT, d} to (gregorianCalendarTriplet from tdnr)
set MTEXT to MTEXT & "Calendar date for Full Moon = Yr " & yr & " mt = " & MT & " d = " & d & return


# LAST QUARTER 

set dnr to meanSyzygyForMoonNr(m, (27.745602 + (1 / exc))) # LastQuarter
# log "mean syzygy = " & dnr
#  Calculations in order to find the true syzygy

set tdnr to (trueSyzygyForLastQuarterMoonNr(dnr, MAM, MAS) + EpochJD2000 - 1)
set tdnr to tdnr - ((deltaT for (decimalYear(sYear, sMonth))) / 86400) # + 0.5 - TGMT
# log "true syzygy = " & tdnr

set {yr, MT, d} to (gregorianCalendarTriplet from tdnr)
set MTEXT to MTEXT & "Calendar date for Last Quarter = Yr " & yr & " mt = " & MT & " d = " & d & return

tell application "TextEdit"
	make new document
	set text of its document 1 to MTEXT
end tell


to MeanAnomalyForSunInRad for t
	# Meeus (45.3)
	return (((357.52772 + (3.599905034E+4 * t) + (1.603E-4 * (t ^ 2)) + ((t ^ 3) / 300000)) mod 360) * pi / 180)
end MeanAnomalyForSunInRad

to MeanAnomalyForMoonInRad for t
	# Meeus ( 45.4)
	return (((134.96298 + (4.77198867398E+5 * t) + (0.0086972 * (t ^ 2)) + ((t ^ 3) / 56250)) mod 360) * pi / 180)
end MeanAnomalyForMoonInRad

to MoonsArgOfLatitude for t
	# Meeus (45.5)
	return (((93.2720993 + (4.832020175273E+5 * t) - (0.0034029 * (t ^ 2)) - ((t ^ 3) / 3526000) + ((t ^ 4) / 8.6331E+8)) mod 360) * pi / 180)
end MoonsArgOfLatitude

to MoonsAmplitude for t
	# Meeus 45.6
	return (1 - (0.002516 * t) - (7.4E-7 * (t ^ 2)))
end MoonsAmplitude
to calculateJDE(JD, yr, mnt)
	# Meus p. 131, but the  deltaT is from Nasa
	
	return (JD + ((deltaT for (decimalYear(yr, mnt))) / 86400))
end calculateJDE

# returns the time in Julian Centuries
# JDE is a JD extrepolated to Dynamical time.

to T_inJulianCenturies for aJDE
	return ((aJDE - 2451545) / 36525)
end T_inJulianCenturies

to fullMoonNumberForCalendricMonth(y, m, d)
	# The Calendric Month is a Gregorian Month.
	
	# get todays julian day number smart to start with the first of the 
	# month, maybe to add the coeffecient of Terrestial Time too.
	set curJD to JDN(y, m, d, 0, 0, 0)
	set JDSinceEpoch to curJD - EpochJD2000
	#	log " Julian days since epoch: " & JDSinceEpoch
	set fullMoonCycles to JDSinceEpoch div daysPerFullmoonCycle
	#	log "fullMoonCycles since EpochJD2000: " & fullMoonCycles
	#		set daysLeftOver to JDSinceEpoch mod fullMoonCycles
	set daysLeftOver to JDSinceEpoch mod daysPerFullmoonCycle
	#	log "daysLeftOver: " & daysLeftOver
	set extraMoonPhases to daysLeftOver div synodicMonth
	#	log "extraMoonPhases: " & extraMoonPhases
	# calulates the moonNumber we are after:
	set moonNumber to (fullMoonCycles * fullMoonsPerFullmoonCycle) + extraMoonPhases
	# So, now we know the moon number we are looking for the full-moon for.
	# I kind of know the current moons "age" now, if it should matter later
	return moonNumber
	
end fullMoonNumberForCalendricMonth

# Simplified: thanks to Nigel Garvey.
to meanSyzygyForMoonNr(moonNumber, primeFactor)
	local d, corrD
	# http://en.wikipedia.org/wiki/Full_moon
	# Returns the mean syzygy, which must be corrected.
	set d to (primeFactor + (29.530588861 * moonNumber) + ((102.026 * 1.0E-11) * (moonNumber ^ 2)))
	#	log "d: " & d
	set corrD to d + (-7.39E-4 - (235 * 1.0E-11) * (moonNumber ^ 2))
	#	log "corrD: " & corrD
	return corrD
end meanSyzygyForMoonNr


to trueSyzygyForNewMoonNr(meansyzygy, MAM, MAS)
	# instant dynamical time, (21.1) Meeus
	# (45.6)
	return (meansyzygy + (-0.4072 * (sin MAM)) + (0.01608 * (sin (2 * MAM))) + (0.17241 * (sin MAS)))
end trueSyzygyForNewMoonNr

to trueSyzygyForFirstQuarterMoonNr(meansyzygy, MAM, MAS)
	# instant dynamical time, (21.1) Meeus
	# (45.6)
	return (meansyzygy + (-0.62801 * (sin MAM)) + (0.00862 * (sin (2 * MAM))) + (0.17172 * (sin MAS)) + w)
end trueSyzygyForFirstQuarterMoonNr

to trueSyzygyForFullMoonNr(meansyzygy, MAM, MAS)
	# instant dynamical time, (21.1) Meeus
	# (45.6)
	return (meansyzygy + (-0.40614 * (sin MAM)) + (0.01614 * (sin (2 * MAM))) + (0.17302 * (sin MAS)))
end trueSyzygyForFullMoonNr

to trueSyzygyForLastQuarterMoonNr(meansyzygy, MAM, MAS)
	# instant dynamical time, (21.1) Meeus
	# (45.6)
	return (meansyzygy + (-0.62801 * (sin MAM)) + (0.00862 * (sin (2 * MAM))) + (0.17172 * (sin MAS)) - w)
end trueSyzygyForLastQuarterMoonNr

on timeInJulianCenturies(JDE)
	# returns the time in Julian Centuries
	# JDE is a JD extrepolated to Dynamical time.
end timeInJulianCenturies

# aDecimalYear = year + (month - 0.5) / 12
to decimalYear(yr, mnth)
	return (yr + (mnth - 0.5) / 12)
end decimalYear

to deltaT for aDecimalYear
	# http://eclipse.gsfc.nasa.gov/SEcat5/deltatpoly.html
	# simplifed, halved the contents of the if tests, 
	# and speeded up the multiplication of powers thanks to Nigel Garvey
	local u
	if aDecimalYear < -500 then
		set u to (aDecimalYear - 1820) / 100
		return (-20 + 32 * (u ^ 2))
	else if aDecimalYear < 500 then
		set u to aDecimalYear / 100
		return (1.05836E+4 - (1014.41 * u) + (33.78311 * (u ^ 2)) - (5.952053 * (u ^ 3)) - (0.1798452 * (u ^ 4)) + (0.022174192 * (u ^ 5)) + (0.0090316521 * (u ^ 6)))
	else if aDecimalYear < 1600 then
		set u to (aDecimalYear - 1000) / 100
		return (1574.2 - (556.01 * u) + (71.23472 * (u ^ 2)) + (0.319781 * (u ^ 3)) - (0.8503463 * (u ^ 4)) - (0.005050998 * (u ^ 5)) + (0.0083572073 * (u ^ 6)))
	else if aDecimalYear < 1700 then
		set u to aDecimalYear - 1600
		return (120 - (0.9808 * u) - (0.01532 * (u ^ 2)) + ((u ^ 3) / 7129))
	else if aDecimalYear < 1800 then
		set u to aDecimalYear - 1700
		return (8.83 + (0.1603 * u) - (0.0059285 * (u ^ 2)) + (1.3336E-4 * (u ^ 3)) - ((u ^ 4) / 1174000))
	else if aDecimalYear < 1860 then
		set u to aDecimalYear - 1800
		return (13.72 - (0.332447 * u) + (0.0068612 * (u ^ 2)) + (0.0041116 * (u ^ 3)) - (3.7436E-4 * (u ^ 4)) + (1.21272E-5 * (u ^ 5)) - (1.69E-7 * (u ^ 6)) + (8.75E-10 * (u ^ 7)))
	else if aDecimalYear < 1900 then
		set u to aDecimalYear - 1860
		return (7.62 + (0.5737 * u) - (0.251574 * (u ^ 2)) + (0.01680668 * (u ^ 3)) - (4.473624E-4 * (u ^ 4)) + ((u ^ 5) / 233174))
	else if aDecimalYear < 1920 then
		set u to aDecimalYear - 1900
		return (-2.79 + (1.494119 * u) - (0.0598939 * (u ^ 2)) + (0.0061966 * (u ^ 3)) - (1.97E-4 * (u ^ 4)))
	else if aDecimalYear < 1941 then
		set u to aDecimalYear - 1920
		return (21.2 + (0.84493 * u) - (0.0761 * (u ^ 2)) + (0.0020936 * (u ^ 3)))
	else if aDecimalYear < 1961 then
		set u to aDecimalYear - 1950
		return (29.07 + (0.407 * u) - ((u ^ 2) / 233) + ((u ^ 3) / 2547))
	else if aDecimalYear < 1986 then
		set u to aDecimalYear - 1975
		return (45.45 + (1.067 * u) - ((u ^ 2) / 260) - ((u ^ 3) / 718))
	else if aDecimalYear < 2005 then
		set u to aDecimalYear - 2000
		return (63.86 + (0.3345 * u) - (0.060374 * (u ^ 2)) + (0.0017275 * (u ^ 3)) + (6.51814E-4 * (u ^ 4)) + (2.373599E-5 * (u ^ 5)))
	else if aDecimalYear < 2050 then
		set u to aDecimalYear - 2000 # (correct!)
		return (62.92 + (0.32217 * u) + (0.005589 * (u ^ 2)))
	else if aDecimalYear < 2150 then
		return (-20 + (32 * (((aDecimalYear - 1820) / 100) ^ 2 - (0.5628 * (2150 - aDecimalYear)))))
	else
		set u to (aDecimalYear - 1820) / 100
		return (-20 + (32 * (u ^ 2)))
	end if
end deltaT

to JDN(_year, _month, _day, _hour, _min, _sec)
	# 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

to 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 as integer
	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