Calculate positions of astronomical objects from a terrestrial observatory

I am currently converting some AstroPy-based Python libraries to Julia and one of them (which is the most important one) is used to calculate the positions of astronomical objects for a given (detector) location and observatory site. I cannot wrap AstroPy because the main purpose of the Julia-port is to make it fast (AstroPy is very slow, at least in this context).

That’s quite an essential calculation and essentially what every planetarium software can do, so I assume that there is at least one package for Julia which can do this but I could not find one yet (I looked e.g. in Julia Astro · GitHub)

Anyways, long story short, what I am currently trying to do is:

UTM coordinates of an observation site -> AltAz (height) coordinates -> sky coordinates of an astronomical object

The first part is easy, using Geodesy.jl:

julia> using Geodesy

julia> utmz = UTMZ(268221.6, 4742381.9, -2500.0, 32, true)
UTMZ(268221.6, 4.7423819e6, -2500.0, zone=32 (north))

julia> LLA(utmz, wgs84)
LLA(lat=42.79891663180618°, lon=6.1657005045321345°, alt=-2500.0)

The second part is to get the AltAz coordinates of e.g. the Sun for a given time and location.

Here is an example of the AstroPy-based code which shows how to get from a time 2015-02-17T08:24:59 and location lat=42.79891663180618°, lon=6.1657005045321345°, alt=-2500.0 the position of the Sun:

>>> import astropy

>>> from astropy.units import deg

>>> from astropy.coordinates import get_body

>>> from km3astro.time import np_to_astrotime

>>> from astropy.coordinates import EarthLocation, Longitude, Latitude, AltAz

>>> d = astropy.time.Time(["2015-02-17T08:24:59"])

>>> sun = get_body('sun', d)

>>> sun
<SkyCoord (GCRS: obstime=['2015-02-17T08:24:59.000'], obsgeoloc=(0., 0., 0.) m, obsgeovel=(0., 0., 0.) m / s): (ra, dec, distance) in (deg, deg, AU)
    [(330.29789687, -12.12304247, 0.98810396)]>

>>> loc = EarthLocation.from_geodetic(
...         lon=Longitude(6.1657005045321345 * deg),
...         lat=Latitude(42.79891663180618 * deg),
...         height=-2500,
...     )

>>> frame = AltAz(obstime=d, location=loc)

>>> sun.transform_to(frame)
<SkyCoord (AltAz: obstime=['2015-02-17T08:24:59.000'], location=(4658222.52798148, 503223.58589033, 4309439.14018905) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt, distance) in (deg, deg, AU)
    [(126.85810735, 17.99016698, 0.98809073)]>

Any ideas which package(s) can do this?

I managed to get the sun position using AstroLib.jl from @giordano but it gives a slightly different initial position (I guess related to some default observation site?) since AstroPy gives [(330.29789687, -12.12304247, 0.98810396)]:

julia> using AstroLib, Dates

julia> d = DateTime("2015-02-17T08:24:59.000")
2015-02-17T08:24:59

julia> sunpos(jdcnv(d))
(330.5004065151536, -12.048606175173527, 328.3455596527475, 23.434817151869996)

and the observatory:

julia> Observatory("ORCA", 42.79891663180618, 6.1657005045321345, -2500, 0)
Observatory: ORCA
latitude:    42.79891663180618°N
longitude:   6.1657005045321345°E
altitude:    -2500.0 m
time zone:   UTC+0

So now I guess I have to utilise SkyCoords.jl to put the puzzle pieces together? :wink:

2 Likes

Here’s one way to do it using GitHub - barrettp/Astrometry: Astrometry is a set of IAU standard algorithms for calculating the time and position of celestial objects. (using ERFA.jl is basically the same):

julia> using Dates

julia> using Astrometry.SOFA: utctai, taitt, epv00, c2s, atco13

julia> lon=deg2rad(6.1657005045321345)
0.10761177449596131

julia> lat=deg2rad(42.79891663180618)
0.7469820115115795

julia> ht=-2500.0
-2500.0

julia> dt=DateTime(2015,2,17,8,24,59)
2015-02-17T08:24:59

julia> jd=datetime2julian(dt), 0.0
(2.4570708506828705e6, 0.0)

julia> tt=taitt(utctai(jd...)...)
(day = 2.4570708506828705e6, fraction = 0.0007775925925925925)

julia> (ph,vh),(pb,vb)=epv00(tt...)
(helio = [[-0.8391921145803088, 0.4786001678353725, 0.20747989429881425], [-0.009368464035175597, -0.013470109416865897, -0.00583911256445573]], bary = [[-0.8361701400394311, 0.4781065457088137, 0.20711139490691688], [-0.009364810692891077, -0.013464963923163561, -0.00583698267640598]])

julia> rc, dc = c2s(-1 .* ph)
(-0.5183029384180501, -0.21155225242615236)

julia> rad2deg(mod2pi(rc)), rad2deg(dc) # Sun RA/dec
(330.3034291194167, -12.12105121050476)

julia> sunobs = atco13(rc, dc, 0.0, 0.0, 0.0, 0.0, jd..., 0.0, lon, lat, ht, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
(azi = 2.214123455663249, zen = 1.2567460055547417, ha = -0.8916527461637037, dec = -0.2102781487602873, ra = 5.76493468419683, eo = -0.0034065630114019875)

julia> rad2deg(sunobs.azi), 90-rad2deg(sunobs.zen) # Sun az/el
(126.85992933042543, 17.99375796178859)
1 Like

You could use EarthOrientation.jl to get values for the dut1 and px/py (polar motion) parameters to atco13 for improved accuracy. I’m assuming neutrinos don’t refract much… :smile:

1 Like

Awesome @DaveMacMahon many thanks for your expertise :smiley: I need to digest all those abbreviations tomorrow :wink: I know TAI at least :laughing:

Years ago I started to try to write a package to output star charts where the user could specify any (lat, lon) + datetime. I may have the functions you need here StarCharts.jl/src/coordinates.jl, but don’t look too closely at the other code because its awful. I copied the equations from an astro textbook and do not guarantee their accuracy. The time.jl file will have the functions for calculations local hour angle and local standard time, but it may be better to use AstroTime.jl as it may contain more accurate functions.

1 Like

I’m assuming neutrinos don’t refract much
As an ex high energy physicist I assume the same thing… then again assumption is the mother of all…
A search reveals nuindex.pdf (bnl.gov)

Hmmm. that paper is prime and ready for translation into Julia
[Well you know what to do then - Ed]

1 Like

Very interesting indeed :wink:

Maybe I should have specified “atmospheric refraction”… :stuck_out_tongue_winking_eye:

1 Like

There’s also GitHub - JuliaAstro/SPICE.jl: Julia wrapper for NASA NAIF's SPICE toolkit that you can use for such problems.