/** 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 }