// Program to calculate sun alignments with Manhattan's East-West streets.

use sun.frink

North = +1
South = -1
West = +1
East = -1

lat = DMS[40,45,261.16] North
long = DMS[73,59,9.02] West

// Approximate direction of the streets.
// see http://en.wikipedia.org/wiki/Manhattanhenge
// http://www.charlespetzold.com/etc/AvenuesOfManhattan/index.html
axis = (270 + 29.008) degrees

// Conversion to coordinate system used by Meeus
axisMeeus = (axis + 180 degrees) mod circle

println[axisMeeus -> degrees]

date = #2016-01-01 12:01 PM#
dateOut = ### yyyy-MM-dd hh:mm:ss a zzz ###

sep = "\t"

while (date <= #2020-12-31#)
{
   date = sunSecantAzimuth[date, lat, long, axisMeeus]
   
   [azimuth, altitude] = refractedSunAzimuthAltitude[date, lat, long]
   actualSunset = sunset[date, lat, long]
   [actualSunsetAz, actualSunsetAlt] = refractedSunAzimuthAltitude[actualSunset, lat, long]
   actualSunsetAz = (actualSunsetAz + 180 degrees) mod circle
   print[(date -> [dateOut, "Eastern"]) + "$sep"]
   print[format[JD[date],day,5] + "$sep"]
   print[format[altitude,degrees,2] + "$sep"]
   print[format[actualSunsetAz, degrees, 2] + "$sep"]
   println[]
   
   date = date + 1 day
}