Solar Equinoxes, Solstices and Quarters

Hello.

This script gives you a report over the First day of the seasons of a year, together with the dates of the Equnioxes and Solstices.

You’ll need Satimage.osax, in order to run this.

Enjoy.

script solarEvents
	property scriptTitle : "Solar Equinoxes, Solstices and Quarters"
	# Compute Quarters:
	property quarterNames : {"First day of Spring", "First Day of Summer", "First Day of Autumn", "First Day of Winter"}
	property equinoxNames : {"Vernal Equinox", "Summer Solstice", "Auturnal Equinox", "Winter Solstice"}
	property quarterAngles : {315.0, 45.0, 135.0, 225.0}
	property debug : true
	property tlvl : me
	property EpochJD2000 : 2.451545E+6
	property JulianDaysInCentury : 36525
	
	script solaris
		# Solar coordinates low accuracy Meeus chapter 24 p. 151
		
		on longitude(JDE_t)
			local T_, L0, m, ecc, equOfC, trueLongitude, trueAnamoly, R, omega, apparent, meanObliguityOfEcliptic
			set T_ to JulianCenturies of tlvl by JDE_t
			
			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)
			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))))
			
			# The suns right ascension, and declination, in equatorial coordinates
			
			set meanObliguityOfEcliptic to (meanObliquityOfEcliptic_Laskar for JDE_t)
			return mapWithinModulus((sunsRightAscension(meanObliguityOfEcliptic, apparentLong) * 180 / pi), 360)
		end longitude
		
		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
		
		on sunsDeclination(trueObliguity, longitude)
			# (24.7)
			return (asin ((sin (trueObliguity * pi / 180)) * (sin (longitude * pi / 180))))
		end sunsDeclination
		
		on meanObliquityOfEcliptic_Laskar for JD
			# Meuus (21.3) p.135
			# Extracted from sidereal Lib.
			local U_
			set U_ to (JD - 2451545) / 3652500
			return (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_))))))))))))))))))))
		end meanObliquityOfEcliptic_Laskar
		on mapWithinModulus(aDecNumber, aModulus)
			set aDecNumber to aDecNumber mod aModulus
			if aDecNumber < 0 then set aDecNumber to aDecNumber + aModulus
			return aDecNumber
		end mapWithinModulus
		
	end script
	
	script equinoxlib
		property s_table : {¬
			{485, 324.96, 1934.136}, ¬
			{203, 337.23, 3.2964467E+4}, ¬
			{199, 342.08, 20.186}, ¬
			{182, 27.85, 4.45267112E+5}, ¬
			{156, 73.14, 4.5036886E+4}, ¬
			{136, 171.52, 2.2518443E+4}, ¬
			{77, 222.54, 6.5928934E+4}, ¬
			{74, 296.72, 3034.906}, ¬
			{70, 243.58, 9037.513}, ¬
			{58, 119.81, 3.3718147E+4}, ¬
			{52, 297.17, 150.678}, ¬
			{50, 21.02, 2281.226}, ¬
			{45, 247.54, 2.9929562E+4}, ¬
			{44, 325.15, 3.1555956E+4}, ¬
			{29, 60.93, 4443.417}, ¬
			{18, 155.12, 6.7555328E+4}, ¬
			{17, 288.79, 4562.452}, ¬
			{16, 198.04, 6.2894029E+4}, ¬
			{14, 199.76, 3.1436921E+4}, ¬
			{12, 95.39, 1.4577848E+4}, ¬
			{12, 287.11, 3.1931756E+4}, ¬
			{12, 320.81, 3.4777259E+4}, ¬
			{9, 227.73, 1222.114}, ¬
			{8, 15.45, 1.6859074E+4}}
		
		on sumOfS(T)
			local sigma
			set sigma to 0
			repeat with i from 1 to (length of s_table)
				set sigma to sigma + ((item 1 of item i of s_table) * (cos (((item 2 of item i of s_table) * pi / 180) + ((item 3 of item i of s_table) * pi / 180 * T))))
			end repeat
			return sigma
		end sumOfS
		on B_year(YR)
			return ((YR - 2000) / 1000)
		end B_year
		
		on jdeforvernalequinox(YR)
			# Returns JDE, that must be corrected to mean time, and civiltime.
			# the correction for civiltime, is done by adding 0.5 to the hours.
			local JDE, T_, w, deltaLambda, s
			set JDE to ApproxVernalEquinox(YR)
			set T_ to (JulianCenturies by JDE) of tlvl
			set w to (3.5999373E+4 * T_ - 2.47) * pi / 180
			set deltaLambda to 1 + (0.00334 * (cos (w))) + 7.0E-4 * (cos (2 * w))
			set s to sumOfS(T_)
			return (JDE + (1.0E-5 * s) / deltaLambda)
		end jdeforvernalequinox
		
		on ApproxVernalEquinox(YR)
			local JDE0
			if YR ≤ 1000 then
				set YR to (YR / 1000)
				set JDE0 to 1.72113929189E+6 + (YR * (3.652421374E+5 + (YR * (0.06134 + (YR * (0.00111 + (YR * 7.1E-4)))))))
			else
				set YR to ((YR - 2000) / 1000)
				set JDE0 to 2.45162380984E+6 + (YR * (3.6524237404E+5 + (YR * (0.05169 + (YR * (-0.00411 + (YR * -5.7E-4)))))))
			end if
			return JDE0
		end ApproxVernalEquinox
		
		on JDEForSummerSolstice(YR)
			# Returns JDE, that must be corrected to mean time, and civiltime.
			# the correction for civiltime, is done by adding 0.5 to the hours.
			local JDE, T_, w, deltaLambda, s
			set JDE to ApproxSummerSolstice(YR)
			set T_ to (JulianCenturies by JDE) of tlvl
			set w to (3.5999373E+4 * T_ - 2.47) * pi / 180
			set deltaLambda to 1 + (0.00334 * (cos (w))) + 7.0E-4 * (cos (2 * w))
			set s to sumOfS(T_)
			return (JDE + (1.0E-5 * s) / deltaLambda)
		end JDEForSummerSolstice
		
		on ApproxSummerSolstice(YR)
			local JDE0
			if YR ≤ 1000 then
				set YR to (YR / 1000)
				set JDE0 to 1.72123325401E+6 + (YR * (3.6524172562E+5 + (YR * (-0.05323 + (YR * (0.00907 + (YR * 2.5E-4)))))))
			else
				set YR to ((YR - 2000) / 1000)
				set JDE0 to 2.45171656767E+6 + (YR * (3.6524162603E+5 + (YR * (0.00325 + (YR * (0.00888 + (YR * -3.0E-4)))))))
			end if
			return JDE0
		end ApproxSummerSolstice
		
		on JDEForAuturnalEquinox(YR)
			# Returns JDE, that must be corrected to mean time, and civiltime.
			# the correction for civiltime, is done by adding 0.5 to the hours.
			local JDE, T_, w, deltaLambda, s
			set JDE to ApproxAuturnalEquinox(YR)
			set T_ to (JulianCenturies by JDE) of tlvl
			set w to (3.5999373E+4 * T_ - 2.47) * pi / 180
			set deltaLambda to 1 + (0.00334 * (cos (w))) + 7.0E-4 * (cos (2 * w))
			set s to sumOfS(T_)
			return (JDE + (1.0E-5 * s) / deltaLambda)
		end JDEForAuturnalEquinox
		
		on ApproxAuturnalEquinox(YR)
			local JDE0
			if YR ≤ 1000 then
				set YR to (YR / 1000)
				set JDE0 to 1.72132570455E+6 + (YR * (3.6524249558E+5 + (YR * (-0.11677 + (YR * (-0.00297 + (YR * 7.4E-4)))))))
			else
				set YR to ((YR - 2000) / 1000)
				set JDE0 to 2.45181021715E+6 + (YR * (3.6524201767E+5 + (YR * (-0.11575 + (YR * (0.00337 + (YR * 7.8E-4)))))))
			end if
			return JDE0
		end ApproxAuturnalEquinox
		
		on jdeforwintersolstice(YR)
			# Returns JDE, that must be corrected to mean time, and civiltime.
			# the correction for civiltime, is done by adding 0.5 to the hours.
			local JDE, T_, w, deltaLambda, s
			set JDE to ApproxWinterSolstice(YR)
			set T_ to (JulianCenturies by JDE) of tlvl
			set w to (3.5999373E+4 * T_ - 2.47) * pi / 180
			set deltaLambda to 1 + (0.00334 * (cos (w))) + 7.0E-4 * (cos (2 * w))
			set s to sumOfS(T_)
			return (JDE + (1.0E-5 * s) / deltaLambda)
		end jdeforwintersolstice
		
		on ApproxWinterSolstice(YR)
			local JDE0
			if YR ≤ 1000 then
				set YR to (YR / 1000)
				set JDE0 to 1.72141439987E+6 + (YR * (3.6524288257E+5 + (YR * (-0.00769 + (YR * (-0.00933 - (YR * 6.0E-5)))))))
			else
				set YR to ((YR - 2000) / 1000)
				set JDE0 to 2.45190005952E+6 + (YR * 3.6524274049E+5 + (YR * (-0.06223 + (YR * (-0.00823 + (YR * 3.2E-4))))))
			end if
			return JDE0
		end ApproxWinterSolstice
	end script
	
	on run
		set suggestion to (current date)'s year as string
		set wanted_year to inputYear(suggestion)
		set {equinoxTable, quarterTable} to {{}, {}}
		# We do compute the previous Winter Solstice.
		# We also keep all times in UT, and converts at the end.
		
		set end of equinoxTable to my equinoxlib's jdeforwintersolstice(wanted_year - 1)
		set end of equinoxTable to my equinoxlib's jdeforvernalequinox(wanted_year)
		set end of equinoxTable to my equinoxlib's JDEForSummerSolstice(wanted_year)
		set end of equinoxTable to my equinoxlib's JDEForAuturnalEquinox(wanted_year)
		set end of equinoxTable to my equinoxlib's jdeforwintersolstice(wanted_year)
		
		repeat with q from 1 to 4
			# compute Arithmetic Mean and prepare delta
			set JDMean to ((item q of equinoxTable) + (item (q + 1) of equinoxTable)) / 2
			set {delta, n} to {60, 1}
			# Use binary search, to iterate towards the solution.
			repeat while true
				set curLong to my solaris's longitude(JDMean)
				if (abs ((item q of quarterAngles) - curLong)) < 0.1 then
					set end of quarterTable to JDMean
					exit repeat
				end if
				if curLong > (item q of quarterAngles) then
					set JDMean to JDMean - (delta / (n * 2))
				else
					set JDMean to JDMean + (delta / (n * 2))
				end if
				set n to n + 1
			end repeat
		end repeat
		_MAKEDOC((scriptTitle & " For " & wanted_year as text) & "." & linefeed & linefeed, 600, 300)
		repeat with i from 1 to 8
			local yr_, mnt_, day_, h_, m_, s_
			if (i mod 2 = 1) then # odd - Quarter!
				set {yr_, mnt_, day_, h_, m_, s_} to adjustJulianDateAccordingToGMT(item (i div 2 + 1) of quarterTable)
				_PRINT(formatDate(yr_, mnt_, day_) & " : " & item (i div 2 + 1) of quarterNames & linefeed)
			else
				set {yr_, mnt_, day_, h_, m_, s_} to adjustJulianDateAccordingToGMT(item (i div 2 + 1) of equinoxTable)
				_PRINT(formatDate(yr_, mnt_, day_) & " : " & item (i div 2) of equinoxNames & linefeed)
			end if
			
		end repeat
		
	end run
	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 lastdayInMonth(YR, mo)
		local daynumbers
		set daynumbers to {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}
		if mo < 1 or mo > 12 then
			error "lastDay: Month out of range..." number 3010
		else if mo = 2 then
			if ((YR mod 4) = 0) and (((YR mod 100) ≠ 0) or ((YR mod 400) = 0)) and (mo > 2) then
				set item 2 of daynumbers to 29
			end if
		end if
		return item mo of daynumbers
	end lastdayInMonth
	
	on set_date_long(YR, mnth, dy, hr, min, sec)
		# NG http://macscripter.net/viewtopic.php?id=41374 there is a short version there too.
		tell (current date) to set {day, year, its month, day, its hours, its minutes, its seconds, cur_date} to {1, YR, mnth, dy, hr, min, sec, it}
		return cur_date
		
	end set_date_long
	
	# returns a time triplet, from a Julian day number,
	# Adjusted to return UT.
	
	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, tempvar, _hours, _minutes, _seconds
		set timeFraction to (JDN mod 1) * 24
		set tempvar to (timeFraction div 1 as integer)
		set _hours to (tempvar + 12) mod 24
		set timeFraction to (timeFraction - tempvar)
		set _minutes to timeFraction * 60 div 1
		set _seconds to (timeFraction - (_minutes / 60)) * 3600
		return {_hours, _minutes, _seconds}
	end timeTriplet
	
	on adjustJulianDateAccordingToGMT(JDE)
		local dt, _yr, mnt, _day, _hr, _min, _sec
		# returns a sextet
		set JDE to JDE + 0.5
		
		# Find day and time of the event so far in ut time.
		
		set {_yr, mnt, _day} to gregorianCalendarTriplet from JDE
		set {_hr, _min, _sec} to timeTriplet from JDE
		# we have to adjust the hours
		if _hr > 12 then
			set _hr to _hr - 12
		else
			set _hr to _hr + 12
		end if
		
		# figure out the TimeToGMT at this instace
		
		set dt to set_date_long(_yr, mnt, _day, _hr, _min, _sec)
		
		set _hr to _hr + (civilTimeToGmtFromLocalTZInSeconds(dt) / 3600)
		
		if _hr > 24 then
			set _hr to _hr - 24
			set _day to _day + 1
			if not _day ≤ lastdayInMonth(_yr, mnt) then
				if mnt < 12 then
					set _day to 1
					set mnt to mnt + 1
				else
					set _yr to _yr + 1
					set mnt to 1
					set _day to 1
				end if
			end if
		else if _hr < 0 then
			set _hr to _hr + 24
			if _day > 1 then
				set _day to _day - 1
			else if mnt > 1 then
				set mnt to mnt - 1
				set _day to lastdayInMonth(_yr, mnt)
			else
				set _yr to _yr - 1
				set mnt to 12
				set _day to 31
			end if
		end if
		return {_yr, mnt, _day, _hr, _min, _sec}
	end adjustJulianDateAccordingToGMT
	
	# 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 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 _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 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 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
	
	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 _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 inputYear(suggestion)
		local stdMsg, probe, theMessage, anicon
		set stdMsg to "Enter the year you want to see the Solar Phases for."
		set theMessage to stdMsg
		set anicon to note
		repeat while true
			tell application (path to frontmost application as text)
				try
					set theDigitstring to text returned of (display dialog theMessage default answer suggestion with title my scriptTitle with icon anicon buttons {"Quit", "Ok"} cancel button 1 default button 2)
				on error
					error number -128
				end try
			end tell
			try
				set probe to theDigitstring as number
				exit repeat
			on error
				set theMessage to "\"" & theDigit & "\" isn't a valid year!" & return & stdMsg
				set anicon to caution
			end try
		end repeat
		return theDigitstring
	end inputYear
	
end script
tell solarEvents to run