Solar Eclipse of
 August 21st, 2017

José Martín-García, Algorithms R&D
Jeff Bryant, Wolfram|Alpha

Join the Conversation #WolframTechConf

Wolfram PrecisionEclipse project

• Improvements on the Astro functionality.

• Improvements on the date/time-related functionality.

• Research on the theory and history of solar eclipses.

• Independent computation of contact times for the eclipse.

• Stephen’s blog on the eclipse. [blog].

• Jeff’s blogs on geographical aspects of the eclipse. [blog1], [blog2].

• The www.precisioneclipse.com microsite.

• Gathering of data and related information for a notebook that can be downloaded from the microsite.

Ingredients for eclipse prediction

1) Good ephemeris data on the position of Sun, Earth, Moon.

2) Good orientation data for the Earth.

3) Careful treatment of timing: time systems, finite speed of light, ...

Positions of Sun, Earth, Moon

• NASA provides “kernels” of data of multiple types: positions of planets, moons, spacecrafts; relations between time systems; structure of spacecrafts and their instruments; and much more. [NASA SPICE spk kernels].

• These are 3D positions in AUs in the J2000 frame centered at the Solar System Barycenter, from 12:00 to 24:00 UTC on August 21st 2017:

sunFunction=InterpolatingFunction[];

earthFunction=InterpolatingFunction[];

moonFunction=InterpolatingFunction[];

• Examples:

sunFunction[20000]

earthFunction[20000]

moonFunction[20000]

• Relevant quantities using astronomical units:

auInMeters=QuantityMagnitude[Quantity[1,"AstronomicalUnit"],"Meters"]

SunRadius=QuantityMagnitude[(Entity["Sun", Entity["Star","Sun"], "Entity[\"Star\", \"Sun\"]", "star"])["Radius"],"AstronomicalUnit"]

MoonRadius=QuantityMagnitude[(Entity["Moon", Entity["PlanetaryMoon","Moon"], "Entity[\"PlanetaryMoon\", \"Moon\"]", "planetary moon"])["Radius"],"AstronomicalUnit"]

EarthSemimajorAxis=QuantityMagnitude[GeodesyData["ITRF00","SemimajorAxis"],"AstronomicalUnit"]

EarthSemiminorAxis=QuantityMagnitude[GeodesyData["ITRF00","SemiminorAxis"],"AstronomicalUnit"]

LightSpeed=N@QuantityMagnitude[(Quantity[, "c", speed of light, "SpeedOfLight"]),"AstronomicalUnit"/"Seconds"]

Some basic physical facts

• Most of the speed of the Moon comes from Earth’s motion:

(Entity["Earth", Entity["Planet","Earth"], "Entity[\"Planet\", \"Earth\"]", "planet"])["AverageOrbitVelocity"]//UnitConvert

(Entity["Moon", Entity["PlanetaryMoon","Moon"], "Entity[\"PlanetaryMoon\", \"Moon\"]", "planetary moon"])["AverageOrbitVelocity"]//UnitConvert

Norm[moonFunction[20001]-moonFunction[20000]](Quantity[, "au", astronomical units, "AstronomicalUnit"])//UnitConvert

Norm[(moonFunction[20001]-earthFunction[20001])-(moonFunction[20000]-earthFunction[20000])](Quantity[, "au", astronomical units, "AstronomicalUnit"])//UnitConvert

• Rotations are mostly aligned.

Time

• Astronomical computations require a homogeneous axis of time. Civil time must be aligned with the rotation of the Earth, which is slowing down.

Due to tidal braking, the duration of one solar day is about 1.70 ± 0.05 ms/century longer, so each “solar century” is about 31 seconds longer than the previous one.

Sum[i dT/n,{i,n}]

UnitConvert[%/.{dT->Quantity[1.7,"Milliseconds"],n->36524}]

Most of Earth’s rotational energy is lost in heat, but part is being transferred to the Moon, which is climbing the Earth’s gravity well and getting further away a few centimeters per year.

• Time Systems:

	TAI: International Atomic Time. Homogeneous time counting SI seconds since 1958.
	
	TT = TAI + 32.184 s: Terrestrial Time. Homogeneous time, extending the homogeneous time scale previously used in Astronomy.
	
	UT1: Time that follows the rotation of the Earth, and hence is not homogeneous.
	
	UTC = TAI - 37 s: Coordinated Universal Time. Follows TAI except for an integer number of leap secs, so Abs[DUT1 = UT1 - UTC] < 1 s.

Then define:

	ΔT = TT - UT1 = 32.184 + 37 + UTC - UT1 = 69.184 - DUT1

The value 32.184 is constant.
The integer 37 changes when new leap seconds are added.
The value of DUT1 changes continuously. See [USNO]. We took the value predicted on August 9th for August 21st: 0.33687 seconds.

ΔT=32.184+37-0.33687

• In this talk all times have been already corrected to be UTC times, and will be given in the GMT time zone unless otherwise indicated.

Umbra, penumbra and antumbra

Depending on when we are we will see a different type of eclipse:

	Umbra: total eclipse.
	Penumbra: partial eclipse.
	Antumbra: annular eclipse.
	Otherwise: no eclipse.

sunUmbra[x2_, {r1_, r2_}] ={r1 (r1-r2),r1 (Sqrt[(x2^2)-(((r1-r2))^2)])}/x2;
moonUmbra[x2_, {r1_, r2_}] ={x2^2+r2(r1-r2) ,r2 (Sqrt[(x2^2)-(((r1-r2))^2)])}/x2;
sunPenumbra[x2_, {r1_, r2_}] ={r1 (r1+r2),r1 (Sqrt[(x2^2)-(((r1+r2))^2)])}/x2;
moonPenumbra[x2_, {r1_, r2_}] ={x2^2-r2 (r1+r2),-r2 (Sqrt[(x2^2)-(((r1+r2))^2)])}/x2;
Manipulate[
If[r1>r2&&r1+r2<x2,
Module[{
mU=moonUmbra[x2, {r1, r2}],
mPU=moonPenumbra[x2, {r1, r2}],
sU=sunUmbra[x2, {r1, r2}],
sPU=sunPenumbra[x2, {r1, r2}],
dirUmbra,dirPenumbra},
dirUmbra = Normalize[mU-sU];
 dirPenumbra = Normalize[mPU-sPU];
Graphics[
 {PointSize[0.01], Yellow,Disk[{0,0}, r1], LightGray,Disk[{x2,0}, r2],
Black, Circle[{0,0}, r1], Circle[{x2,0}, r2], Gray, 
Line[{sU-2r1 dirUmbra,mU+8r2 dirUmbra}], 
Line[{sPU-2r1 dirPenumbra,mPU+8r2 dirPenumbra}], 
Line[{{1,-1}sU-2r1 {1,-1} dirUmbra,{1,-1}mU+8r2 {1,-1} dirUmbra}], Line[{{1,-1}sPU-2r1 {1,-1} dirPenumbra,{1,-1}mPU+8r2 {1,-1} dirPenumbra}], 
  Purple, Line[{sU,mU}], Point[sU],Point[mU],Line[{{1,-1} sU,{1,-1} mU}], Point[{1,-1} sU],Point[{1,-1} mU],
  Blue, Line[{sPU,mPU}], Point[sPU],Point[mPU],Line[{{1,-1} sPU,{1,-1} mPU}],Point[{1,-1} sPU],Point[{1,-1} mPU]},  Axes -> True,Ticks->False, ImageSize -> 400]],
Text["parameters out of bounds"]],
{{x2, 10,"distance"}, 2, 100},
{{r1, 3,Style["radius",Darker[Yellow]]}, 0, 6},
{{r2, 1,Style["radius",Gray]}, 0, r1},SaveDefinitions->True]

Totality and partiality

• Polygons of totality and partiality for the eclipse of August 21st. The totality polygon is as given by SolarEclipse. The partiality polygon has been recomputed.

totality=Polygon[GeoPosition[]];

partiality=Polygon[GeoPosition[]];

• They have nontrivial shapes.

GeoGraphics[{Red,partiality,Black,GeoStyling[Opacity[0.6]],totality},GeoProjection->"Equirectangular",GeoRange->"World",GeoCenter->-90,GeoZoomLevel->2]

GeoGraphics[{Red,partiality,Black,GeoStyling[Opacity[0.6]],totality},GeoProjection->"Orthographic",GeoRange->"World",GeoCenter->{35,-90},GeoZoomLevel->2,GeoGridLines->Automatic]

• We can understand the geometry if we show the path of the umbra and penumbra cones over the surface of the Earth, as seen from the Sun. In this case the size of the penumbra was roughly double the size of the Moon.

2MoonRadius/EarthSemimajorAxis

Manipulate[
With[{
pos={-9.845+0.00015t,3.08-0.00004t},do=DateObject[{2017,8,21,0,0,t},TimeZone->0]
},
GeoGraphics[{
Text[do,{0.55,1.13}],
Arrow[{{-1.1,0.748},{1.1,0.161}}],
totality,
Disk[pos,0.01],
Opacity[0.5],Disk[pos,0.54],
Red,Point@@partiality
},
GeoProjection->"Orthographic",
GeoRange->"World",
GeoCenter->GIS`AstroGeoPosition["Sun",do],GeoGridLines->Quantity[15,"AngularDegrees"],
PlotRange->{{-1.1,1.1},{-1.,1.2}}]
],
{t,56000,76000},
SaveDefinitions->True]

Orientation of the Earth

• We need the orientation of the Earth in the J2000 frame. This means a rotation matrix from ITRF00 to J2000 coordinates.

• Data obtained from [IERS RS/PC] for the period between 12:00 and 24:00 UTC of August 21st.

angleFunction=InterpolatingFunction[];

vFunction=InterpolatingFunction[];

orientationFunction[secs_]:=RotationMatrix[angleFunction[secs],vFunction[secs]]

• Example of computation of maximum alignment for our location:

geoxyz=GeoPositionXYZ[(Entity["Champaign", Entity["City",{"Champaign","Illinois","UnitedStates"}], "Entity[\"City\", {\"Champaign\", \"Illinois\", \"UnitedStates\"}]", "city"])]

xyz=First[geoxyz]/auInMeters

angle[t_]:=With[{p=earthFunction[t]+orientationFunction[t].xyz},VectorAngle[sunFunction[t]-p,moonFunction[t]-p]/Degree]

ListPlot[{#,angle[#]}&/@Range[20000,26000,60]]

This computation is wrong by a few tens of seconds, because it assumes infinite propagation speed.

Take into account the speed of light

• Compute the center of the umbra at a given UTC time (in seconds since 12:00 GMT).

UmbraCenter[tE_?NumericQ,maxits_:Infinity]:=Module[{tS,tM,sun,earth,moon,vertex,x,y,z,λ,sol},
earth=earthFunction[tE];
tS=tE;
tM=tE;
FixedPoint[(
sun=sunFunction[tS];
moon=moonFunction[tM];
vertex=sun+(moon-sun)(SunRadius/(SunRadius-MoonRadius));
tS=tE-Sqrt[Norm[sun-vertex]^2-SunRadius^2]/LightSpeed;
tM=tE-Sqrt[Norm[moon-vertex]^2-MoonRadius^2]/LightSpeed;
vertex
)&,
0,
maxits
];
{x,y,z}=Inverse[orientationFunction[tE]].(vertex+λ(sun-moon)-earth);
(* The first solution is on the other side of the Earth. Take the second *)
sol=NSolve[((x^2+y^2)/(EarthSemimajorAxis^2))+((z^2)/(EarthSemiminorAxis^2))==1,λ][[2]];
GeoPosition[GeoPositionXYZ[auInMeters{x,y,z}/.sol]]
];

This is the time of maximum eclipse that we reported for Carbondale:

t=QuantityMagnitude[(DateObject["Mon 21 Aug 2017 13:21:24""GMT"-5., DateObject[{2017,8,21,13,21,24},"Instant","Gregorian",-5.]])-(DateObject["Mon 21 Aug 2017 12:00:00""GMT", DateObject[{2017,8,21,12,0,0},"Instant","Gregorian",0.`]]),"Seconds"]

p=UmbraCenter[t]

p0=UmbraCenter[t,1]

GeoDistance[p,p0]

GeoGraphics[{totality,GeoMarker[{(Entity["Carbondale", Entity["City",{"Carbondale","Illinois","UnitedStates"}], "Entity[\"City\", {\"Carbondale\", \"Illinois\", \"UnitedStates\"}]", "city"]),p,p0}]},GeoCenter->(Entity["Carbondale", Entity["City",{"Carbondale","Illinois","UnitedStates"}], "Entity[\"City\", {\"Carbondale\", \"Illinois\", \"UnitedStates\"}]", "city"]),GeoRange->Quantity[40,"Miles"],GeoScaleBar->{"Imperial", "Metric"}]

• Compute many centers:

centers=GeoPosition[UmbraCenter/@Range[17400,28800,300]]

GeoGraphics[{Point[centers],Red,totality},GeoProjection->"Orthographic",GeoGridLines->Automatic]

precisioneclipse.com

• Given a specified location with known latitude and longitude, and the output of the computations previously discussed, generating the results of the precisioneclipse.com website is just a matter of some relatively simple visualization.

View from Champaign

• Champaign was not in the path of totality, but has about 93% coverage. As a result, at it’s peak, the Sun was not totally obscured.

View from Totality

• Carbondale, IL was in the path of totality, nearly right on the center line.

Maps Panel

• First, we need a specified location:

userlocation=GeoPosition[(Entity["Carbondale", Entity["City",{"Carbondale","Illinois","UnitedStates"}], "Entity[\"City\", {\"Carbondale\", \"Illinois\", \"UnitedStates\"}]", "city"])]

• Then, we can use the Wolfram Language functions GeoGraphics and SolarEclipse to generate a map of the United States with the specified user location as a GeoMarker.

GeoGraphics[{GeoMarker[userlocation],Black,SolarEclipse[DateObject[{2017,8,21}],"TotalPhasePolygon"]},GeoRange->(Entity["United States", Entity["Country","UnitedStates"], "Entity[\"Country\", \"UnitedStates\"]", "country"])]

• To get the exact results of the website, it’s just a matter of customizing the GeoMarker.

mypin=Graphics[GraphicsGroup[{FaceForm[Orange],FilledCurve[{{Line[Join[{{0,0}},({Cos[#1],3+Sin[#1]}&)/@Range[-((2 π)/20),π+((2 π)/20),(π/20)],{{0,0}}]]},{Line[(0.5 {Cos[#1],6+Sin[#1]}&)/@Range[0,2 π,(π/20)]]}}]}]]

• With the graphic defined, we can use it as a second argument to GeoMarker and use the “Alignment” option.

GeoGraphics[{GeoMarker[userlocation,mypin,"Alignment"->Bottom],Black,SolarEclipse[DateObject[{2017,8,21}],"TotalPhasePolygon"]},GeoRange->(Entity["United States", Entity["Country","UnitedStates"], "Entity[\"Country\", \"UnitedStates\"]", "country"])]

• The second panel is just a zoomed-in version of the above using a specific GeoCenter.

GeoGraphics[{GeoMarker[userlocation,mypin,"Alignment"->Bottom,"Scale"->Scaled[.1]],Black,SolarEclipse[DateObject[{2017,8,21}],"TotalPhasePolygon"]},GeoRange->Quantity[50,"Miles"],GeoCenter->userlocation]

Contact Points

One of the critical details of eclipse prediction is knowing when the contact points are. During a total solar eclipse, there are 4 contact points denoted C1, C2, C3, and C4. As the Moon starts to make contact with the Sun’s disk, this begins C1. As the eclipse progresses and the Moon’s disk initially covers the Sun’s disk, this denotes C2. Shortly afterwards, just before the Sun’s disk re-emerges, this marks C3. Finally, as the Moon’s disk makes final contact with the Sun’s disk and the eclipse ends, this denotes C4. Only total solar eclipses have C2 and C3. A partial solar eclipse will only have C1 and C4 since the Moon never totally covers the Sun’s disk.

Eclipse Phases - Importance of Projection

• Given the altitude and azimuth of the Sun and Moon at a given time, it sounds straight forward to simply plot 2 disks at those coordinates as follows. The positions chosen are just after first contact so we would expect to see the Moon just starting to cover the Sun.

do=DateObject[{2017,8,21,12,20,0},TimeZone->-5]

moonAltAz=Reverse@QuantityMagnitude[MoonPosition[userlocation,do],"AngularDegrees"]

sunAltAz=Reverse@QuantityMagnitude[SunPosition[userlocation,do],"AngularDegrees"]

sunRadius=Quantity[.2634,"AngularDegrees"];
moonRadius=Quantity[.26763,"AngularDegrees"];

eclipsegraphicSimple[sunAltAz_,moonAltAz_]:=Graphics[{
{Hue[.14,1],Disk[sunAltAz,QuantityMagnitude@sunRadius]},
{Black,Opacity[1],Disk[moonAltAz,QuantityMagnitude@moonRadius]}
},Background->White,PlotRange->All]

eclipsegraphicSimple[sunAltAz,moonAltAz]

• The reason for this discrepancy is because Graphics generates its contents on a flat plane, but the sky is not flat. It’s more like a dome so you must take the projection into account, especially for extended objects like the Sun and Moon. GeoGraphics, although designed originally for Earth-based maps, has support for the necessary projections so we can use it to get the results we need but just turning off the normal background. Since we don’t want to see the entire Lunar disk, we use a white circle as a mask to only see the part that obscures the Sun.

eclipsegraphicWhite[sunAltAz_,moonAltAz_]:=GeoGraphics[{
{GeoStyling[Hue[.14,1]],GeoDisk[sunAltAz,sunRadius]},{GeoStyling[Black,Opacity[1]],GeoDisk[moonAltAz,moonRadius]},{White,Opacity[1],Thickness[.6],GeoCircle[sunAltAz,sunRadius+Quantity[.3,"AngularDegrees"]]}},GeoProjection->{"Orthographic","ReferenceModel"->180/Pi,"Centering"->sunAltAz},GeoBackground->None,Background->White,PlotRange->.5,GeoCenter->sunAltAz]

eclipsegraphicWhite[sunAltAz,moonAltAz]

• For the longer time sequence, we want to show the scene a bit more like what you would see with your eyes. During an eclipse, the sky is dark so we use a dark gray coloring and in this case we don’t obscure the Moon. Otherwise, the scene is nearly identical to the white variant above.

eclipseGraphic[sunAltAz_,moonAltAz_]:=GeoGraphics[{
{GeoStyling[Hue[.14,1]],GeoDisk[sunAltAz,sunRadius]},
{GeoStyling[Black,Opacity[1]],GeoDisk[moonAltAz,moonRadius]}
},GeoProjection->{"Orthographic","ReferenceModel"->180/Pi,"Centering"->sunAltAz},GeoBackground->None,Background->GrayLevel[.1],PlotRange->.5,GeoCenter->sunAltAz]

eclipseGraphic[sunAltAz,moonAltAz]

Related Links

Website: precisioneclipse.com

Wolfram Blog: Get Ready for the Total Solar Eclipse of 2017

Wolfram Blog: Double Eclipse! Or Why Carbondale, Illinois, Is Special

Stephen Wolfram Blog: When Exactly Will the Eclipse Happen? A Multimillenium Tale of Computation

Future

• We are now very interested in eclipses. We will do better for future solar eclipses.

• Continue improving our support for astronomy-related and time-related functionality.

• AstroPosition & AstroGraphics.