"""A set of routines to assist in the WCS transforms used in the ``cube_build`` step."""
import math
import numpy as np
__all__ = ["radec2std", "std2radec", "v2v32radec_estimate", "radec2v2v3_estimate"]
# _______________________________________________________________________
[docs]
def radec2std(crval1, crval2, ra, dec, rot_angle=None):
"""
Compute the tangent projection coordinates ``(xi,eta)`` from ``ra,dec`` using ``crval1`` and ``crval2``.
Parameters
----------
crval1 : float
RA value of tangent point
crval2 : float
Dec value of tangent point
ra : ndarray or float
A list (or single value) of RA points to convert
dec : ndarray or float
A list (or single value) of Dec points to convert
rot_angle : float or None, optional
Rotation angle given in degrees
Returns
-------
xi : float
X-axis tangent plane coordinate of ``ra,dec``
eta : float
Y-axis tangent plane coordinate of ``ra,dec``
""" # noqa: E501
if np.isscalar(ra):
ra = np.asarray([ra])
dec = np.asarray([dec])
rad2arcsec = (180.0 * 3600.0) / math.pi
deg2rad = math.pi / 180.0
ra0 = crval1 * deg2rad
dec0 = crval2 * deg2rad
radiff = ra * deg2rad - ra0
decr = dec * deg2rad
h = np.sin(decr) * math.sin(dec0) + np.cos(decr) * math.cos(dec0) * np.cos(radiff)
xi = np.cos(decr) * np.sin(radiff) / h
eta = (np.sin(decr) * math.cos(dec0) - np.cos(decr) * math.sin(dec0) * np.cos(radiff)) / h
xi = -xi
xi = xi * rad2arcsec
eta = eta * rad2arcsec
# xi is made negative so it increases in the opposite direction
# of ra to match the images the Parity of the ifu_cube.
if rot_angle is not None:
temp1 = xi * np.cos(-rot_angle * deg2rad) - eta * np.sin(-rot_angle * deg2rad)
temp2 = xi * np.sin(-rot_angle * deg2rad) + eta * np.cos(-rot_angle * deg2rad)
xi = temp1
eta = temp2
return xi, eta
# ________________________________________________________________________________
[docs]
def std2radec(crval1, crval2, xi, eta):
"""
Compute ``ra,dec`` from the tangent plane rectangular coordinates.
Compute the ``ra,dec`` values of tangent plane rectangular coordinates using
``crval1, crval2`` (the tangent point). This routine takes the rectangular
plane and projects it onto the spherical plane using ``crval1, crval2`` as
the tangent plane.
Parameters
----------
crval1 : float
RA value of tangent point
crval2 : float
Dec value of tangent point
xi : float
Xi rectangular coordinate of tangent plane projected ``ra,dec``
eta : float
Eta rectangular coordinate of tangent plane projected ``ra,dec``
Returns
-------
ra : float
List (or single value) of RA computed values
dec : float
List (or single value) of Dec computed values
"""
if np.isscalar(xi):
eta = np.asarray([eta])
xi = np.asarray([xi])
rad2arcsec = (180.0 * 3600.0) / math.pi
deg2rad = math.pi / 180.0
ra0 = crval1 * deg2rad
dec0 = crval2 * deg2rad
# tangent projection
xi = xi / rad2arcsec
eta = eta / rad2arcsec
xi = -xi # xi is made negative so it increases in the opposite direction
# of ra to match the images the Parity of the ifu_cube is
# for ra is PC1_1 = -1
ra0 = crval1 * deg2rad
dec0 = crval2 * deg2rad
beta = math.cos(dec0) - eta * math.sin(dec0)
angle = xi / beta
ra = np.arctan(angle) + ra0
gamma = np.sqrt(xi * xi + beta * beta)
angle = eta * math.cos(dec0) + math.sin(dec0)
angle = angle / gamma
dec = np.arctan(angle)
ra = ra / deg2rad
dec = dec / deg2rad
mask = ra < 0
ra[mask] += 360.0
mask = ra > 360.0
ra[mask] -= 360.0
return ra, dec
# _______________________________________________________________________
[docs]
def v2v32radec_estimate(ra_ref, dec_ref, roll_ref, v2_ref, v3_ref, v2, v3):
"""
Estimation of RA and Dec from the ``v2, v3`` coordinates.
This routine is used for debugging purposes. It is not actually used
in the ``cube_build`` step for routine IFU cube building.
The conversion from V2,V3 to ``ra,dec`` is handled more accurately by
the transforms provided by :ref:`assign_wcs <assign_wcs_step>`.
Parameters
----------
ra_ref : float
RA of reference point given in arcseconds
dec_ref : float
Dec of reference point given in arcseconds
roll_ref : float
Roll angle given in degrees
v2_ref : float
V2 coordinate of reference point given in arcseconds
v3_ref : float
V3 coordinate of reference point given in arcseconds
v2 : float
V2 coordinate given in arcseconds
v3 : float
V3 coordinate given in arcseconds
Returns
-------
ra : float
Calculated RA from ``v2, v3``
dec : float
Calculate Dec from ``v2, v3``
Notes
-----
It is assumed that the ``v2,v3`` coordinates have the effects of dithering included.
"""
d2r = math.pi / 180.0
v2deg = v2.copy() / 3600.0 # convert to degrees
v3deg = v3.copy() / 3600.0 # convert to degrees
v2_ref = v2_ref / 3600.0 # convert to degrees
v3_ref = v3_ref / 3600.0 # convert to degrees
v3_ref_rad = v3_ref * d2r
roll_ref_rad = roll_ref * d2r
delta_v2 = (v2deg - v2_ref) * math.cos(v3_ref_rad)
delta_v3 = v3deg - v3_ref
delta_ra = delta_v2 * math.cos(roll_ref_rad) + delta_v3 * math.sin(roll_ref_rad)
delta_dec = -delta_v2 * math.sin(roll_ref_rad) + delta_v3 * math.cos(roll_ref_rad)
ra = ra_ref + delta_ra / math.cos(dec_ref * d2r)
dec = delta_dec + dec_ref
return ra, dec
# _______________________________________________________________________
[docs]
def radec2v2v3_estimate(ra_ref, dec_ref, roll_ref, v2_ref, v3_ref, ra, dec):
"""
Convert ``ra,dec`` to ``v2, v3``.
This routine is used for debugging purposes. It is not actually used
in the ``cube_build`` step for routine IFU cube building.
The conversion from RA,Dec to V2,V3 is handled more accurately by
the transforms provided by :ref:`assign_wcs <assign_wcs_step>`.
Parameters
----------
ra_ref : float
RA of reference point given in degrees
dec_ref : float
Dec of reference point given in degrees
roll_ref : float
Roll angle given in degrees
v2_ref : float
V2 coordinate of reference point given in arcseconds
v3_ref : float
V3 coordinate of reference point given in arcseconds
ra : float
RA coordinate given in degrees
dec : float
Dec coordinate given in degrees
Returns
-------
v2 : float
V2 coordinate in arcseconds
v3 : float
V3 coordinate in arcseconds
"""
d2r = math.pi / 180.0
r2d = 180.0 / math.pi
v3_ref_rad = (v3_ref) / 3600.0 * d2r # convert from arcseconds to radians
v2_ref_rad = (v2_ref) / 3600.0 * d2r # convert from arcseconds to radians
roll_ref_rad = roll_ref * d2r
ra_ref_rad = ra_ref * d2r
dec_ref_rad = dec_ref * d2r
this_ra = ra * d2r
this_dec = dec * d2r
delta_ra = (this_ra - ra_ref_rad) * math.cos(dec_ref_rad)
delta_dec = this_dec - dec_ref_rad
dv2 = delta_ra * math.cos(roll_ref_rad) - delta_dec * math.sin(roll_ref_rad)
dv3 = delta_ra * math.sin(roll_ref_rad) + delta_dec * math.cos(roll_ref_rad)
v2 = v2_ref_rad + (dv2 / math.cos(v3_ref_rad))
v3 = v3_ref_rad + dv3
v2 = v2 * r2d * 3600.0 # convert from radians to arcseconds
v3 = v3 * r2d * 3600.0 # convert from radians to arcseconds
return v2, v3