--[[
Celestia CELX script by Mark Meretzky http://oit2.scps.nyu.edu/~meretzkm/
Show the sidereal time and view south towards the meridian
from the given location at hour o'clock in the time zone on the current day.
Apple screen capture: shift apple 3 to create .png.

Apple window capture:
Grab.app, Capture -> Window -> Choose Window
Click on Celestia
File -> Save As...
Double click on saved file to open it; then print it.
]]

--If it's December 1941 in Casablanca, what time is it in New York?
UT = 0		--Universal Time
EST = -5	--Eastern Standard Time is 5 hours behind UT
EDT = -4	--Eastern Daylight Time is 4 hours behind UT

sun = {}
sun.object = celestia:find("Sol")

earth = {}
earth.object = celestia:find("Sol/Earth")
earth.frame = celestia:newframe("planetographic", earth.object)

moon = {}
moon.object = celestia:find("Sol/Earth/Moon")

observer = {}
observer.object = celestia:getobserver()
observer.fov = 32	--field of view, in degrees
observer.object:setfov(math.rad(observer.fov))

--[[
Adjust loc.zone and loc.hour seasonally.
loc.hour should be when it gets dark.
]]

loc = {}	--location on earth
loc.zone = EDT
loc.hour = 12 + 9	--9:00 p.m. in the above time zone

--Andrus Planetarium, Hudson River Museum, Yonkers, New York 10701-1801 USA
loc.altitude = 30	--meters above sea level

function dms_to_r(degrees, minutes, seconds)	--to radians
	return math.rad(degrees + (minutes + seconds / 60) / 60)
end

loc.latitude = dms_to_r(40, 57, 15)	--positive north, negative south
loc.longitude = dms_to_r(-73, -53, -48)	--positive east, negtive west
loc.standard = dms_to_r(-75, 0, 0)	--standard meridian for our time zone

--Speed of light is defined to be 299,792,458 meters per second.
--See macro KM_PER_LY in src/celengine/astro.h.
micro_lightyear = 9460730.4725808	--kilometers

--Distance in micro lightyears from center of Earth
--to location on surface.
loc.r = (earth.object:radius() + loc.altitude / 1000) / micro_lightyear

--Convert one Unicode character to UTF-8.

function utf8(unicode)
	if unicode < 127 then		--0x007F
		return string.char(unicode)
	elseif unicode < 2047 then	--0x07FF
		return string.char(
			192 + math.floor(unicode / 64),
			128 + math.mod(unicode, 64)
		)
	elseif unicode < 32767 then	--0xFFFF
		return string.char(
			224 + math.floor(unicode / 4096),
			128 + math.mod(math.floor(unicode / 64), 64),
			128 + math.mod(unicode, 64)
		)
	else
		celestia:flash("Bad character " .. unicode .. " to utf8.")
	end
end

--Unicode characters encoded in UTF-8
degree_char = utf8(176)		--\u00B0	\194\176
minute_char = utf8(8242)	--\u2032	\226\128\178
second_char = utf8(8243)	--\u2033	\226\128\179

--Round to closest integer.  If tied, round up.

function round(x)
	return math.floor(.5 + x)
end

function signum(x)
	if x > 0 then
		return 1
	end
	if x < 0 then
		return -1
	end
	return 0
end

--[[
Return an x for which math.abs(f(x)) < delta.
Ideally, f(x) == 0.
The argument x is a good guess.
Start the interpolation at the two points x - epsilon, x + epsilon.
]]

function zero(f, x, delta, epsilon)
	if delta == nil then
		--Assume that delta and the return value of f are in radians.
		--within half a second of arc
		delta = math.pi / 180 / 60 / 60 / 2
	end

	if delta <= 0 then
		celestia:flash(string.format("delta %f must be positive.",
			delta), 5)
		celestia:wait(5)
		return 0
	end

	if math.abs(f(x)) < delta then
		return x
	end

	if epsilon == nil then
		--Assume that x is a julian date.
		epsilon = 2.0 / 24 / 60	--2 minutes of time
	end

	local a = x - epsilon
	local fa = f(a)
	if math.abs(fa) < delta then
		return a
	end

	local b = x + epsilon
	local fb = f(b)
	if math.abs(fb) < delta then
		return b
	end

	for i = 0, 2 do
		if fa == fb then
			celestia:flash(string.format(
				"can't interpolate: f(%f) == f(%f) == %f",
				a, b, fa), 10)
			celestia:wait(10)
			return x
		end

		--linear interpolation
		local c = (a * fb - b * fa) / (fb - fa)

		local fc = f(c)
		if math.abs(fc) < delta then
			return c
		end

		if signum(a) == signum(b) then
			if math.abs(fa) < math.abs(fb) then
				b = c
				fb = fc
			else
				a = c
				fa = fc
			end
		else
			if signum(fa) == signum(fc) then
				a = c
				fa = fc
			else
				b = c
				fb = fc
			end
		end
	end

	celestia:flash(string.format(
		"zero: more than 3 iterations, f(%f) == %f, f(%f) == %f",
		a, fa, b, fb), 10)
		celestia:wait(10)
	return c
end

--Convert radians to degrees, minutes, seconds.
--Result may be a latitude, which can be negative.

directionNS = 0	--North/South
directionEW = 1	--East/West
directionNone = 2

function dms(r, direction)
	local d = math.deg(r)
	if d < 0 then
		d = -d
	end

	local degrees = math.floor(d)
	d = d - degrees

	local minutes = math.floor(60 * d)
	d = d - minutes / 60

	local seconds = round(60 * 60 * d)
	if seconds == 60 then
		seconds = 0
		minutes = minutes + 1
		if minutes == 60 then
			minutes = 0
			degrees = degrees + 1
			if degrees == 360 then
				degrees = 0
			end
		end
	end

	local news = ""

	if hours ~= 0 or minutes ~= 0 or seconds ~= 0 then
		if r > 0 then
			if direction == directionNS then
				news = " N"
			elseif direction == directionEW then
				news = " E"
			end
		end

		if r < 0 then
			if direction == directionNS then
				news = " S"
			elseif direction == directionEW then
				news = " W"
			end
		end
	end

	return string.format("%02d%s %02d%s %02d%s%s",
		degrees, degree_char,
		minutes, minute_char,
		seconds, second_char,
		news)
end

--Convert radians to hours, minutes, seconds of right ascension.
--Right ascension is always >= 0 and < 360.

function hms(r)
	local d = math.deg(r)

	while d < 0 do
		d = d + 360
	end

	while d >= 360 do
		d = d - 360
	end

	--1 hour of right ascension is 360 / 24 degrees of arc.
	local hours = math.floor(d * 24 / 360)
	d = d - hours * 360 / 24

	--1 minute of right ascension is 360 / (24 * 60) degrees of arc.
	local minutes = math.floor(d * 24 * 60 / 360)
	d = d - minutes * 360 / (24 * 60)

	--1 second of right ascension is 360 / (24 * 60 * 60) degrees of arc.
	local seconds = round(d * 24 * 60 * 60 / 360)

	if seconds == 60 then
		seconds = 0
		minutes = minutes + 1
		if minutes == 60 then
			minutes = 0
			hours = hours + 1
			if hours == 24 then
				hours = 0
			end
		end
	end

	return string.format("%02dh %02dm %02ds", hours, minutes, seconds)
end

--[[
Convert latitude and longitude in radians
to unit vector in planetographic coordinates.

Planetographic coordinates are right-handed:
X axis points from center of planet toward latitude 0, longitude 180 E.
Y axis points from center of planet toward the north pole.
Z axis points from center of planet toward latitude 0, longitude 90 E.

Normally, the X axis points toward theta == 0.  Here the X axis
points towards theta == pi, which is why the local x has a negative sign.
]]

function ll_to_v(latitude, longitude)	--in radians
	local cosine = math.cos(latitude)
	local x = -math.cos(longitude) * cosine
	local y =  math.sin(latitude)
	local z =  math.sin(longitude) * cosine
	return celestia:newvector(x, y, z)
end

--To convert a vector to a position, add zeropos to it.
zeropos = celestia:newposition(0, 0, 0)

--in planetographic frame
loc.vector = ll_to_v(loc.latitude, loc.longitude)
loc.position = zeropos + loc.r * loc.vector

zname = {}
zname[UT]  = "UT"
zname[EST] = "EST"
zname[EDT] = "EDT"

mname = {}	--meridian
mname[false] = "a.m."
mname[true]  = "p.m."

function getlocaltime()
	return celestia:tdbtoutc(celestia:gettime() + loc.zone / 24)
end

function setlocaltime(year, month, day, hour, minute, seconds)
	local julian =
		celestia:utctotdb(year, month, day, hour, minute, seconds)
		- loc.zone / 24
	celestia:settime(julian)
	return julian
end

monthname = {
	"January",	--has subscript 1
	"February",
	"March",
	"April",
	"May",
	"June",
	"July",
	"August",
	"September",
	"October",
	"November",
	"December"
}

cumulative = {
	0,	--January
	31,	--February
	31 + 28,	--March
	31 + 28 + 31,	--April
	31 + 28 + 31 + 30,	--May
	31 + 28 + 31 + 30 + 31,	--June
	31 + 28 + 31 + 30 + 31 + 30,	--July
	31 + 28 + 31 + 30 + 31 + 30 + 31,	--August
	31 + 28 + 31 + 30 + 31 + 30 + 31 + 31,	--September
	31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30,	--October
	31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31,	--November
	31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31 + 30,	--December
}

--Julian day -0.5 to 0.5 (January 1, 4713 BC) was a Monday.

dayname = {
	"Monday",
	"Tuesday",
	"Wednesday",
	"Thursday",
	"Friday",
	"Saturday",
	"Sunday"
}

function date_to_s(julian, zone)

	--Add .0001 seconds because julian date is slightly truncated.
	local julian = julian + zone / 24 + .0001 / (24 * 60 * 60)
	local t = celestia:tdbtoutc(julian)
	local pm = t.hour >= 12

	if pm then
		t.hour = t.hour - 12
	end

	if t.hour == 0 then
		t.hour = 12
	end

	return string.format("%d:%02d:%02d %s %s, %s %d, %d %s",
		t.hour, t.minute, t.seconds, mname[pm],
		dayname[1 + math.mod(round(julian), 7)],
		monthname[t.month], t.day, t.year, zname[zone])
end

--[[
Convert cartesian to spherical coordinates.
Y axis points up.
theta is the longitude, phi is the latitude (not the colatitude).
]]

function xyz_to_tp(x, y, z)
	local r = math.sqrt(x * x + z * z)
	local phi = math.atan2(y, r)
	local theta = math.atan2(z, x)
	if theta < 0 then
		theta = theta + 2 * math.pi
	end
	return theta, phi
end

--Unit vector in universal coordinates from center of earth towards celestial
--north pole.  Planetographic Y axis is earth's axis.

polaris_up     = celestia:newposition(0, 1, 0)
polaris_center = celestia:newposition(0, 0, 0)

function polaris(julian)
	return earth.frame:from(polaris_up,     julian)
	     - earth.frame:from(polaris_center, julian)
end

--[[
Convert a vector in universal coordinates
to right ascension, declination in radians in equatorial coordinates.

Universal coordinates are right-handed:
X axis points from origin through the sun toward the vernal equinox in Pisces.
Y axis points from origin toward north pole of ecliptic in Draco.
Z axis points from origin toward the winter solstice in Sagittarius.

Normally, the Y axis points toward theta == pi/2.  Here we have the Z axis
pointing towards RA 18h, which is why the local z has a negative sign.

Equatorial coordinates are right-handed:
X axis points from origin toward the vernal equinox in Pisces.
Y axis points from origin toward the north celesial pole in Ursa Minor.
Z axis points from origin toward the RA 18h Dec 0.
]]

function v_to_rd(v, julian)
	local pole = polaris(julian)

	--Tilt of earth's axis in radians (about 23 degrees).
	--See "Obliquity" in data/solarsys.ssc file.
	local tilt = math.atan2(-pole:getz(), pole:gety())

	--Transform v from universal coordinates to equatorial coordinates.
	local rot = celestia:newrotation(celestia:newvector(-1, 0, 0), tilt)
	local v = rot:transform(v)

	return xyz_to_tp(v:getx(), v:gety(), -v:getz())
end

--[[
Return the right ascension and declination in radians of the body
as seen from position p in earth.frame at Julian date julian.
]]

function rd(body, p, julian)
	local bodypos = body:getposition(julian)
	local upos = earth.frame:from(p, julian)
	return v_to_rd(bodypos - upos, julian)
end

--[[
Return the altitude and azimuth (eastward from the due north point on the
horizon) in radians of the body as seen from loc on the Julian date.
At the north pole of the earth, the altitude and azimuth of a body are
the same as the planetographic coordinates of the body.
If our location is not at latitude 90 N, rot1 takes us to the correct latitude.
If our location is not at longitude 0, rot2 takes us to the correct longitude.
]]

latty = loc.latitude

rot1 = celestia:newrotation(celestia:newvector(0, 0, 1), math.pi / 2 - latty)

rot2 = celestia:newrotation(
	celestia:newvector(math.cos(latty), math.sin(latty), 0),
	loc.longitude
)

aarot = rot1 * rot2	--product of two quaternions

function aa(body, julian)
	local bodypos = earth.frame:to(body:getposition(julian), julian)
	local v = aarot:transform(bodypos - loc.position)
	local azimuth, altitude = xyz_to_tp(v:getx(), v:gety(), v:getz())
	return altitude, azimuth
end

--Return the sidereal time (the hour angle of the vernal equinox) in radians
--at the given longitude and julian date.

vernal = celestia:newvector(1, 0, 0)

function sidereal(longitude, julian)		--in radians
	local position = earth.object:getposition(julian) + vernal
	position = earth.frame:to(position, julian)
	local angle = longitude - math.atan2(position:getz(), -position:getx())
	while angle < 0 do
		angle = angle + 2 * math.pi
	end
	while angle >= 2 * math.pi do
		angle = angle - 2 * math.pi
	end
	return angle
end

forward = false
backward = false

function callback(str)
	if str == "F" then
		forward = true
		return true
	end

	if str == "B" then
		backward = true
		return true
	end

	return false
end

celestia_keyboard_callback = callback

--[[
Return negative if Sun is to left of meridian, positive if to right.
]]

function sun_radians_before_meridian(julian)
	local altitude, azimuth = aa(sun.object, julian)
	return azimuth - math.pi
end

function sun_page()

	--Guess that sun is on the meridian at noon of local apparent time
	--on the current day.

	local t = getlocaltime()
	t.hour = 12
	if loc.zone == EDT then
		t.hour = t.hour + 1
	end

	--equation of time
	local n = cumulative[t.month] + t.day
	local b = 2 * math.pi * (n - 81) / 364
	local minutes =
		9.87 * math.sin(2 * b) - 7.53 * math.cos(b) - 1.5 * math.sin(b)

	local julian =
		celestia:utctotdb(t.year, t.month, t.day, t.hour, 0, 0.0)
		+ (loc.standard - loc.longitude) / (2 * math.pi)
		- minutes / 60 / 24
		- loc.zone / 24

	julian = zero(sun_radians_before_meridian, julian)
	celestia:settime(julian)

	--Go up a few kilometers up to get a dark sky.
	local clear = (earth.object:radius() + 6) / micro_lightyear

	observer.object:setposition(earth.frame:from(
		zeropos + clear * loc.vector, julian))

	--position of target at which we are looking

	altitude, azimuth = aa(sun.object, julian)

	local targetpos = zeropos + clear * loc.vector
		+ ll_to_v(
			loc.latitude - math.pi / 2 + altitude,
			loc.longitude
		)

	observer.object:lookat(
		earth.frame:from(targetpos, julian),
		polaris(julian)
	)

	--Get the sidereal time.
	local ra, dec = rd(sun.object, loc.position, julian)

	local s = string.format(
		"%s on meridian at RA %s, Dec %s; alt %s, az %s\n"
		.. "at sidereal time %s, "
		.. "Julian date %f,\n%s\n%s.\n\n"
		.. "Zeiss setup, page %d of %d.  "
		.. "Uppercase F and B to go forward and back.\n"
		.. "View is due south toward the meridian\n"
		.. "from Andrus, "
		.. "latitude %s, longitude %s, "
		.. "alt %d m.",

		"Sun",
		hms(ra), dms(dec, directionNS),
		dms(altitude, directionNone), dms(azimuth, directionNone),

		hms(sidereal(loc.longitude, julian)),
		julian,
		date_to_s(julian, loc.zone), date_to_s(julian, UT),

		page, npages,
		dms(loc.latitude, directionNS),
		dms(loc.longitude, directionEW),
		loc.altitude
	)

        --[[
        Output file deposited in
        Unix: /home/m/meretzkm/public_html/planetarium/c/share/celestia
	Vista: C:\Users\mark\AppData\Local\VirtualStore\Program Files
		(x86)\Celestia
        ]]

	--[[

	local handle, message =
		io.open("/Users/itslogin/Desktop/sunpage.txt", "a")	--Mac
	if handle == nil then
		celestia:flash("could not open file: " .. message)
	else
		handle:write(s)
		handle:close()
	end
	]]

	celestia:print(
		s,
		2147483647,	--duration: max value for signed 32-bit integer
		-1, -1,		--origin at lower left
		1, 10   	--indent one character
	)
end

--[[
Return negative if Moon is to left of meridian, positive if to right.
]]

function moon_radians_before_meridian(julian)
	local altitude, azimuth = aa(moon.object, julian)
	return azimuth - math.pi
end

function moon_page()

	--Go to noon of the current day.

	local t = getlocaltime()
	t.hour = 12
	if loc.zone == EDT then
		t.hour = t.hour + 1
	end
	local julian = setlocaltime(t.year, t.month, t.day, t.hour, 0, 0.0)

	--Guess the next time that the moon will be on the meridian.

	local right, dec = rd(moon.object, zeropos, julian)	--geocentric
	local st = sidereal(loc.longitude, julian)
	local month = 27.321582		--days per sidereal month
	local theta = right - st
	if theta < 0 then
		theta = theta + 2 * math.pi
	end
	if theta >= 2 * math.pi then
		theta = theta - 2 * math.pi
	end

	julian = julian + (theta / (2 * math.pi)) *
		(1 + (1 + (1 + 1 / month) / month) / month)

	julian = zero(moon_radians_before_meridian, julian)
	celestia:settime(julian)

	--Go up a few kilometers up to get a dark sky.
	local clear = (earth.object:radius() + 6) / micro_lightyear

	observer.object:setposition(earth.frame:from(
		zeropos + clear * loc.vector, julian))

	--position of target at which we are looking

	altitude, azimuth = aa(moon.object, julian)

	local targetpos = zeropos + clear * loc.vector
		+ ll_to_v(
			loc.latitude - math.pi / 2 + altitude,
			loc.longitude
		)

	observer.object:lookat(
		earth.frame:from(targetpos, julian),
		polaris(julian)
	)

	local ra, dec = rd(moon.object, loc.position, julian)

	local s = string.format(
		"%s on meridian at RA %s, Dec %s; alt %s, az %s\n"
		.. "at sidereal time %s, "
		.. "Julian date %f,\n%s\n%s.\n\n"
		.. "Zeiss setup, page %d of %d.  "
		.. "Uppercase F and B to go forward and back.\n"
		.. "View is due south toward the meridian\n"
		.. "from Andrus, "
		.. "latitude %s, longitude %s, "
		.. "alt %d m.",

		"Moon",
		hms(ra), dms(dec, directionNS),
		dms(altitude, directionNone), dms(azimuth, directionNone),

		hms(sidereal(loc.longitude, julian)),
		julian,
		date_to_s(julian, loc.zone), date_to_s(julian, UT),

		page, npages,
		dms(loc.latitude, directionNS),
		dms(loc.longitude, directionEW),
		loc.altitude
	)

	celestia:print(
		s,
		2147483647,	--duration: max value for signed 32-bit integer
		-1, -1,		--origin at lower left
		1, 10   	--indent one character
	)
end

planet = {
	"Sol",
	"Sol/Earth/Moon",

	"Sol/Mercury",
	"Sol/Venus",
	"Sol/Mars",
	"Sol/Jupiter",
	"Sol/Saturn",
	"Sol/Uranus",
	"Sol/Neptune",
	"Sol/Pluto"
}

constellation = {
	"Pisces",
	"Aries",
	"Taurus",
	"Gemini",
	"Cancer",
	"Leo",
	"Virgo",
	"Libra",
	"Scorpius",
	"Sagittarius",
	"Capricornus",
	"Aquarius"
}

function dark_page()
	local t = getlocaltime()
	local julian = setlocaltime(t.year, t.month, t.day, loc.hour, 0, 0.0)
	local position = zeropos + loc.r * loc.vector;
	local v = ll_to_v(0, loc.longitude)

	--Look at the point where the celestial equator crosses the meridian.
	observer.object:setposition(earth.frame:from(position, julian))

	observer.object:lookat(
		earth.frame:from(position + v, julian),
		polaris(julian)
	)

	local s = string.format(
		"Sidereal time %s, "
		.. "Julian date %f,\n%s\n%s.\n\n",

		hms(sidereal(loc.longitude, julian)),
		julian, date_to_s(julian, loc.zone), date_to_s(julian, UT)
	)

	for i, name in ipairs(planet) do
		local ra, dec = rd(celestia:find(name), loc.position, julian)
		local c = 1 + math.floor(12 * ra / (2 * math.pi))

		first, last, shortname = string.find(name, "(%a+)$")
		s = s .. string.format("RA %s, Dec %s, %s in %s\n",
			hms(ra), dms(dec, directionNS),
			shortname, constellation[c]
		)
	end

	s = s .. string.format(
		  "\nZeiss setup, page %d of %d.  "
		.. "Uppercase F and B to go forward and back.\n"
		.. "View is due south toward the meridian\n"
		.. "from Andrus, "
		.. "latitude %s, longitude %s, "
		.. "alt %d m.",

		page, npages,
		dms(loc.latitude, directionNS),
		dms(loc.longitude, directionEW),
		loc.altitude
	)

	celestia:print(
		s,
		2147483647,	--duration: max value for signed 32-bit integer
		-1, -1,		--origin at lower left
		1, 20   	--indent one character
	)
end

npages = 3
page = 1

function main()
	celestia:show(
		"boundaries",
		"constellations",
		"galaxies",
		"grid",
		"orbits",
		"planets",
		"stars"
	)
	celestia:showlabel("constellations", "galaxies", "planets", "stars")
	celestia:select(earth.object)
	observer.object:setframe(earth.frame)

	--[[
	celestia:requestsystemaccess()
	wait(0)
	]]
	celestia:requestkeyboard(true)

	local oldpage = page - 1
	local real_julian = celestia:gettime()
	celestia:settimescale(0)		--freeze time

	while true do
		if forward then
			forward = false		--Don't let the key repeat.
			page = page + 1
			if page > npages then
				page = 1
			end
		end

		if backward then
			backward = false
			page = page - 1
			if page <= 0 then
				page = npages
			end
		end

		if page ~= oldpage then
			celestia:settime(real_julian)
			if page == 1 then
				sun_page()
			elseif page == 2 then
				moon_page()
			else
				dark_page()
			end
		end
		oldpage = page

		wait(1)
	end
end

--Call the main function.
main()