// This program draws the transit of Venus (easily adapted for Mercury)
// across the sun, taking into account a lot of effects.

use planets.frink
use sun.frink

timezone = "Mountain"
lat = 39.58560 degrees North
long = 104.89598 degrees West

/*lat = DMS[40, 33, 35.84]  North
long = DMS[105,04,57.64]  West*/


df = ### yyyy-MM-dd-HH:mm ###

g = new graphics

//for planet = [Planet.Mercury, Planet.Venus]
planet = Planet.Venus
startdate = #2012-06-05 04:04 PM#
enddate =   sunset[startdate, lat, long]
println["Sunset is $enddate"]

g = new graphics
g.backgroundColor[.8, .8, 1]
// Draw the sun
g.color[1,1,0]  // Yellow
azsun = 40 deg
altsun = 10 deg

[azsun, altsun] = refractedSunAzimuthAltitude[startdate, lat, long]
sunRadius =   sunRadiusAngle[startdate]
g.fillEllipseCenter[(azsun-azsun) cos[-altsun], -(altsun-altsun), 2 sunRadius, 2 sunRadius]
g.color[0,0,0]
win = g.show[]
g.font["SansSerif", 0.01 degree]
dateformat = ###h:mm a###

for d = startdate to enddate step 20 min
{
   [azsun, altsun] = refractedSunAzimuthAltitude[d, lat, long]
   [azp, altp]     = planet.refractedAzimuthAltitude[d, lat, long]

   separation = angularSeparation[azsun, altsun, azp, altp]

   sunRadius =   sunRadiusAngle[d]
   venusRadius = planet.radiusAngle[d]
   //println[venusRadius -> [0, "deg", "arcmin", "arcsec"]]

   flag =""
   if separation < sunRadius + venusRadius
   {
      flag = "*"
      
      if separation < sunRadius - venusRadius
         flag = "**"
   }

   azp = (azp + 180 degrees) mod circle
   azsun = (azsun + 180 degrees) mod circle
   azsep = (azp-azsun) mod (360 degrees)
   if azsep > 180 degrees
      azsep = azsep - circle
   // Draw Venus
   // How much more black could it be?  The answer is none.  None more black.
   x = azsep cos[-altp]
   y = -(altp-altsun)
   g.fillEllipseCenter[x, y, 2 venusRadius, 2 venusRadius]
   
   g.text[d -> [dateformat, timezone], x -  venusRadius, y +  venusRadius, "right", "center", 45 deg]
//      g.add[planet.drawPolygonRelativeToZenith[d, lat, long, azp, -altp, venusRadius, true]]
   win.repaint[]
   
   println[(d->timezone) + "\t" + format[separation,"deg", 5] + "\t" + format[azsep,"deg",2] + "\t$flag" ]
}

g.write["venusChart.svg", 1024, 768]
g.write["venusChart.png", 1024, 768]
g.write["venusChart.html", 1024, 768]