# 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?

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

0.10761177449596131

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)

(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)

(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…

1 Like

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

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

Maybe I should have specified “atmospheric refraction”…

1 Like

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