Mar-21-2019, 08:14 PM
There is no error(inconsistent use of tabs and space) in code you posted.
As i did install astropy as test in your other thread i can run code.
Code is from a Jupyter Notebook,so running as script need to add
As i did install astropy as test in your other thread i can run code.
Code is from a Jupyter Notebook,so running as script need to add
print()
to get output.# ---------bryna.py ------------ import numpy as np from astropy.coordinates import SkyCoord, Angle, EarthLocation from astropy.time import Time from astropy import units # mjd = 55780.1 latitude = Angle('-26d42m11.94986s') longitude = Angle('116d40m14.93485s') # obs_time = Time(mjd, format='mjd', location = (longitude, latitude)) lst_apparent = obs_time.sidereal_time('apparent') print(lst_apparent) # Output1: 7h11m46.2716s # # The frame radio astronomers call the apparent or current epoch is the # "true equator & equinox" frame, notated E_gamma in the USNO circular 179 # (http://aa.usno.navy.mil/publications/docs/Circular_179.php) # astropy doesn't have this frame but it's pretty easy to adapt the CIRS frame # by modifying the ra to reflect the difference between # GAST (Grenwich Apparent Sidereal Time) and the earth rotation angle (theta). # The earth rotation angle is most clearly explained in the 2nd paragraph of Chapter 6 # of the USNO circular (page 50) # The relationship between E_gamma and CIRS is shown in the figure on page 57 of the circular def egamma_to_cirs_ra(egamma_ra, time): from astropy import _erfa as erfa from astropy.coordinates.builtin_frames.utils import get_jd12 era = erfa.era00(*get_jd12(Time(mjd, format='mjd'), 'ut1')) theta_earth = Angle(era, unit='rad') assert(isinstance(time, Time)) gast = time.sidereal_time('apparent', longitude=0) cirs_ra = egamma_ra - (gast-theta_earth) return cirs_ra # If we define a source at RA = apparent LST and Dec = latitude it should be at zenith loc_obj = EarthLocation.from_geodetic(lon=longitude, lat=latitude) # use CIRS but with a different ra to account for the difference between LST & earth's rotation angle cirs_ra = egamma_to_cirs_ra(lst_apparent, Time(mjd, format='mjd')) egamma_zenith_coord = SkyCoord(ra=cirs_ra, dec=latitude, frame='cirs', obstime=Time(mjd, format='mjd'), location = loc_obj) print('-' * 20) print(egamma_zenith_coord) # Output2 #<SkyCoord (CIRS: obstime=55780.1): (ra, dec) in deg # (107.78961213, -26.70331941)> # # check where it is in altaz (should be near zenith): egamma_zenith_altaz = egamma_zenith_coord.transform_to('altaz') print('-' * 20) print(egamma_zenith_altaz) # Output3 #<SkyCoord (AltAz: obstime=55780.1, location=(-2559302.5737783727, 5095070.526830904, -2848887.400942108) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0, obswl=1.0 micron): (az, alt) in deg # (28.60732437, 89.99985797)> # n = (egamma_zenith_altaz.alt - Angle('90d')).to_string(unit=units.degree, sep=('deg', 'm', 's')) print('-' * 20) print(n) # Output4 # '-0deg00m00.5113s' # # ------ EOF: bryna.py ----------Test:
Output:(del_env) E:\div_code\del_env
λ python as.py
7h11m46.2716s
--------------------
<SkyCoord (CIRS: obstime=55780.1): (ra, dec) in deg
(107.78961179, -26.70331941)>
--------------------
<SkyCoord (AltAz: obstime=55780.1, location=(-2559302.57377837, 5095070.5268309, -2848887.40094211) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt) in deg
(28.60357927, 89.99985796)>
--------------------
-0deg00m00.5114s