More on Astronomy:
Solar Eclipses

José Martín-García & Jeff Bryant,
Wolfram Research

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

ABSTRACT

Wolfram Language 14 introduces new functionality to compute everything about solar eclipses, like the upcoming total eclipse of April 8, 2024. The updated function SolarEclipse can compute lots of things about more than 72 thousand solar eclipses in 30 thousand years, with high precision. This includes both global properties of the events and local properties as observed by particular observers anywhere on Earth. Analyzing eclipses requires very good control of the position of the Moon, so we have also added new functions to report lunation numbers, phases of the moon and more. This talk will demonstrate these new capabilities of the WL Astro functionality, with multiple visualization examples to help understanding the time evolution and geometry of any type of solar eclipse.

solarEclipseAnimate[eclipse_,loc_]:=DynamicModule[{elems,T0,Tmax,hmax,x,y,dec,mu,l1,l2,tanf1,tanf2,b},
elems=SolarEclipse[eclipse,"BesselianElementsFunctions"];
T0=elems["ReferenceDate"];
Tmax=SolarEclipse[eclipse,"MaximumEclipseDate"];
hmax=QuantityMagnitude[Tmax-T0,"Hours"];
{x,y,dec,mu,l1,l2,tanf1,tanf2}=Rest[Values[elems]];
b=1-GeodesyData["ITRF00","Flattening"];
Manipulate[GraphicsRow[{Graphics[{
LightBlue,Disk[{0,0},{1,Sqrt[Sin[dec[t]]^2+b^2 Cos[dec[t]]^2]}],
Black,Point[{x[t],y[t]}],
Opacity[0.1],Disk[{x[t],y[t]},l1[t]],
Opacity[0.5],Disk[{x[t],y[t]},Abs@l2[t]]
},PlotRange->1.5{{-1,1},{-1,1}},ImageSize->1000],
Graphics[{
LightBlue,Disk[{0,0},{1,Sqrt[Sin[dec[t]]^2+b^2 Cos[dec[t]]^2]}],
Black,Point[{x[t],y[t]}],
Opacity[0.1],Disk[{x[t],y[t]},l1[t]],
Opacity[0.5],Disk[{x[t],y[t]},Abs@l2[t]]
},PlotRange->{x[t],y[t]}+0.05{{-1,1},{-1,1}},ImageSize->1000],
AstroGraphics[(Entity["Sun", Entity["Star","Sun"], "Entity[\"Star\", \"Sun\"]", "star"]),AstroReferenceFrame->{"Equatorial",If[loc==="ShadowAxis",Replace[SolarEclipse[eclipse,{"ShadowAxisPosition",eclipse+Quantity[t-hmax,"Hours"]}],_Missing:>Here],loc],eclipse+Quantity[t-hmax,"Hours"]},AstroRange->Quantity[1, "AngularDegrees"],ImageSize->1000]},ImageSize->1000],{{t,hmax},-3+hmax,3+hmax},SaveDefinitions->True]
];
$Line=-1;
""

Summary

AstroPosition (precise ephemeris) and AstroGraphics (sky maps).

Update of SolarEclipse.

Besselian Elements.

Conclusions and future work.

1. AstroPosition & AstroGraphics

1.1. Precise Ephemeris Data

• Our solar system computations use NASA/JPL DE441 data:

AstroPosition["Sun"]

Dynamic[AstroPosition[(Entity["Sun", Entity["Star","Sun"], "Entity[\"Star\", \"Sun\"]", "star"])],UpdateInterval->0.1]

AstroPosition[(Entity["Moon", Entity["PlanetaryMoon","Moon"], "Entity[\"PlanetaryMoon\", \"Moon\"]", "planetary moon"]),"Equatorial"]

Dynamic[AstroPosition["Moon","Equatorial"],UpdateInterval->0.1]

1.2. Sky Maps with AstroGraphics

• Map of the region around Orion:

AstroGraphics[(Entity["Betelgeuse", Entity["Star","Betelgeuse"], "Entity[\"Star\", \"Betelgeuse\"]", "star"]),AstroReferenceFrame->"Equatorial",AstroRange->(QuantityPostfix[25, "°", degrees, "AngularDegrees"])]

• Add constellation illustrations:

AstroGraphics[(Entity["Betelgeuse", Entity["Star","Betelgeuse"], "Entity[\"Star\", \"Betelgeuse\"]", "star"]),AstroReferenceFrame->"Equatorial",AstroRange->(QuantityPostfix[25, "°", degrees, "AngularDegrees"]),AstroBackground->AstroStyling[{"BlackSky","ShowConstellationIllustrations"->True}]]

• Choose a different area of the sky:

AstroGraphics[{White,PointSize[Large],Point@AstroPosition[{0,0},"Galactic"]},AstroRange->(QuantityPostfix[25, "°", degrees, "AngularDegrees"]),AstroReferenceFrame->"Galactic",AstroBackground->"GalacticSky",AstroGridLines->Automatic]

1.3. Solar Eclipses

• Visualize a solar eclipse from a location on Earth:

date=SolarEclipse[]

Manipulate[AstroGraphics[(Entity["Sun", Entity["Star","Sun"], "Entity[\"Star\", \"Sun\"]", "star"]),AstroReferenceFrame->{"Equatorial",Here,date+Quantity[h,"Hours"]},AstroRange->Quantity[1, "AngularDegrees"]],{h,-1,3},SaveDefinitions->True]

• Observe from the Sun:

Manipulate[AstroGraphics[(Entity["Earth", Entity["Planet","Earth"], "Entity[\"Planet\", \"Earth\"]", "planet"]),AstroReferenceFrame->{"Ecliptic",(Entity["Sun", Entity["Star","Sun"], "Entity[\"Star\", \"Sun\"]", "star"]),date+Quantity[h,"Hours"]},AstroRange->Quantity[.005, "AngularDegrees"]],{h,-2,2},SaveDefinitions->True]

1.4. New Moons

• Find the next new and full moons:

NewMoon[]

CalendarConvert[%,"Jewish"]

FullMoon[]

• Find the current instant in the continuous numbering of new moons, having 0 as the first new moon of year 2000:

LunationNumber[]

Dynamic[LunationNumber[]//InputForm,UpdateInterval->0.1]

• Every solar eclipse is close to a new moon, but not every new moon is a solar eclipse:

newmoons2023=FromLunationNumber/@Range[285,296]

Select each one and see what happens around the instant of new moon:

Manipulate[AstroGraphics[(Entity["Sun", Entity["Star","Sun"], "Entity[\"Star\", \"Sun\"]", "star"]),AstroRange->Quantity[6,"AngularDegrees"],AstroReferenceFrame->{"Ecliptic",date+Quantity[h,"Hours"]},AstroGridLines->Automatic],{date,newmoons2023},{{h,0},-12,12},SaveDefinitions->True]

2. Rewrite of SolarEclipse

2.1. General Comments

• There are 72K+ solar eclipses between years -13200 and +17200 approx.

• SolarEclipse can now compute 80+ properties about each of them.

• Global or local properties. Dates, positions, curves, polygons, magnitudes, types, ...

• Terminology (diagram not to scale!):

In a solar eclipse the penumbra intersects the Earth, but not necessarily the umbra or the shadow axis.

• This is to scale, in millions of kilometers, the second plot being enlarged by a factor close to 400:

With[{
au=QuantityMagnitude["AstronomicalUnit","Gigameters"],
mO=QuantityMagnitude[AstroDistance[(Entity["Earth", Entity["Planet","Earth"], "Entity[\"Planet\", \"Earth\"]", "planet"]),(Entity["Moon", Entity["PlanetaryMoon","Moon"], "Entity[\"PlanetaryMoon\", \"Moon\"]", "planetary moon"])],"Gigameters"],
sR=QuantityMagnitude[(Entity["Sun", Entity["Star","Sun"], "Entity[\"Star\", \"Sun\"]", "star"])["Radius"],"Gigameters"],
mR=QuantityMagnitude[(Entity["Moon", Entity["PlanetaryMoon","Moon"], "Entity[\"PlanetaryMoon\", \"Moon\"]", "planetary moon"])["Radius"],"Gigameters"],
eR=QuantityMagnitude[(Entity["Earth", Entity["Planet","Earth"], "Entity[\"Planet\", \"Earth\"]", "planet"])["Radius"],"Gigameters"]
},
Column[{
Graphics[{
{RGBColor[1,0.7,0.1],Disk[{0,0},sR]},
{RGBColor[0.3,0.4,1],Disk[{au,0},eR]},
{Gray,Circle[{au,0},mO]},
{Gray,Disk[{au-mO,0},mR]}
},
ImageSize->Full,Frame->True],
Graphics[{
{RGBColor[0.3,0.4,1],Disk[{au,0},eR]},
{Gray,Circle[{au,0},mO]},
{Gray,Disk[{au-mO,0},mR]}
},
PlotRange->{au+{-1.025mO,0.025mO},{-0.025 mO,0.025mO}},
ImageSize->Full,Frame->True]
}
]
]

2.2. Global Properties

• Eclipses are identified by the date of maximum eclipse:

eclipse=SolarEclipse[]

• Then we can query properties about them:

SolarEclipse[eclipse,"Type"]

• Search from other dates, and in either direction:

SolarEclipse[DateObject[{1900,1,1}]]

SolarEclipse[DateObject[{1900,1,1}],TimeDirection->-1]

• Search by type:

SolarEclipse[Now,EclipseType->"Total"]

SolarEclipse[Now,EclipseType->"Annular"]

SolarEclipse[Now,EclipseType->"Hybrid"]

SolarEclipse[Now,EclipseType->"Partial"]

2.3. Eclipse Maps

2.3.1. Example of a Total Eclipse

• Map for the solar eclipse of April 8, 2024:

SolarEclipse["EclipseMap"]

• For every eclipse there is an instant of “maximum eclipse” when the shadow axis passes closest to Earth’s center.

eclipse=SolarEclipse[]

SolarEclipse[eclipse,"MaximumEclipseDate"]

%===eclipse

mloc=SolarEclipse[eclipse,"MaximumEclipsePosition"]

• Look at Earth from the Sun:

DynamicModule[{T0,xF,yF,decF,muF,l1F,l2F,tanf1F,tanf2F,P1loc,P2loc,P3loc,P4loc,C1loc,C2loc,line},
{T0,xF,yF,decF,muF,l1F,l2F,tanf1F,tanf2F}=Values@SolarEclipse[eclipse,"BesselianElementsFunctions"];
P1loc=SolarEclipse[eclipse,"ContactPointP1Position"];
P2loc=SolarEclipse[eclipse,"ContactPointP2Position"];
P3loc=SolarEclipse[eclipse,"ContactPointP3Position"];
P4loc=SolarEclipse[eclipse,"ContactPointP4Position"];
C1loc=SolarEclipse[eclipse,"ContactPointC1Position"];
C2loc=SolarEclipse[eclipse,"ContactPointC2Position"];
line=Line[{{xF[-5],yF[-5]},{xF[5],yF[5]}}];
Manipulate[With[{date=T0+Quantity[h,"Hours"]},GeoGraphics[{White,SolarEclipse[eclipse,"ShadowAxisLine"],SolarEclipse[eclipse,"UmbraEnvelopePolygon"],Opacity[0.5,Gray],Disk[{xF[h],yF[h]},l1F[h]],Opacity[1],PointSize[Large],Red,Point[{P1loc,P2loc,P3loc,P4loc}],Orange,Point[{C1loc,C2loc,mloc}],Black,line,Opacity[0.95,Black],Disk[{xF[h],yF[h]},Abs[SolarEclipse[eclipse,{"BesselianCapitalL2",date,SolarEclipse[eclipse,{"ShadowAxisPosition",date}]}]]]},GeoRange->"World",GeoProjection->"Orthographic",GeoCenter->AstroSubpoint["Moon",date],GeoBackground->"Satellite",GeoGridRange->{{-2.2,2.2},{-1.2,1.2}},GeoGridLines->Automatic,ImageSize->800]],{h,-2.5,3}],SaveDefinitions->True
]

• Simplification: Evolution of the solar eclipse during a few hours around maximum eclipse:

solarEclipseAnimate[eclipse,mloc]

2.3.2. Example of an Annular Eclipse

eclipse=SolarEclipse[DateObject[{2030}],EclipseType->"Annular"]

SolarEclipse[eclipse,"EclipseMap"]

mloc=SolarEclipse[eclipse,"MaximumEclipsePosition"]

• Evolution of the solar eclipse during ±3 hours from maximum eclipse:

solarEclipseAnimate[eclipse,mloc]

2.3.3. Example of a Hybrid Eclipse

eclipse=SolarEclipse[DateObject[{1940}],EclipseType->"Hybrid"]

SolarEclipse[eclipse,"EclipseMap"]

mloc=SolarEclipse[eclipse,"MaximumEclipsePosition"]

• Eclipse magnitude as observed along the shadow axis path:

magnitude[h_?NumericQ]:=With[{date=eclipse+Quantity[h,"Hours"]},SolarEclipse[eclipse,{"LocalMagnitude",date,SolarEclipse[eclipse,{"ShadowAxisPosition",date}]}]]

Plot[magnitude[h],{h,-2,2}]

• Evolution of the solar eclipse:

solarEclipseAnimate[eclipse,"ShadowAxis"]

2.3.4. Example of a Partial Eclipse

eclipse=SolarEclipse[EclipseType->"Partial"]

SolarEclipse[eclipse,"EclipseMap"]

mloc=SolarEclipse[eclipse,"MaximumEclipsePosition"]

• Evolution of the solar eclipse:

solarEclipseAnimate[eclipse,mloc]

2.3.5. Example of a “Grazing” Eclipse

• There can be total or annular solar eclipses that are not central. The last one happened in 2014:

eclipse=SolarEclipse[DateObject[{2014}]]

SolarEclipse[eclipse,"Type"]

• There are no contact points with the shadow axis, so it’s not a central eclipse:

SolarEclipse[eclipse,"Central"]

{SolarEclipse[eclipse,"ContactPointC1Date"],SolarEclipse[eclipse,"ContactPointC2Date"]}

But there are exterior contact points with the umbra, so it is an umbral eclipse:

SolarEclipse[eclipse,"Umbral"]

{SolarEclipse[eclipse,"ContactPointU1Date"],SolarEclipse[eclipse,"ContactPointU4Date"]}

SolarEclipse[eclipse,"EclipseMap"]

mloc=SolarEclipse[eclipse,"MaximumEclipsePosition"]

• Evolution of the solar eclipse:

solarEclipseAnimate[eclipse,mloc]

2.4. Eclipse Cycles

2.4.1. Saros

• Every solar eclipse belongs to a Saros series:

eclipse=SolarEclipse[]

SolarEclipse[eclipse,"SarosSeries"]

These are the eclipses of Saros series 139:

saros139=SolarEclipse[{"Saros",139}]

Length[saros139]

TimelinePlot[saros139]

Eclipses of the same Saros series are separated by a Saros

Quantity[1,"Saros"]

UnitConvert[%,"Days"]

Around[Differences[saros139]]

The Saros is the result of this near-coincidence:

%/Quantity[1,"SynodicMonths"]

%%/Quantity[1,"DraconicMonths"]

Our eclipse is number 30 in the series:

SolarEclipse[eclipse,"IndexInSaros"]

• Show the magnitudes of the eclipses in the Saros series:

magnitudes=SolarEclipse[saros139,"Magnitude"]

DateListPlot[magnitudes,Joined->False]

• Compare the central paths of the eclipses of the series:

geographics[eclipse_]:=With[{asp=AstroSubpoint[(Entity["Sun", Entity["Star","Sun"], "Entity[\"Star\", \"Sun\"]", "star"]),eclipse]},Labeled[GeoGraphics[{GeoStyling[],Replace[SolarEclipse[eclipse,"UmbraEnvelopePolygon"],_Missing:>{}],Orange,PointSize[Large],Point[asp]},GeoRange->"World",GeoProjection->"Orthographic",GeoCenter->asp,GeoBackground->None,GeoGridLines->Automatic],eclipse]]
ListAnimate[geographics/@saros139,AnimationRate->4]

2.4.2. The Saros-Inex Panorama

• Every eclipse belongs to an Inex series:

eclipse=SolarEclipse[]

SolarEclipse[eclipse,"InexSeries"]

si=SolarEclipse[eclipse,"SarosInex"]

• Find the {Saros, Inex} pairs for thousands of eclipses:

sis=SolarEclipse[{DateObject[{1000}],DateObject[{3000}],All},"SarosInex"];

Dimensions[sis]

Graphics[{Point[sis],Red,PointSize[Large],Point[si]},Frame->True]

2.5. Local Circumstances

• Let us observe the next eclipse from our location in Champaign, IL:

eclipse=SolarEclipse[]

• We can define local maximum eclipse, the instant in which the shadow axis is closest to us:

date=SolarEclipse[eclipse,{"LocalMaximumEclipseDate",Here}]

• The magnitude of the eclipse is 0.98, not total for us. Observers in Champaign will need eclipse glasses:

SolarEclipse[eclipse,{"LocalMagnitude",date,Here}]

• Compute the contact times:

p1date=SolarEclipse[eclipse,{"LocalPenumbraContact1Date",Here}]

p2date=SolarEclipse[eclipse,{"LocalPenumbraContact2Date",Here}]

{p1date,p2date}-date

Plot the evolution of the magnitude:

magnitude[h_?NumericQ]:=SolarEclipse[eclipse,{"LocalMagnitude",date+Quantity[h,"Hours"],Here}]

Plot[magnitude[h],{h,-1.5,1.5}]

• Find where is the umbra at the time of local maximum eclipse:

date+Quantity[{-1,0,1},"Minutes"]

SolarEclipse[eclipse,{"UmbraPolygon",#}]&/@%

GeoGraphics[{GeoMarker[],%}]

3. Besselian Elements

3.1. Main Ideas

• Method introduced by Friedrich Bessel in 1824, and still regularly used today.

• 1) Construct a frame with origin at Earth’s center, Z axis aligned with the shadow axis, and X axis on Earth’s equatorial plane.

• 2) Model Earth as an ellipsoid of revolution, of semiaxes {1, b}.

• 3) Model the shadow cones by these variables:
	- Equatorial declination dec of the shadow axis.
	- Besselian position {x, y, z} of the shadow vertex.
	- Aperture angle f of the cone.

Manipulate[Graphics3D[{
{Opacity[0.9],Ellipsoid[{0,0,0},{1,1,b}]},
ShadowCone[dec,{x,y},{z,f}],
ShadowAxis[dec,{x,y}],
BesselianFrame[dec]
},PlotRange->1.2,ImageSize->Large],
{{b,0.5},0,1},
{{dec,0},-Pi/2,Pi/2},
{{x,-0.3},-1,1},
{{y,0.2},-1,1},
{{z,2},-10,10},
{{f,0.04},0,Pi/2},
Initialization:>(
BesselianFrame[d_]:={
InfinitePlane[{0,0,0},{{-Sin[d],0,Cos[d]},{0,1,0}}],
Thick,{
Arrow[1.2{-{Cos[d],0,Sin[d]},{Cos[d],0,Sin[d]}}],
Arrow[1.2{-{0,1,0},{0,1,0}}],
Arrow[1.2{-{-Sin[d],0,Cos[d]},{-Sin[d],0,Cos[d]}}]
}};ShadowAxis[dec_,{x_,y_}]:={Purple,Thick,Line[{{-1.2,x,y},{1.2,x,y}}.RotationMatrix[dec,{0,1,0}]]};ShadowCone[dec_,{x_,y_},{H_,ap_}]:={Opacity[0.5,Purple],
Cone[{{H-20,x,y},{H,x,y}}.RotationMatrix[dec,{0,1,0}],20 Tan[ap]],
Cone[{{H+20,x,y},{H,x,y}}.RotationMatrix[dec,{0,1,0}],20 Tan[ap]],
Opacity[1],PointSize[0.02],Point[{H,x,y}.RotationMatrix[dec,{0,1,0}]]
};)
]

• 4) One more number is needed: Earth’s rotation as a function of time. Concretely the Greenwich hour angle μ of the shadow axis.

• Most computations are related to the problem of solving the intersection of the ellipsoid and the cone, as well as various discriminants, time-derivative conditions, etc.

• Note that the frame moves in time:

extendedLine[p1_,p2_,e_]:=Line[{p1,p2+e(p2-p1)}]
axis[sun_,moon_]:=extendedLine[sun,moon,1];
shadow[{sun_,sunR_},{moon_,moonR_}]:=With[{dsm=Sqrt[(sun-moon).(sun-moon)]},
With[{vL=RotationMatrix[-ArcTan[(sun-moon)[[1]]/dsm]].{-(Sqrt[(dsm^2)-(((sunR-moonR))^2)])/dsm,-(sunR-moonR)/dsm},
vR=RotationMatrix[-ArcTan[(sun-moon)[[1]]/dsm]].{(Sqrt[(dsm^2)-(((sunR-moonR))^2)])/dsm,-(sunR-moonR)/dsm}},
{Opacity[0.3],
extendedLine[sun+sunR vL,moon+moonR vL,0.75],
extendedLine[sun+sunR vR,moon+moonR vR,0.75]
}
]
]
body[{pos_,r_},color_]:={Circle[pos,r],Point[pos],Opacity[0.2,color],Disk[pos,r]};
BesselianPlane[sun_,moon_,earth_]:=Line[With[{shift=3Normalize@Cross[sun-moon]},{earth-shift,earth+shift}]]
Manipulate[
sun={0,10};
sunR=2;
moon={moonx,3};
moonR=0.5;
earth={0,0};
earthR=1;
Graphics[{
body[{sun,sunR},Yellow],
body[{moon,moonR},Gray],
body[{earth,earthR},Blue],
axis[sun,moon],
shadow[{sun,sunR},{moon,moonR}],
shadow[{sun,sunR},{moon,-moonR}],
BesselianPlane[sun,moon,earth]
},
PlotRange->{{-3.5,3.5},{-1.1,12.1}},
ImageSize->300
],
{moonx,2,-2},
SaveDefinitions->True
]

3.2. Example Computation with Besselian Elements

• Compute polynomial approximations for the relevant functions:

eclipse=SolarEclipse[]

These are functions in hours with respect to a given TT reference date:

SolarEclipse[eclipse,"BesselianElementsFunctions"]

{T0,x,y,dec,mu,l1,l2,tanf1,tanf2}=Values[%];

• From these functions, find for example the instant of maximum eclipse:

{min,trule}=FindMinimum[Sqrt[x[t]^2+y[t]^2],{t,0}]

T0+Quantity[t/.trule,"Hours"]

Compare to the reported value after changing to the "TT" time system:

TimeSystemConvert[SolarEclipse[eclipse,"MaximumEclipseDate"],"TT"]

Conclusions

• Wolfram Language 14.0 adds several functions about the Moon (NewMoon, FullMoon, MoonPhaseDate, LunationNumber, FromLunationNumber), as well as AstroSubPoint.

• SolarEclipse has been rewritten, now handling more than 72200 eclipses, from -13200 to +17200 approx., the coverage of DE441.

• We will update our web [precisioneclipse.com] to support eclipses in the near future and past.

• The functionality will be extended in the future to work with lunar eclipses (LunarEclipse), as well as general occultation/transit scenarios.