eclipse.frink

View or download eclipse.frink in plain text format


/** This program calculates lunar and solar eclipses for the specified place
    on earth. .  It first performs a rough calculation to see if an eclipse
    could possibly occur, and then calculates the eclipse at higher quality.
*/


use sun.frink

startDate = #2017-08-01#
endDate   = #2017-08-31#

timezone = "US/Mountain"
lat =   39.60 degrees North
long = 104.98 degrees West

date = startDate

while date < endDate
{
   // Find lunar eclipses
   date = fullMoon[date]  // Find full moon
   phase = moonPhaseAngle[date]
   if phase < 1.59 degrees
   {
      [az, alt] = refractedMoonAzimuthAltitude[date, lat, long]
      print["Lunar: " + (date->timezone) + "\tAngle=" + format[phase, "deg", 4]]
      if alt > 0 deg
         println["\t+"]
      else
         println["\t-"]
   }

   // Find solar eclipses
   date = date + 1/2 lunarmonth
   date = newMoon[date]  // Find new moon
   phase = moonPhaseAngle[date]
   if phase > 178.41 degrees
   {
      [az, alt] = refractedSunAzimuthAltitude[date, lat, long]
      print["Solar: " + (date->timezone) + "\tAngle=" + format[180 degrees-phase, "deg", 4]]
      if alt > 0 deg
         println["\t+"]
      else
         println["\t-"]

//      ti = new interval[date-2 hours, date+2 hours]
//      println["Angular separation " + (sunMoonAngularSeparation[date, lat, long]->"deg")]
//      println["Angular separation " + (sunMoonAngularSeparation[ti, lat, long]->"deg")]
      sunrise = sunrise[date, lat, long]
      sunset = sunset[date, lat, long]
//      f = {|date| sunMoonAngularSeparation[date, lat, long]}
//      min = secantMinimize[f, sunrise, sunset, 1 s]
//      println[min]

      mintime = now[]
      minsep = 180 degrees
      maxPercent = 0
      contactMade = false

      anim = new Animation[1/30 s]
      
      DAY:
      for t = sunrise to sunset step 20 s
      {
         sunRadius = sunRadiusAngle[t]
         moonRadius = moonRadiusAngle[t]
         sep = sunMoonAngularSeparation[t, lat, long]
         percentEclipse = solarEclipseRatio[sunRadius, moonRadius, sep]
         if sep < minsep
         {
            minsep = sep
            mintime = t
            maxPercent = percentEclipse
         }
         
         if sep < sunRadius+moonRadius
         {
            contactMade = true
            [az, alt] = refractedSunAzimuthAltitude[t, lat, long]
            az = (az+180 degrees) mod circle
            print[(t -> timezone) + "\t" + format[sep, "deg", 5] + "\t" + format[az, "deg", 2] + "\t" + format[alt, "deg", 2] + "\t"]
            print[format[percentEclipse, percent, 1] + "%\t"]
            if (sep < abs[moonRadius-sunRadius])
               if (moonRadius < sunRadius)
                  print["* Annular"]
               else
                  print["* Total"]
            else
               print["Partial"]
            println[]
            anim.add[drawSolarEclipse[t,lat,long,true]]
         } else
            // We made contact previously this day, but not any more?  Bail out
            if contactMade == true
               break DAY
      }
      println["Minimum separation " + format[minsep, "degrees", 5] + " at " + (mintime -> timezone)]
      println["Maximum eclipse: " + format[maxPercent, percent, 2] + "%"]
      g = drawSolarEclipse[mintime,lat,long]
      g.show[]
      g.write["eclipse.png",480,480]
      
      print["Writing animation..."]
      anim.write["eclipse.gif", 320, 320]
      println["done."]
   }
   date = date + 1/2 lunarmonth
}

drawSolarEclipse[date,lat,long,forAnimation=false] :=
{
   g = new graphics

   [skyr,skyg,skyb] = skyDarkness[date, lat, long, [[.1,.1,.1], [.1,.1,.2], [.1,.1,.3], [.1,.1,.4], [.6,.6,1]]]
   g.backgroundColor[skyr, skyg, skyb]

   sunRadius = sunRadiusAngle[date]
   [sunaz, sunalt] = refractedSunAzimuthAltitude[date, lat, long]

   moonRadius = moonRadiusAngle[date]
   [moonaz, moonalt] = refractedMoonAzimuthAltitude[date, lat, long]

   if forAnimation
   {
      g.color[0,0,0,0]
      diam = 2 sunRadius + 4 moonRadius
      g.fillRectCenter[0 deg, 0 deg, diam, diam]
   }
   
   // Draw the sun
   g.color[1,1,0]  // Yellow
   g.fillEllipseCenter[(sunaz-sunaz) cos[-sunalt], -(sunalt-sunalt), 2 sunRadius, 2 sunRadius]
   
   g.color[0,0,0,.8]  
   g.fillEllipseCenter[(moonaz-sunaz) cos[-moonalt], -(moonalt-sunalt), 2 moonRadius, 2 moonRadius]
   
//   println["$sunaz\t$sunalt\t$moonaz\t$moonalt"]
   return g
}


/** This calculates the area of intersection.  See:
    http://mathworld.wolfram.com/Circle-CircleIntersection.html
*/

intersectionArea[sunRadius, moonRadius, angularSeparation] :=
{
   R = max[sunRadius, moonRadius]
   r = min[sunRadius, moonRadius]
   d = angularSeparation

   if d > R+r     // No overlap
      return 0 r^2

   if d < R-r     // Smaller circle completely inside larger?
      return pi r^2  // Return full area of smaller

   A = r^2 arccos[(d^2 + r^2 - R^2)/(2 d r)] +
       R^2 arccos[(d^2 + R^2 - r^2)/(2 d R)] -
       1/2 sqrt[(-d + r + R)(d + r - R)(d - r + R)(d + r + R)]

   return A
}

/** Calculates the percentage of total solar eclipse. */
solarEclipseRatio[sunRadius, moonRadius, angularSeparation] :=
{
   intersection = intersectionArea[sunRadius, moonRadius, angularSeparation]
   sunArea = pi sunRadius^2
   if (intersection >= sunArea)
      return 100 percent
   else
      return intersection / sunArea
}


View or download eclipse.frink in plain text format


This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 17592 days, 15 hours, 33 minutes ago.