The Astro Project

José Martín-García, Wolfram Research

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

Abstract

This talk introduces the new Astro framework for astronomical computation and visualization in the Wolfram Language. Modelled on our successful Geo functionality, AstroPosition and AstroGraphics will allow us to find and visualize many interesting facts about our Solar System and beyond.

Introduction

New in 13.2: Computational Astronomy in the Wolfram Language.

We already have Astro functions (MoonPosition, SiderealTime, ...) and data collections (entities of domains “Star”, “Galaxy”, ...).

Introduce a first-class AstroPosition in WL.

Produce maps of the sky with AstroGraphics, reusing our projections.

Provide fast and precise computation. Take relativity into account.

Use modern data. NASA “kernel” data for the Solar System.

Current Astro Team

Jose Martin-Garcia               Jeff Bryant                           Truman Tapia

Outline

1) Reference frames and time systems

2) AstroPosition

3) Astro computation

4) AstroGraphics

5) Application examples

6) Constellations

Conclusions

1. Reference Frames and Time Systems

A position in space/time is always given wrt something else:

Physics: (reference) frame, time system. Datum in Geo.

Maths: coordinate system, calendar. Projection in Geo.

Time is important in Astro: everything moves. Relativity.

All spherical computation (no ellipsoids).

Precise results only possible using precise frames.

1.1. Main Four Families of Frames

“Equatorial”: use Earth’s equatorial plane.

Variants: “ICRS”, “J2000”, “GCRS”, “TETE”, “MEME”, “CIRS”.

“Ecliptic”: use Earth’s orbital plane.

Variants: “TrueEcliptic”, “MeanEcliptic”.

“Galactic”: use Milky Way’s plane.

“Horizon” frames: use local horizontal plane.

1.2. Time Systems (recap)

Atomic time at the surface of the Earth

Variants: “TAI”, “TT”, “UTC”

Time determined by the rotation of the Earth

Variants: “UT1”

Time at places other than Earth surface.

Variants: “TCG”, “TCB”, “TDB”, ...

dates=DateRange[DateObject[{1960,1,1}],Today,"Month"];
dt=GeoOrientationData[dates,"TAIMinusUT1"];
leap=GeoOrientationData[dates,"TAIMinusUTC"];

Orange: TT - UT1 ± 0.9s,  Blue: TT - UTC

diff=Quantity[32.184,"Seconds"];
Show[DateListPlot[{Transpose[{dates,dt+diff-Quantity[0.9,"Seconds"]}],Transpose[{dates,dt+diff+Quantity[0.9,"Seconds"]}]},PlotStyle->Orange],DateListPlot[Transpose[{dates,leap+diff}]],ImageSize->Large,PlotRange->{32,72}]

2. AstroPosition

2.1. Main ideas

• Follow GeoPosition:

AstroPosition[coords,frame,csys]

coords is a coordinate array.

frame is “Equatorial” or {“Equatorial”, date, loc, ...}, etc.

csys is “Cartesian”, “RADec”, “LonLat”, ... or a projection.

• Example: Observe the Moon’s center from your location:

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

• Angles are stored as degrees, distances as astronomical units:

%//InputForm

• Coordinate conversions for that same date and location:

AstroPosition[%,"Equatorial"]

AstroPosition[%,"Ecliptic"]

AstroPosition[%,"Ecliptic","Cartesian"]

• Dynamic motion of the Moon:

Dynamic[AstroPosition[(Entity["Moon", Entity["PlanetaryMoon","Moon"], "Entity[\"PlanetaryMoon\", \"Moon\"]", "planetary moon"])],UpdateInterval->0.1]

Dynamic[AstroPosition[(Entity["Moon", Entity["PlanetaryMoon","Moon"], "Entity[\"PlanetaryMoon\", \"Moon\"]", "planetary moon"]),"GCRS"],UpdateInterval->0.1]

• Extract properties:

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

%["Altitude"]

DMSString[%,"+-"]

• Directions in space:

AstroPosition["Zenith"]

AstroPosition["NorthCelestialPole"]

The true equinox is the X axis of the TETE and TrueEcliptic frames:

AstroPosition["TrueEquinox","TETE"]

This is a consequence of precession:

AstroPosition["TrueEquinox","J2000"]

This is a consequence of nutation:

AstroPosition["TrueEquinox","MEME"]

2.2. Example: Observe from the Moon

• Observe our location from a point on the Moon:

loc=GeoPosition[(Entity["Mare Tranquillitatis", Entity["SolarSystemFeature","MareTranquillitatisMoon"], "Entity[\"SolarSystemFeature\", \"MareTranquillitatisMoon\"]", "solar system feature"])]

AstroPosition[Here,loc]

AstroPosition[Here,GeoPosition[{0,0},"Moon"]]

Define this function that takes any day d in October 2022:

pos[obj_,d_]:=AstroPosition[obj,{"Horizon",DateObject[{2022,10,d,0,0,0}],GeoPosition[{0,0},"Moon"]}]

Plot[pos[Here,d]["Altitude"],{d,0,60},PlotPoints->100]

Plot[pos["Earth",d]["Altitude"],{d,0,60},PlotPoints->100]

Plot[pos[Here,d]["Altitude"]-pos["Earth",d]["Altitude"],{d,0,60},PlotPoints->100]

2.3. Example: Centers and Barycenters

• Compute the ecliptic Z coordinate of the Earth over a few months:

earthEclipticZ[year_]:=AstroPosition[(Entity["Earth", Entity["Planet","Earth"], "Entity[\"Planet\", \"Earth\"]", "planet"]),{"Ecliptic",DateObject[{year},CalendarType->"GregorianYear"],"SSB","LightDeflection"->False},"Cartesian"]["Data"][[3]]

au=QuantityMagnitude["AstronomicalUnit","Kilometers"]//N

Plot[earthEclipticZ[year]au,{year,2020,2021}]

embEclipticZ[year_]:=AstroPosition[(Entity["Sun", Entity["Star","Sun"], "Entity[\"Star\", \"Sun\"]", "star"]),{"Ecliptic",DateObject[{year},CalendarType->"GregorianYear"],"SSB","LightDeflection"->False},"Cartesian"]["Data"][[3]]

Plot[embEclipticZ[year]au,{year,2020,2021}]

• Compute the position of Jupiter’s barycenter wrt Jupiter’s center:

posJ[d_?NumericQ]:=AstroPosition["JupiterBarycenter",{"Equatorial",DateObject[{2022,10,d},"Instant"],"Jupiter"},"Cartesian"]["Data"][[{1,2}]]au

plot=ParametricPlot[posJ[d],{d,0,100},PlotRange->250]

Manipulate[Show[{plot,Graphics[{Red,PointSize[0.03],Point[posJ[d]]}]}],{d,0,100},]

3. Astro Computation

3.1. AstroDistance

• Compute the current distance to an object from your location:

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

UnitConvert[%,"Kilometers"]

Plot[AstroDistance[(Entity["Moon", Entity["PlanetaryMoon","Moon"], "Entity[\"PlanetaryMoon\", \"Moon\"]", "planetary moon"]),DateObject[{2022,10,d},"Instant"]],{d,0,31},TargetUnits->"Kilometers"]

• Approach to Jupiter last month:

jdist[year_?NumericQ]:=AstroDistance[(Entity["Jupiter", Entity["Planet","Jupiter"], "Entity[\"Planet\", \"Jupiter\"]", "planet"]),DateObject[{year},CalendarType->"GregorianYear"]]

Plot[jdist[year],{year,2010,2030}]

Plot[jdist[year],{year,2022,2023}]

Find the minimum:

FindArgMin[QuantityMagnitude[jdist[year],"Meters"],{year,2022.5,2022,2023}]

DateObject[%,CalendarType->"GregorianYear"]

CalendarConvert[%]

3.2. AstroAngularSeparation

• Angular separation between two objects as seen from a third (your local position by default):

AstroAngularSeparation[(Entity["Venus", Entity["Planet","Venus"], "Entity[\"Planet\", \"Venus\"]", "planet"]),(Entity["Mars", Entity["Planet","Mars"], "Entity[\"Planet\", \"Mars\"]", "planet"])]

• Find the recent date of closest approach between Jupiter and Saturn:

sep[year_?NumericQ]:=AstroAngularSeparation[(Entity["Jupiter", Entity["Planet","Jupiter"], "Entity[\"Planet\", \"Jupiter\"]", "planet"]),(Entity["Saturn", Entity["Planet","Saturn"], "Entity[\"Planet\", \"Saturn\"]", "planet"]),DateObject[{year},CalendarType->"GregorianYear"]]

Plot[sep[year],{year,2000,2040}]

FindArgMin[QuantityMagnitude[sep[year],"Radians"],{year,2021,2020}]

DateObject[%,CalendarType->"GregorianYear"]

CalendarConvert[%]

4. AstroGraphics

4.1. Maps of the whole sky

• By default, show your current local sky. Color code:
	- Green: Horizon plane
	- Purple: Equator plane
	- Blue: Ecliptic plane
	- Brown: Galactic plane

AstroGraphics[]

Manipulate[AstroGraphics[AstroFrame->{"Horizon",DateObject[{2022,10,19,h,0,0}],GeoPosition[{lat,lon}]},ImageSize->Large],{h,0,24},{{lat,40},-90,90},{{lon,-90},-180,180}]

• Use other frames:

AstroGraphics[AstroFrame->{"Equatorial",Now,Here}]

AstroGraphics[AstroFrame->{"Ecliptic",Now,Here}]

AstroGraphics[AstroFrame->{"Galactic",Now,Here}]

• Choose your own projection for any frame. Also add grid lines:

AstroGraphics[AstroProjection->"WinkelTripel",AstroRange->All,AstroFrame->{"Equatorial",Now,Here},AstroGridLines->{"Galactic"->Automatic,"Horizon"->Quantity[15,"AngularDegrees"]}]

• And a background image:

AstroGraphics[AstroFrame->"Equatorial",AstroBackground->"VisibleSky"]

AstroGraphics[AstroFrame->{"Equatorial",Now,Here},AstroBackground->"VisibleSky",AstroProjection->"LambertAzimuthal"]

4.2. Maps of local regions of the sky

• Specify an object and give an angular radius:

AstroGraphics[(Entity["Rigel", Entity["Star","Rigel"], "Entity[\"Star\", \"Rigel\"]", "star"]),AstroRange->Quantity[20,"AngularDegrees"],AstroFrame->"Equatorial"]

Manipulate[AstroGraphics[(Entity["Alnilam", Entity["Star","Alnilam"], "Entity[\"Star\", \"Alnilam\"]", "star"]),AstroRange->Quantity[s,"AngularDegrees"],AstroFrame->"Equatorial",ImageSize->Large],{{s,30},.1,30}]

• Again, we can add the background image:

AstroGraphics[{(Entity["Shaula", Entity["Star","Shaula"], "Entity[\"Star\", \"Shaula\"]", "star"]),White,PointSize[0.02],Point[AstroPosition[{0,0},"Galactic"]]},AstroRange->(QuantityPostfix[10, "°", degrees, "AngularDegrees"]),AstroBackground->"VisibleSky",AstroFrame->"J2000",AstroProjection->"LambertAzimuthal",AstroGridLines->{"Galactic"->Automatic}]

• Show the area of the Andromeda galaxy:

AstroGraphics[(Entity["M31", Entity["Galaxy","M31"], "Entity[\"Galaxy\", \"M31\"]", "galaxy"]),AstroRange->(QuantityPostfix[3, "°", degrees, "AngularDegrees"]),AstroFrame->"Equatorial"]

5. Application Examples

5.1. Zoom in around Jupiter

AstroPosition[(Entity["Jupiter", Entity["Planet","Jupiter"], "Entity[\"Planet\", \"Jupiter\"]", "planet"]),"Equatorial"]

Manipulate[AstroGraphics[(Entity["Jupiter", Entity["Planet","Jupiter"], "Entity[\"Planet\", \"Jupiter\"]", "planet"]),AstroRange->Quantity[2^s,"AngularDegrees"],AstroFrame->{"Equatorial",Now+Quantity[d,"Days"]},ImageSize->Large],{s,3,-3},{d,0,31}]

5.2. Retrograde motion of Mars

• Take a collection of weekly positions of Mars over a year, with a background of today:

dates=DateRange[DateObject[{2022,7,1}],DateObject[{2023,7,1}],7];

mars=Tooltip[Point@AstroPosition[(Entity["Mars", Entity["Planet","Mars"], "Entity[\"Planet\", \"Mars\"]", "planet"]),{"Equatorial",#}],#]&/@dates;

AstroGraphics[{White,mars},AstroFrame->"Equatorial"]

5.3. Solar eclipses

• Eddington's observations of gravitational light deflection:

date=SolarEclipse[DateObject[{1919}]]

location=GeoPosition[(Entity["Principe Island", Entity["Island","PrincipeIsland"], "Entity[\"Island\", \"PrincipeIsland\"]", "island"])]

Visualize what happened during the eclipse, centering around the position of the Sun:

Manipulate[
AstroGraphics[(Entity["Sun", Entity["Star","Sun"], "Entity[\"Star\", \"Sun\"]", "star"]),AstroFrame->{"Equatorial",date+Quantity[t,"Hours"],location},AstroRange->Quantity[1,"AngularDegrees"],ImageSize->Large],{t,-1,6.5},]

Observed position of the star HIP20842, also known as HR1403:

ap[t_,ld_]:=AstroPosition[(Entity["HIP 20842", Entity["Star","HIP20842"], "Entity[\"Star\", \"HIP20842\"]", "star"]),{"Equatorial",date+t,location,"LightDeflection"->ld}]

Compute the separation in arcseconds between the positions corrected and not corrected for gravitational deflection:

sep[t_]:=AstroAngularSeparation[ap[t,True],ap[t,False]]

Plot[sep[Quantity[t,"Hours"]],{t,0,18},TargetUnits->"Arcseconds"]

This is approximately what Eddington measured:

sep[Quantity[1.1,"Hours"]]//UnitConvert[#,"Arcseconds"]&

This is the value for the deflection of a grazing star:

sep[Quantity[6.5,"Hours"]]//UnitConvert[#,"Arcseconds"]&

• Annular eclipse next year. Notice how in this case the Moon looks smaller than the Sun.

Manipulate[AstroGraphics[{(Entity["Sun", Entity["Star","Sun"], "Entity[\"Star\", \"Sun\"]", "star"])},AstroFrame->{"Equatorial","Date"->date,"Location"->(Entity["Eugene", Entity["City",{"Eugene","Oregon","UnitedStates"}], "Entity[\"City\", {\"Eugene\", \"Oregon\", \"UnitedStates\"}]", "city"])},AstroRange->(QuantityPostfix[1, "°", degree, "AngularDegrees"]),PlotLabel->date],{{date,(DateObject["Sat 14 Oct 2023 09:20:00""PDT", DateObject[{2023,10,14,9,20,0},"Instant","Gregorian","America/Los_Angeles"]])},(DateObject["Sat 14 Oct 2023 08:00:00""PDT", DateObject[{2023,10,14,8,0,0},"Instant","Gregorian","America/Los_Angeles"]]),(DateObject["Sat 14 Oct 2023 11:00:00""PDT", DateObject[{2023,10,14,11,0,0},"Instant","Gregorian","America/Los_Angeles"]]),(Quantity[1, "min", minute, "Minutes"])}]

5.4. Aberration

• Stars appear to move in small ellipses due to Earth’s orbit around the Sun.

Manipulate[AstroGraphics[AstroPosition[(Entity["Eltanin", Entity["Star","Etamin"], "Entity[\"Star\", \"Etamin\"]", "star"]),DateObject[{2022,1,d},"Instant"]],AstroRange->Quantity[40,"Arcseconds"],AstroFrame->{"Equatorial",DateObject[{2022,1,1},"Instant"]}],{d,1,365}]

• Take a star of high ecliptic latitude:

star=(Entity["Eltanin", Entity["Star","Etamin"], "Entity[\"Star\", \"Etamin\"]", "star"]);

AstroPosition[star,"Ecliptic"]

Compute values of right ascension and declination at the beginning of 2022. Use the “GCRS” frame, whose axes do not rotate wrt distant stars.

pos[star_,d_,prop_]:=AstroPosition[star,{"GCRS",DateObject[{2022,1,d}]}][prop]

ra=pos[star,1,"RightAscension"]

dec=pos[star,1,"Declination"]

Plot the changes wrt those values for a year, in arcseconds:

ParametricPlot[{
(pos[star,d,"RightAscension"]-ra)Cos[dec Degree],(pos[star,d,"Declination"]-dec)
}3600,{d,1,360}]

• Do the same for a star of lower ecliptic latitude:

star=(Entity["Betelgeuse", Entity["Star","Betelgeuse"], "Entity[\"Star\", \"Betelgeuse\"]", "star"]);

AstroPosition[star,"Ecliptic"]

ra=pos[star,1,"RightAscension"]

dec=pos[star,1,"Declination"]

We get a more eccentric ellipse, but of the same major axis:

ParametricPlot[{(pos[star,d,"RightAscension"]-ra)Cos[dec Degree],pos[star,d,"Declination"]-dec}3600,{d,1,360}]

• The aberration major semiaxis has approximately value v/c :

Quantity[(Entity["Earth", Entity["Planet","Earth"], "Entity[\"Planet\", \"Earth\"]", "planet"])[(EntityProperty["average orbit velocity", EntityProperty["Planet","AverageOrbitVelocity"], "EntityProperty[\"Planet\", \"AverageOrbitVelocity\"]"])]/(Quantity[, "c", speed of light, "SpeedOfLight"]),"Radians"]

UnitConvert[2%,"Arcseconds"]

6. Constellations

• This is the sky around the North Celestial Pole:

AstroGraphics[(Entity["Polaris", Entity["Star","Polaris"], "Entity[\"Star\", \"Polaris\"]", "star"]),AstroRange->Quantity[35,"AngularDegrees"],AstroFrame->"Equatorial"]

• Activate constellation illustrations:

AstroGraphics[(Entity["Polaris", Entity["Star","Polaris"], "Entity[\"Star\", \"Polaris\"]", "star"]),AstroRange->Quantity[35,"AngularDegrees"],AstroFrame->"Equatorial",AstroBackground->AstroStyling[{"DarkSky","ShowConstellationIllustrations"->True,"ShowMainPlanes"->False}]]

AstroGraphics[(Entity["Rigel", Entity["Star","Rigel"], "Entity[\"Star\", \"Rigel\"]", "star"]),AstroRange->Quantity[35,"AngularDegrees"],AstroFrame->"Equatorial",AstroBackground->AstroStyling[{"DarkSky","ShowConstellationIllustrations"->True,"ShowMainPlanes"->False}],AstroProjection->"LambertAzimuthal"]

• Use styling options for full control:

style={
"ShowConstellationIllustrations"->True,
"ShowConstellationRegions"->True,"ConstellationRegionStyle"->Darker[Blue,.2],"ConstellationLineStyle"->Opacity[0.6,LightBlue],
"ShowConstellationBoundaries"->True,"ConstellationBoundaryStyle"->Directive[Dashing[{0,Tiny}],Yellow],"ConstellationNameStyle"->GrayLevel[0.7]
};

AstroGraphics[(Entity["Scorpius", Entity["Constellation","Scorpius"], "Entity[\"Constellation\", \"Scorpius\"]", "constellation"]),AstroBackground->AstroStyling[{"VisibleSky","ShowConstellations"->{(Entity["Scorpius", Entity["Constellation","Scorpius"], "Entity[\"Constellation\", \"Scorpius\"]", "constellation"])},style}],AstroFrame->"J2000"]

AstroGraphics[{(Entity["Ophiuchus", Entity["Constellation","Ophiuchus"], "Entity[\"Constellation\", \"Ophiuchus\"]", "constellation"]),(Entity["Serpens", Entity["Constellation","Serpens"], "Entity[\"Constellation\", \"Serpens\"]", "constellation"])},AstroBackground->AstroStyling[{"VisibleSky","ShowConstellations"->{(Entity["Ophiuchus", Entity["Constellation","Ophiuchus"], "Entity[\"Constellation\", \"Ophiuchus\"]", "constellation"]),(Entity["Serpens", Entity["Constellation","Serpens"], "Entity[\"Constellation\", \"Serpens\"]", "constellation"])},style}],AstroFrame->"J2000"]

AstroGraphics[(Entity["Aries", Entity["Constellation","Aries"], "Entity[\"Constellation\", \"Aries\"]", "constellation"]),AstroBackground->AstroStyling[{"VisibleSky","ShowConstellations"->{(Entity["Aries", Entity["Constellation","Aries"], "Entity[\"Constellation\", \"Aries\"]", "constellation"])},style}],AstroFrame->"J2000"]

Conclusions

New Astro project in 13.2. Experimental status

AstroPosition describes observations of anything from anywhere any time.

AstroGraphics shows maps of the sky in several frames and projections.

There is still much to do:

Velocities, orbits and orbital elements.

All about eclipses, occultations, transits and other events.

More data, more types of charts, more computational functions, ...

Thank you!

PacletDirectoryLoad["~/Devel/Repositories/AstronomyConvenienceFunctions"];
$Line=0;