// This program animates the occultation of Mars behind the Moon

use planets.frink
use sun.frink

timezone = "US/Mountain"
lat = 40 degrees North
long = 105 degrees West

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

// This graphics object is the animated version
g = new graphics
g.font["Monospaced", .008 degree]
g.antialiasedText[false]
anim = new Animation[.3 s]
win = g.show[800,800]

//for planet = [Planet.Mars, Planet.Venus]
planet = Planet.Mars
startdate = #2025-01-13 05:00 PM#
startdate = beginningOfMinute[sunset[startdate, lat, long]]
enddate = startdate + 5 hours

println["Sunset is " + (startdate->timezone)]

frame = 0
staticsteps = 4
g2 = new graphics
g2.font["Monospaced", .008 degree]
g2.backgroundColor[.8, .8, 1]
// Draw the moon
g2.color[.9,.9,.9] 
moonRadius =   moonRadiusAngle[startdate]
g2.fillEllipseCenter[0, 0, 2 moonRadius, 2 moonRadius]

for d = startdate to enddate step 30 s
{
   [azmoon, altmoon] = refractedMoonAzimuthAltitude[d, lat, long]
   [azp, altp]     = planet.refractedAzimuthAltitude[d, lat, long]

   separation = angularSeparation[azmoon, altmoon, azp, altp]

   moonRadius =   moonRadiusAngle[d]
   marsRadius = planet.radiusAngle[d]
   //println[marsRadius -> [0, "deg", "arcmin", "arcsec"]]

   flag = ""
   if separation < (moonRadius + marsRadius) * 1.15 // Near moon?
   {
      flag = "*"
      
      if separation < moonRadius + marsRadius  // First contact
         flag = "**"

      g = new graphics
      [skyr,skyg,skyb] = skyDarkness[d, lat, long, [[.1,.1,.1], [.1,.1,.2], [.1,.1,.3], [.1,.1,.4], [.8,.8,1]]]
      g.backgroundColor[skyr, skyg, skyb]
      
      // Draw the moon
      g.color[1,1,0]  // Yellow
      g.fillEllipseCenter[(azmoon-azmoon) cos[-altmoon], -(altmoon-altmoon), 2 moonRadius, 2 moonRadius]

      // Draw Mars
      g.color[.5,.5,.5]
      g.fillEllipseCenter[(azp-azmoon) cos[-altp], -(altp-altmoon), 2 marsRadius, 2 marsRadius]
//      g.add[planet.drawPolygonRelativeToZenith[d, lat, long, azp, -altp, marsRadius, true]]
      g.text[d->[tf, timezone], (azp-azmoon) cos[-altp], -(altp-altmoon), "center", "top"]
      win.replaceGraphics[g]
      anim.add[g]
      
      if frame mod staticsteps == 0  OR  flag == "*"   // Step every 1 steps or when nearly touching
      {
         g2.color[0,0,0]
         g2.fillEllipseCenter[(azp-azmoon) cos[-altp], -(altp-altmoon), 2 marsRadius, 2 marsRadius]
         g2.text[d->[tf, timezone], (azp-azmoon) cos[-altp], -(altp-altmoon), "center", "top"]
      }
   }
   
   println[(d->timezone) + "\t" + format[separation,"deg", 5] + "\t$flag" ]

   frame = frame + 1
}

g2.show[]
filename = "MarsOccultation" + (startdate -> datef) + "_$timezone"
filename =~ %s/\//_/g  // Turn slashes into underscores in filemane.
print["Rendering animation..."]
anim.write["$filename.gif", 700, 700]
println["done."]

g2.write["$filename.png", 1024, 1024]
g2.write["$filename.svg", 1024, 1024]
g2.write["$filename.html", 1024, 1024]
browse["file:$filename.gif"]