skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.4k stars 211 forks source link

Apparent RA and DEC #694

Closed kpachar closed 2 years ago

kpachar commented 2 years ago

Hello,

I am able to find the apparent position of a body expressed as altitude and azimuth while accounting for refraction using the local atmospheric temperature and pressure, like this:

def find_altaz(body, t):
    earth = planets['earth']
    observer = earth + wgs84.latlon(12.971599 * N, 77.594566 * E, elevation_m=914)
    alt, az, distance = observer.at(t)\
         .observe(body)\
         .apparent()\
         .altaz(temperature_C=27.0, pressure_mbar=913)
    return {'alt': alt, 'az': az}

However, I can't pass the temperature and pressure to the method that calculates RA and DEC, as the radec method does not accept those arguments.

def find_radec(body, t):
    earth = planets['earth']
    observer = earth + wgs84.latlon(12.971599 * N, 77.594566 * E, elevation_m=914)
    # print(t.utc)
    ra, decl, distance = observer.at(t)\
         .observe(body)\
         .apparent()\
         .radec('date') # <===== Can't pass temperature and pressure_mbar
    return {'ra': ra, 'dec': decl}

Is there any way to adjust the returned right ascention and declenation values to account for the atmospheric temperature and pressure like the way it is done for the azimuth and altitude?

Thank you!

brandon-rhodes commented 2 years ago

That's a good question — so that I can think about where in the documentation this might belong, and how Skyfield might in the future make a prettier API for this, could you share whether you're aiming a visual telescope or doing something like aiming a radio antenna?

The way to approach the problem with Skyfield's current routines would be to ask for the refracted Alt/Az, turn them back into rectangular coordinates, and rotate them back into the main ICRS (J2000) reference frame. It takes a few steps, and requires dropping into lower-level Skyfield code:

import numpy as np
from skyfield.functions import from_spherical
from skyfield.api import load, wgs84, N, E, Distance, Velocity
from skyfield.positionlib import ICRS

ts = load.timescale()
t = ts.utc(2020, 4, 22, 19)
eph = load('de421.bsp')
earth = eph['earth']
body = eph['Jupiter Barycenter']

place = wgs84.latlon(12.971599 * N, 77.594566 * E, elevation_m=914)
observer = earth + place
apparent = observer.at(t).observe(body).apparent()

alt, az, distance = apparent.altaz()
print('Unrefracted Alt:', alt, '  Az:', az)
alt, az, distance = apparent.altaz(temperature_C=27.0, pressure_mbar=913)
print('  Refracted Alt:', alt, '  Az:', az)
ra, dec, distance = apparent.radec('date')
print('Unrefracted  RA:', ra, '  Dec:', dec)

r = from_spherical(distance.au, alt.radians, az.radians)
d = Distance(r)
v = Velocity(r * np.nan)
p = ICRS.from_time_and_frame_vectors(t, place, d, v)

ra, dec, distance = p.radec('date')
print('  Refracted  RA:', ra, '  Dec:', dec)

This is what I get out:

Unrefracted Alt: 00deg 15' 00.5"   Az: 111deg 35' 56.4"
  Refracted Alt: 00deg 38' 21.6"   Az: 111deg 35' 56.4"
Unrefracted  RA: 19h 54m 10.80s   Dec: -20deg 57' 39.3"
  Refracted  RA: 19h 52m 33.81s   Dec: -20deg 51' 58.5"

(I deliberately chose a body and time when it would be close to the horizon, so the Alt difference would be big enough to matter.)

First, take a look at these numbers and confirm that they are the ones you want! I could of course have accidentally computed the wrong thing here without meaning to.

Second, see if some other bodies and times put into the above code result in numbers that look good to you as well.

Then we can think about whether a future version of Skyfield could make things more convenient.

kpachar commented 2 years ago

Thanks a lot. I'm new to Astronomy and it has been hardly a few weeks since I first stumbled upon terms like Right Ascention etc, and was just curious (I don't have any telescopes), just wanted to figure out where the Moon or Sun may be in those coordinates.

brandon-rhodes commented 2 years ago

In the hope that my calculation above answers the question, I'm going to go ahead and close this issue. If you do get the chance to compare the result to observation, then feel free to comment further here if there is any discrepancy!