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?