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 19980 days, 11 hours, 40 minutes ago.