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 20139 days, 8 hours, 4 minutes ago.