# Calculate sunrise and sunset times based on an algorithm that
# originally appeared in QST, June 84, page 56.  Transcribed from
# a BASIC program by N4TL based on calculations originally done
# by WS4G and SM0CCN.

# Usage:  sun [ month day [latitude longitude]]

# If latitude or longitude are omitted, this script will look for
# global variables named LATITUDE and LONGITUDE.  If the month
# or day is omitted, today's date will be used.

proc declination(month, day)
	# Approximation for sun declination
	local dec
	@ dec = day + 30.3*(month - 1)
	@ dec = 10 + .988*dec
	@ dec = -23.5*cos(dec*pi/180)
	return dec
end

proc fnacs(x)
	return pi/2 - atan(x/sqrt(1-x*x))
end

proc printtime(name, timeofday, suffix)
	local hour, min
	if (timeofday < 0) @ timeofday += 24
	if (timeofday > 24) @ timeofday -= 24
	@ hour = floor(timeofday)
	@ min = round((timeofday - hour)*60)
	if (min == 60) then
		@ min = 0
		@ hour++
	end
	calc printf("%s is %d:%02d%s", name, hour, min, suffix)
end

proc sunrise(month, day, latitude, longitude)
	local i, x, y, sunrise, sunset, length, dec, Month
	set Month = january february march april may june ^
		july august september october november december
	if (latitude == "" || latitude == "") then
		if ($?LONGITUDE && $?LATITUDE) then
			@ latitude = LATITUDE
			@ longitude = LONGITUDE
		else
			echo Longitude and Latitude must be specified
			return
		end
	end
	if (month == "" || day == "") then
		local dt
		set dt = `dt`
		set month = $dt:1
		set day = $dt:2
	end
	if (month !~ '[0-9]*') then
		local m
		@ m = "$lower(month)*"
		for i = 0 to 11 do
			if (Month[i] =~ m) break;
		end
		@ month = i + 1
	end
	@ dec = declination(month, day)
	@ x = longitude/15
	@ y = tan(dec*pi/180)*tan(latitude*pi/180)
	if (y > 1 || y < -1) then
		echo No sunrise or sunset
	else
		@ z = fnacs(y)*180/pi/15
		@ sunrise = x + z
		@ sunset = x - z
		@ length = sunset - sunrise
		printtime Sunrise $sunrise " UTC"
		printtime Sunset $sunset " UTC"
		printtime "Length of day" $length
	end
end

sunrise $argv
