eclipse.frink

Download or view 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.

    TODO:  Use procedure in Meeus to find candidate times
*/


use sun.frink

// 2023 eclipse as seen from Albuquerque, NM
startDate = #2023-10-14#
endDate   = #2023-10-15#
timezone = "US/Mountain"
lat =   35.106766 degrees North
long = 106.629181 degrees West

// 2024 eclipse as seen from Dallas, TX
// startDate = #2024-04-05#
// endDate   = #2024-04-10#
// timezone = "US/Central"
// lat =   32.779167 degrees North
// long =  96.808889 degrees West

date = startDate

first = true
while date < endDate
{
   // Find lunar eclipses
   date = fullMoon[date]  // Find full moon
   phase = moonPhaseAngle[date]

   println[date + "\t" + (phase->"degrees")]
   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
   if first
      date = startDate
   else
      date = date + 1/2 lunarmonth
   
   date = newMoon[date]  // Find new moon
   phase = moonPhaseAngle[date]
   first = false
   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."]
   }

   if phase < 12 degrees
      date = date + 1 hour
   else
      date = date + 1/2 lunarmonth

   println["Next date is $date"]
}

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
}


Download or view 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 19944 days, 15 hours, 9 minutes ago.