[AstroPy] ICRS, CIRS, AltAz, difference 16 arcmin

[Markus Wildi]

Hello

I'm quite new to AstroPy and I'm using AstroPy 1.2.1,
NumPy 1.11.1, Python 3.5 on Ubuntu 16.04.

I wrote a script that calculates the AltAz() coordinates
from SkyCoord(). To double check the result I calculated
AltAz coordinates using a simple matrix transforming
CIRS eq. coordinates to AltAz.

I checked

1) the value of the sidereal time (good within less
than 1 sec).

2) the matrix transformation using ha=dec=0. with no
intermediate apparent place calculation. This yields
Az=180 deg and the co-latitude of the observatory.

While the Alt coordinates in this example are in good
agreement the Az show a difference of 16 arcmin although
I set the relevant parameters to 0. (should not influence
Az anyway):

ap_aa=astr_eq.transform_to(AltAz(
  obstime=dtm,
  location=obs,
  obswl=0.5*u.micron,
  pressure=0.*u.hPa,
  temperature= 0.*u.deg_C,
  relative_humidity=0.))

Using different ha values, e.g. ha=15 deg, Alt difference
is about 3 arcmin.

I can not explain this difference and since AstroPy
has been existing for a long time I got the impression
I missed something (in the documentation).

I appreciate any help.

Kind regards, Markus Wildi

The script at the bottom produces the following
(commented) output:

obs    lon: 7d34m59.88s
obs    lat: 47d47m27.96s
obs co-lat: 42d12m32.04s
sid       : 2h03m23.7265s

Apparent position for ha=dec=0.,
ra=Longitude(sid.radian-ha,u.radian):

CIRS   ra : 2h03m23.1977s # ok, sid
CIRS   dec: 0d04m22.9374s

Transformation to AltAz from ICRS (astropy):

AltAz  az : 179d44m02.1999s
AltAz  alt: 42d16m53.7929s # see result rotmat

Transformation to AltAz from CIRS (rotmat):

rotmat az : 180d00m10.7215s
rotmat alt: 42d16m54.9773s # see astropy

Direct transformation of ha=dec=0. (no
intermediate apparent place), seems to
be ok:

astrometic (ha=dec=0., cross check for rotmat)
rotmat az : 180d00m00.0009s
rotmat alt: 42d12m32.04s



#!/usr/bin/env python
import numpy as np
from astropy import units as u
from astropy.time import Time
from astropy.coordinates import SkyCoord,EarthLocation
from astropy.coordinates import AltAz,CIRS
from astropy.coordinates import Angle,Longitude,Latitude

import numpy as np
from numpy import matrix

def y_rot(rads=None):
  return np.matrix([
  [ np.cos(rads), 0., np.sin(rads)],
  [ 0.,           1., 0.],
  [-np.sin(rads), 0., np.cos(rads)]])

def __xyz(lat=None,lon=None):
    s=np.array([np.cos(lat)*np.cos(lon),
                np.cos(lat)*np.sin(lon),
                np.sin(lat)]).transpose()
    s.shape=(3,1)
    return s

def EQtoAltAz(ra=None,dec=None,sid=None,lat=None):
    ha=sid.radian-ra
    mt=y_rot(lat-np.pi/2.)
    t = mt * __xyz(lat=dec,lon=-ha)
    tlat=np.arcsin(float(t[2]))
    tlon=-np.arctan2(float(t[1]),float(t[0]))
    return tlon,tlat


obs=EarthLocation(lon=7.5833*u.deg,lat=+47.7911*u.deg,height=280.*u.m)
print('obs    lon: {}'.format(obs.longitude.to_string(unit=u.deg)))
print('obs    lat: {}'.format(obs.latitude.to_string(unit=u.deg)))

cl=Latitude(np.pi/2. - obs.latitude.radian, u.radian)

print('obs co-lat: {}'.format(cl.to_string(unit=u.deg)))

dtm=Time('2015-06-30T07:00:56.05',format='isot', scale='utc',location=obs)
sid=dtm.sidereal_time(kind='apparent')

ha=dec=0.
ra=Longitude(sid.radian-ha,u.radian)
radec_str='{} {}'.format(ra.to_string(unit=u.hour),dec)
astr_eq=SkyCoord(radec_str, unit=(u.hourangle,u.deg ))

print('sid       :{}'.format(
                     Angle(sid.radian,u.radian).to_string(u.hour)))

ap_eq=astr_eq.transform_to(CIRS(obstime=dtm))

print('CIRS   ra : {}'.format(ap_eq.ra.to_string(unit=u.hour)))
print('CIRS   dec: {}'.format(ap_eq.dec.to_string(unit=u.deg)))

ap_aa=astr_eq.transform_to(AltAz(obstime=dtm,
                                 location=obs,
                                 obswl=0.5*u.micron,
                                 pressure=0.*u.hPa,
                                 temperature= 0.*u.deg_C,
                                 relative_humidity=0.))

print('AltAz  az : {}'.format(ap_aa.az.to_string(unit=u.deg)))
print('AltAz  alt: {}'.format(ap_aa.alt.to_string(unit=u.deg)))

az,alt=EQtoAltAz(ra=ap_eq.ra.radian,
                 dec=ap_eq.dec.radian,
                 sid=sid,
                 lat=obs.latitude.radian)
aaz=Longitude(az+np.pi, u.radian) # North Az=0., E=np.pi/2.
aalt=Latitude(alt, u.radian)

print('rotmat az : {}'.format(aaz.to_string(unit=u.degree)))
print('rotmat alt: {}'.format(aalt.to_string(unit=u.degree)))

az,alt=EQtoAltAz(ra=astr_eq.ra.radian,
                 dec=astr_eq.dec.radian,
                 sid=sid,
                 lat=obs.latitude.radian)
aaz=Longitude(az+np.pi, u.radian)
aalt=Latitude(alt, u.radian)

print('astrometic (ha=dec=0., cross check for rotmat)')
print('rotmat az : {}'.format(aaz.to_string(unit=u.degree)))
print('rotmat alt: {}'.format(aalt.to_string(unit=u.degree)))