# Calculate positions of astronomical objects from a terrestrial observatory

**URL:** <https://discourse.julialang.org/t/calculate-positions-of-astronomical-objects-from-a-terrestrial-observatory/107772>\
**Category:** Astro/Space\
**Tags:** question\
**Created:** [December 18, 2023, 2:17pm UTC](https://discourse.julialang.org/t/calculate-positions-of-astronomical-objects-from-a-terrestrial-observatory/107772 "2023-12-18T14:17:49Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![tamasgal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamasgal/32/27946_2.png) [@tamasgal](https://discourse.julialang.org/u/tamasgal)\
**Post date:** [December 18, 2023, 2:17pm UTC](https://discourse.julialang.org/t/calculate-positions-of-astronomical-objects-from-a-terrestrial-observatory/107772/1 "2023-12-18T14:17:49Z")

</div>

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](https://github.com/orgs/JuliaAstro/repositories?type=all))

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
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:

```python
>>> 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
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
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? 😉

---

<div class="post-metadata">

**Author:** ![DaveMacMahon](https://avatars.discourse-cdn.com/v4/letter/d/9f8e36/32.png) [@DaveMacMahon](https://discourse.julialang.org/u/DaveMacMahon)\
**Post date:** [December 19, 2023, 1:00am UTC](https://discourse.julialang.org/t/calculate-positions-of-astronomical-objects-from-a-terrestrial-observatory/107772/2 "2023-12-19T01:00:43Z")

</div>

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.](https://github.com/barrettp/Astrometry) (using ERFA.jl is basically the same):

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![DaveMacMahon](https://avatars.discourse-cdn.com/v4/letter/d/9f8e36/32.png) [@DaveMacMahon](https://discourse.julialang.org/u/DaveMacMahon)\
**Post date:** [December 19, 2023, 1:04am UTC](https://discourse.julialang.org/t/calculate-positions-of-astronomical-objects-from-a-terrestrial-observatory/107772/3 "2023-12-19T01:04:38Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![tamasgal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamasgal/32/27946_2.png) [@tamasgal](https://discourse.julialang.org/u/tamasgal)\
**Post date:** [December 19, 2023, 1:11am UTC](https://discourse.julialang.org/t/calculate-positions-of-astronomical-objects-from-a-terrestrial-observatory/107772/4 "2023-12-19T01:11:21Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mihalybaci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mihalybaci/32/13528_2.png) [@mihalybaci](https://discourse.julialang.org/u/mihalybaci)\
**Post date:** [December 19, 2023, 12:48pm UTC](https://discourse.julialang.org/t/calculate-positions-of-astronomical-objects-from-a-terrestrial-observatory/107772/5 "2023-12-19T12:48:44Z")

</div>

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](https://github.com/mihalybaci/StarCharts.jl/blob/master/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](https://github.com/mihalybaci/StarCharts.jl/blob/master/src/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](https://github.com/JuliaAstro/AstroTime.jl) as it may contain more accurate functions.

---

<div class="post-metadata">

**Author:** ![johnh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnh/32/3615_2.png) [@johnh](https://discourse.julialang.org/u/johnh)\
**Post date:** [December 19, 2023, 12:55pm UTC](https://discourse.julialang.org/t/calculate-positions-of-astronomical-objects-from-a-terrestrial-observatory/107772/6 "2023-12-19T12:55:01Z")

</div>

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)](https://www.phy.bnl.gov/~diwan/nuindex.pdf)

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

---

<div class="post-metadata">

**Author:** ![tamasgal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamasgal/32/27946_2.png) [@tamasgal](https://discourse.julialang.org/u/tamasgal)\
**Post date:** [December 19, 2023, 1:20pm UTC](https://discourse.julialang.org/t/calculate-positions-of-astronomical-objects-from-a-terrestrial-observatory/107772/7 "2023-12-19T13:20:25Z")

</div>

Very interesting indeed 😉

---

<div class="post-metadata">

**Author:** ![DaveMacMahon](https://avatars.discourse-cdn.com/v4/letter/d/9f8e36/32.png) [@DaveMacMahon](https://discourse.julialang.org/u/DaveMacMahon)\
**Post date:** [December 20, 2023, 6:54am UTC](https://discourse.julialang.org/t/calculate-positions-of-astronomical-objects-from-a-terrestrial-observatory/107772/8 "2023-12-20T06:54:50Z")

</div>

> [@johnh](#):
>
> 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)](https://www.phy.bnl.gov/~diwan/nuindex.pdf)

Maybe I should have specified “atmospheric refraction”… 😜

---

<div class="post-metadata">

**Author:** ![fgerick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fgerick/32/13228_2.png) [@fgerick](https://discourse.julialang.org/u/fgerick)\
**Post date:** [December 20, 2023, 9:12am UTC](https://discourse.julialang.org/t/calculate-positions-of-astronomical-objects-from-a-terrestrial-observatory/107772/9 "2023-12-20T09:12:34Z")

</div>

There’s also [GitHub - JuliaAstro/SPICE.jl: Julia wrapper for NASA NAIF's SPICE toolkit](https://github.com/JuliaAstro/SPICE.jl) that you can use for such problems.