Source code for galpy.util.bovy_coords

###############################################################################
#
#   bovy_coords: module for coordinate transformations between the equatorial
#                and Galactic coordinate frame
#
#
#      Main included functions:
#            radec_to_lb
#            lb_to_radec
#            radec_to_custom
#            custom_to_radec
#            lbd_to_XYZ
#            XYZ_to_lbd
#            rectgal_to_sphergal
#            sphergal_to_rectgal
#            vrpmllpmbb_to_vxvyvz
#            vxvyvz_to_vrpmllpmbb
#            pmrapmdec_to_pmllpmbb
#            pmllpmbb_to_pmrapmdec
#            pmrapmdec_to_custom
#            custom_to_pmrapmdec
#            cov_pmrapmdec_to_pmllpmbb
#            cov_dvrpmllbb_to_vxyz
#            XYZ_to_galcenrect
#            XYZ_to_galcencyl
#            galcenrect_to_XYZ
#            galcencyl_to_XYZ
#            rect_to_cyl
#            cyl_to_rect
#            rect_to_cyl_vec
#            cyl_to_rect_vec
#            vxvyvz_to_galcenrect
#            vxvyvz_to_galcencyl
#            galcenrect_to_vxvyvz
#            galcencyl_to_vxvyvz
#            dl_to_rphi_2d
#            rphi_to_dl_2d
#            Rz_to_coshucosv
#            Rz_to_uv
#            uv_to_Rz
#            Rz_to_lambdanu
#            Rz_to_lambdanu_jac
#            Rz_to_lambdanu_hess
#            lambdanu_to_Rz
#
##############################################################################
#############################################################################
#Copyright (c) 2010 - 2019, Jo Bovy
#All rights reserved.
#
#Redistribution and use in source and binary forms, with or without 
#modification, are permitted provided that the following conditions are met:
#
#   Redistributions of source code must retain the above copyright notice, 
#      this list of conditions and the following disclaimer.
#   Redistributions in binary form must reproduce the above copyright notice, 
#      this list of conditions and the following disclaimer in the 
#      documentation and/or other materials provided with the distribution.
#   The name of the author may not be used to endorse or promote products 
#      derived from this software without specific prior written permission.
#
#THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
#"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
#LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
#A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
#HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
#INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
#BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
#OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
#AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
#LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY
#WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
#POSSIBILITY OF SUCH DAMAGE.
#############################################################################
from functools import wraps
import numpy
from ..util import _rotate_to_arbitrary_vector
from ..util.config import __config__
_APY_COORDS= __config__.getboolean('astropy','astropy-coords')
_APY_LOADED= True
try:
    import astropy.coordinates as apycoords
    from astropy import units
except ImportError:
    _APY_LOADED= False
_APY_COORDS*= _APY_LOADED
_DEGTORAD= numpy.pi/180.
if _APY_LOADED:
    _K= (1.*units.mas/units.yr).to(units.km/units.s/units.kpc,
                                   equivalencies=units.dimensionless_angles())\
                                   .value
else:
    _K=4.74047
# numpy 1.14 einsum bug causes astropy conversions to fail in py2.7 -> turn off
if _APY_COORDS:
    ra, dec= numpy.array([192.25*_DEGTORAD]), numpy.array([27.4*_DEGTORAD])
    c= apycoords.SkyCoord(ra*units.rad,dec*units.rad,
                          equinox='B1950',frame='fk4')
    # This conversion fails bc of einsum bug
    try:
        c= c.transform_to(apycoords.Galactic)
    except TypeError: # pragma: no cover
        _APY_COORDS= False
def scalarDecorator(func):
    """Decorator to return scalar outputs as a set"""
    @wraps(func)
    def scalar_wrapper(*args,**kwargs):
        if numpy.array(args[0]).shape == ():
            scalarOut= True
            newargs= ()
            for ii in range(len(args)):
                newargs= newargs+(numpy.array([args[ii]]),)
            args= newargs
        else:
            scalarOut= False
        result= func(*args,**kwargs)
        if scalarOut:
            out= ()
            for ii in range(result.shape[1]):
                out= out+(result[0,ii],)
            return out
        else:
            return result
    return scalar_wrapper


def degreeDecorator(inDegrees,outDegrees):
    """
    NAME:

       degreeDecorator

    PURPOSE:

       Decorator to transform angles from and to degrees if necessary

    INPUT:

       inDegrees - list specifiying indices of angle arguments (e.g., [0,1])
       outDegrees - list, same as inDegrees, but for function return

    HISTORY:

       ____-__-__ - Written - Bovy

       2019-03-02 - speedup - Nathaniel Starkman (UofT)

    """
    # (modified) old degree decorator
    def wrapper(func):
        @wraps(func)
        def wrapped(*args, **kwargs):
            isdeg = kwargs.get('degree', False)
            if isdeg:
                args = [arg * numpy.pi / 180 if i in inDegrees else arg
                        for i, arg in enumerate(args)]
            out = func(*args, **kwargs)
            if isdeg:
                for i in outDegrees:
                    out[:, i] *= 180. / numpy.pi
            return out
        return wrapped
    return wrapper


[docs]@scalarDecorator @degreeDecorator([0,1],[0,1]) def radec_to_lb(ra,dec,degree=False,epoch=2000.0): """ NAME: radec_to_lb PURPOSE: transform from equatorial coordinates to Galactic coordinates INPUT: ra - right ascension dec - declination degree - (Bool) if True, ra and dec are given in degree and l and b will be as well epoch - epoch of ra,dec (right now only 2000.0 and 1950.0 are supported when not using astropy's transformations internally; when internally using astropy's coordinate transformations, epoch can be None for ICRS, 'JXXXX' for FK5, and 'BXXXX' for FK4) OUTPUT: l,b For vector inputs [:,2] HISTORY: 2009-11-12 - Written - Bovy (NYU) 2014-06-14 - Re-written w/ numpy functions for speed and w/ decorators for beauty - Bovy (IAS) 2016-05-13 - Added support for using astropy's coordinate transformations and for non-standard epochs - Bovy (UofT) """ if _APY_COORDS: epoch, frame= _parse_epoch_frame_apy(epoch) c= apycoords.SkyCoord(ra*units.rad,dec*units.rad, equinox=epoch,frame=frame) c= c.transform_to(apycoords.Galactic) return numpy.array([c.l.to(units.rad).value,c.b.to(units.rad).value]).T #First calculate the transformation matrix T theta,dec_ngp,ra_ngp= get_epoch_angles(epoch) T= numpy.dot(numpy.array([[numpy.cos(theta),numpy.sin(theta),0.],[numpy.sin(theta),-numpy.cos(theta),0.],[0.,0.,1.]]),numpy.dot(numpy.array([[-numpy.sin(dec_ngp),0.,numpy.cos(dec_ngp)],[0.,1.,0.],[numpy.cos(dec_ngp),0.,numpy.sin(dec_ngp)]]),numpy.array([[numpy.cos(ra_ngp),numpy.sin(ra_ngp),0.],[-numpy.sin(ra_ngp),numpy.cos(ra_ngp),0.],[0.,0.,1.]]))) #Whether to use degrees and scalar input is handled by decorators XYZ= numpy.array([numpy.cos(dec)*numpy.cos(ra), numpy.cos(dec)*numpy.sin(ra), numpy.sin(dec)]) galXYZ= numpy.dot(T,XYZ) galXYZ[2][galXYZ[2] > 1.]= 1. galXYZ[2][galXYZ[2] < -1.]= -1. b= numpy.arcsin(galXYZ[2]) l= numpy.arctan2(galXYZ[1]/numpy.cos(b),galXYZ[0]/numpy.cos(b)) l[l<0.]+= 2.*numpy.pi out= numpy.array([l,b]) return out.T
[docs]@scalarDecorator @degreeDecorator([0,1],[0,1]) def lb_to_radec(l,b,degree=False,epoch=2000.0): """ NAME: lb_to_radec PURPOSE: transform from Galactic coordinates to equatorial coordinates INPUT: l - Galactic longitude b - Galactic lattitude degree - (Bool) if True, l and b are given in degree and ra and dec will be as well epoch - epoch of ra,dec (right now only 2000.0 and 1950.0 are supported when not using astropy's transformations internally; when internally using astropy's coordinate transformations, epoch can be None for ICRS, 'JXXXX' for FK5, and 'BXXXX' for FK4) OUTPUT: ra,dec For vector inputs [:,2] HISTORY: 2010-04-07 - Written - Bovy (NYU) 2014-06-14 - Re-written w/ numpy functions for speed and w/ decorators for beauty - Bovy (IAS) 2016-05-13 - Added support for using astropy's coordinate transformations and for non-standard epochs - Bovy (UofT) """ if _APY_COORDS: epoch, frame= _parse_epoch_frame_apy(epoch) c= apycoords.SkyCoord(l*units.rad,b*units.rad,frame='galactic') if not epoch is None and 'J' in epoch: c= c.transform_to(apycoords.FK5(equinox=epoch)) elif not epoch is None and 'B' in epoch: c= c.transform_to(apycoords.FK4(equinox=epoch)) else: c= c.transform_to(apycoords.ICRS) return numpy.array([c.ra.to(units.rad).value,c.dec.to(units.rad).value]).T #First calculate the transformation matrix T' theta,dec_ngp,ra_ngp= get_epoch_angles(epoch) T= numpy.dot(numpy.array([[numpy.cos(ra_ngp),-numpy.sin(ra_ngp),0.],[numpy.sin(ra_ngp),numpy.cos(ra_ngp),0.],[0.,0.,1.]]),numpy.dot(numpy.array([[-numpy.sin(dec_ngp),0.,numpy.cos(dec_ngp)],[0.,1.,0.],[numpy.cos(dec_ngp),0.,numpy.sin(dec_ngp)]]),numpy.array([[numpy.cos(theta),numpy.sin(theta),0.],[numpy.sin(theta),-numpy.cos(theta),0.],[0.,0.,1.]]))) #Whether to use degrees and scalar input is handled by decorators XYZ= numpy.array([numpy.cos(b)*numpy.cos(l), numpy.cos(b)*numpy.sin(l), numpy.sin(b)]) eqXYZ= numpy.dot(T,XYZ) dec= numpy.arcsin(eqXYZ[2]) ra= numpy.arctan2(eqXYZ[1],eqXYZ[0]) ra[ra<0.]+= 2.*numpy.pi return numpy.array([ra,dec]).T
[docs]@scalarDecorator @degreeDecorator([0,1],[]) def lbd_to_XYZ(l,b,d,degree=False): """ NAME: lbd_to_XYZ PURPOSE: transform from spherical Galactic coordinates to rectangular Galactic coordinates (works with vector inputs) INPUT: l - Galactic longitude (rad) b - Galactic lattitude (rad) d - distance (arbitrary units) degree - (bool) if True, l and b are in degrees OUTPUT: [X,Y,Z] in whatever units d was in For vector inputs [:,3] HISTORY: 2009-10-24- Written - Bovy (NYU) 2014-06-14 - Re-written w/ numpy functions for speed and w/ decorators for beauty - Bovy (IAS) """ #Whether to use degrees and scalar input is handled by decorators return numpy.array([d*numpy.cos(b)*numpy.cos(l), d*numpy.cos(b)*numpy.sin(l), d*numpy.sin(b)]).T
[docs]def rectgal_to_sphergal(X,Y,Z,vx,vy,vz,degree=False): """ NAME: rectgal_to_sphergal PURPOSE: transform phase-space coordinates in rectangular Galactic coordinates to spherical Galactic coordinates (can take vector inputs) INPUT: X - component towards the Galactic Center (kpc) Y - component in the direction of Galactic rotation (kpc) Z - component towards the North Galactic Pole (kpc) vx - velocity towards the Galactic Center (km/s) vy - velocity in the direction of Galactic rotation (km/s) vz - velocity towards the North Galactic Pole (km/s) degree - (Bool) if True, return l and b in degrees OUTPUT: (l,b,d,vr,pmll x cos(b),pmbb) in (rad,rad,kpc,km/s,mas/yr,mas/yr) HISTORY: 2009-10-25 - Written - Bovy (NYU) """ lbd= XYZ_to_lbd(X,Y,Z,degree=degree) vrpmllpmbb= vxvyvz_to_vrpmllpmbb(vx,vy,vz,X,Y,Z,XYZ=True) if numpy.array(X).shape == (): return numpy.array([lbd[0],lbd[1],lbd[2],vrpmllpmbb[0],vrpmllpmbb[1],vrpmllpmbb[2]]) else: out=numpy.zeros((len(X),6)) out[:,0:3]= lbd out[:,3:6]= vrpmllpmbb return out
[docs]def sphergal_to_rectgal(l,b,d,vr,pmll,pmbb,degree=False): """ NAME: sphergal_to_rectgal PURPOSE: transform phase-space coordinates in spherical Galactic coordinates to rectangular Galactic coordinates (can take vector inputs) INPUT: l - Galactic longitude (rad) b - Galactic lattitude (rad) d - distance (kpc) vr - line-of-sight velocity (km/s) pmll - proper motion in the Galactic longitude direction (mu_l*cos(b) ) (mas/yr) pmbb - proper motion in the Galactic lattitude (mas/yr) degree - (bool) if True, l and b are in degrees OUTPUT: (X,Y,Z,vx,vy,vz) in (kpc,kpc,kpc,km/s,km/s,km/s) HISTORY: 2009-10-25 - Written - Bovy (NYU) """ XYZ= lbd_to_XYZ(l,b,d,degree=degree) vxvyvz= vrpmllpmbb_to_vxvyvz(vr,pmll,pmbb,l,b,d,XYZ=False,degree=degree) if numpy.array(l).shape == (): return numpy.array([XYZ[0],XYZ[1],XYZ[2],vxvyvz[0],vxvyvz[1],vxvyvz[2]]) else: out=numpy.zeros((len(l),6)) out[:,0:3]= XYZ out[:,3:6]= vxvyvz return out
[docs]@scalarDecorator @degreeDecorator([3,4],[]) def vrpmllpmbb_to_vxvyvz(vr,pmll,pmbb,l,b,d,XYZ=False,degree=False): """ NAME: vrpmllpmbb_to_vxvyvz PURPOSE: Transform velocities in the spherical Galactic coordinate frame to the rectangular Galactic coordinate frame (can take vector inputs) INPUT: vr - line-of-sight velocity (km/s) pmll - proper motion in the Galactic longitude (mu_l * cos(b))(mas/yr) pmbb - proper motion in the Galactic lattitude (mas/yr) l - Galactic longitude b - Galactic lattitude d - distance (kpc) XYZ - (bool) If True, then l,b,d is actually X,Y,Z (rectangular Galactic coordinates) degree - (bool) if True, l and b are in degrees OUTPUT: (vx,vy,vz) in (km/s,km/s,km/s) For vector inputs [:,3] HISTORY: 2009-10-24 - Written - Bovy (NYU) 2014-06-14 - Re-written w/ numpy functions for speed and w/ decorators for beauty - Bovy (IAS) """ #Whether to use degrees and scalar input is handled by decorators if XYZ: #undo the incorrect conversion that the decorator did if degree: l*= 180./numpy.pi b*= 180./numpy.pi lbd= XYZ_to_lbd(l,b,d,degree=False) l= lbd[:,0] b= lbd[:,1] d= lbd[:,2] R=numpy.zeros((3,3,len(l))) R[0,0]= numpy.cos(l)*numpy.cos(b) R[1,0]= -numpy.sin(l) R[2,0]= -numpy.cos(l)*numpy.sin(b) R[0,1]= numpy.sin(l)*numpy.cos(b) R[1,1]= numpy.cos(l) R[2,1]= -numpy.sin(l)*numpy.sin(b) R[0,2]= numpy.sin(b) R[2,2]= numpy.cos(b) invr= numpy.array([[vr,vr,vr], [d*pmll*_K,d*pmll*_K,d*pmll*_K], [d*pmbb*_K,d*pmbb*_K,d*pmbb*_K]]) return (R.T*invr.T).sum(-1)
[docs]@scalarDecorator @degreeDecorator([3,4],[]) def vxvyvz_to_vrpmllpmbb(vx,vy,vz,l,b,d,XYZ=False,degree=False): """ NAME: vxvyvz_to_vrpmllpmbb PURPOSE: Transform velocities in the rectangular Galactic coordinate frame to the spherical Galactic coordinate frame (can take vector inputs) INPUT: vx - velocity towards the Galactic Center (km/s) vy - velocity in the direction of Galactic rotation (km/s) vz - velocity towards the North Galactic Pole (km/s) l - Galactic longitude b - Galactic lattitude d - distance (kpc) XYZ - (bool) If True, then l,b,d is actually X,Y,Z (rectangular Galactic coordinates) degree - (bool) if True, l and b are in degrees OUTPUT: (vr,pmll x cos(b),pmbb) in (km/s,mas/yr,mas/yr); pmll = mu_l * cos(b) For vector inputs [:,3] HISTORY: 2009-10-24 - Written - Bovy (NYU) 2014-06-14 - Re-written w/ numpy functions for speed and w/ decorators for beauty - Bovy (IAS) """ #Whether to use degrees and scalar input is handled by decorators if XYZ: #undo the incorrect conversion that the decorator did if degree: l*= 180./numpy.pi b*= 180./numpy.pi lbd= XYZ_to_lbd(l,b,d,degree=False) l= lbd[:,0] b= lbd[:,1] d= lbd[:,2] R=numpy.zeros((3,3,len(l))) R[0,0]= numpy.cos(l)*numpy.cos(b) R[0,1]= -numpy.sin(l) R[0,2]= -numpy.cos(l)*numpy.sin(b) R[1,0]= numpy.sin(l)*numpy.cos(b) R[1,1]= numpy.cos(l) R[1,2]= -numpy.sin(l)*numpy.sin(b) R[2,0]= numpy.sin(b) R[2,2]= numpy.cos(b) invxyz= numpy.array([[vx,vx,vx], [vy,vy,vy], [vz,vz,vz]]) vrvlvb= (R.T*invxyz.T).sum(-1) vrvlvb[:,1]/= d*_K vrvlvb[:,2]/= d*_K return vrvlvb
[docs]@scalarDecorator @degreeDecorator([],[0,1]) def XYZ_to_lbd(X,Y,Z,degree=False): """ NAME: XYZ_to_lbd PURPOSE: transform from rectangular Galactic coordinates to spherical Galactic coordinates (works with vector inputs) INPUT: X - component towards the Galactic Center (in kpc; though this obviously does not matter)) Y - component in the direction of Galactic rotation (in kpc) Z - component towards the North Galactic Pole (kpc) degree - (Bool) if True, return l and b in degrees OUTPUT: [l,b,d] in (rad or degree,rad or degree,kpc) For vector inputs [:,3] HISTORY: 2009-10-24 - Written - Bovy (NYU) 2014-06-14 - Re-written w/ numpy functions for speed and w/ decorators for beauty - Bovy (IAS) """ #Whether to use degrees and scalar input is handled by decorators d= numpy.sqrt(X**2.+Y**2.+Z**2.) b=numpy.arcsin(Z/d) cosl= X/d/numpy.cos(b) sinl= Y/d/numpy.cos(b) l= numpy.arcsin(sinl) l[cosl < 0.]= numpy.pi-l[cosl < 0.] l[(cosl >= 0.)*(sinl < 0.)]+= 2.*numpy.pi out= numpy.empty((len(l),3)) out[:,0]= l out[:,1]= b out[:,2]= d return out
[docs]@scalarDecorator @degreeDecorator([2,3],[]) def pmrapmdec_to_pmllpmbb(pmra,pmdec,ra,dec,degree=False,epoch=2000.0): """ NAME: pmrapmdec_to_pmllpmbb PURPOSE: rotate proper motions in (ra,dec) into proper motions in (l,b) INPUT: pmra - proper motion in ra (multplied with cos(dec)) [mas/yr] pmdec - proper motion in dec [mas/yr] ra - right ascension dec - declination degree - if True, ra and dec are given in degrees (default=False) epoch - epoch of ra,dec (right now only 2000.0 and 1950.0 are supported when not using astropy's transformations internally; when internally using astropy's coordinate transformations, epoch can be None for ICRS, 'JXXXX' for FK5, and 'BXXXX' for FK4) OUTPUT: (pmll x cos(b),pmbb) for vector inputs [:,2] HISTORY: 2010-04-07 - Written - Bovy (NYU) 2014-06-14 - Re-written w/ numpy functions for speed and w/ decorators for beauty - Bovy (IAS) """ theta,dec_ngp,ra_ngp= get_epoch_angles(epoch) #Whether to use degrees and scalar input is handled by decorators dec[dec == dec_ngp]+= 10.**-16 #deal w/ pole. sindec_ngp= numpy.sin(dec_ngp) cosdec_ngp= numpy.cos(dec_ngp) sindec= numpy.sin(dec) cosdec= numpy.cos(dec) sinrarangp= numpy.sin(ra-ra_ngp) cosrarangp= numpy.cos(ra-ra_ngp) #These were replaced by Poleski (2013)'s equivalent form that is better at the poles #cosphi= (sindec_ngp-sindec*sinb)/cosdec/cosb #sinphi= sinrarangp*cosdec_ngp/cosb cosphi= sindec_ngp*cosdec-cosdec_ngp*sindec*cosrarangp sinphi= sinrarangp*cosdec_ngp norm= numpy.sqrt(cosphi**2.+sinphi**2.) cosphi/= norm sinphi/= norm return (numpy.array([[cosphi,-sinphi],[sinphi,cosphi]]).T\ *numpy.array([[pmra,pmra],[pmdec,pmdec]]).T).sum(-1)
[docs]@scalarDecorator @degreeDecorator([2,3],[]) def pmllpmbb_to_pmrapmdec(pmll,pmbb,l,b,degree=False,epoch=2000.0): """ NAME: pmllpmbb_to_pmrapmdec PURPOSE: rotate proper motions in (l,b) into proper motions in (ra,dec) INPUT: pmll - proper motion in l (multplied with cos(b)) [mas/yr] pmbb - proper motion in b [mas/yr] l - Galactic longitude b - Galactic lattitude degree - if True, l and b are given in degrees (default=False) epoch - epoch of ra,dec (right now only 2000.0 and 1950.0 are supported when not using astropy's transformations internally; when internally using astropy's coordinate transformations, epoch can be None for ICRS, 'JXXXX' for FK5, and 'BXXXX' for FK4) OUTPUT: (pmra x cos(dec),pmdec), for vector inputs [:,2] HISTORY: 2010-04-07 - Written - Bovy (NYU) 2014-06-14 - Re-written w/ numpy functions for speed and w/ decorators for beauty - Bovy (IAS) """ theta,dec_ngp,ra_ngp= get_epoch_angles(epoch) #Whether to use degrees and scalar input is handled by decorators radec = lb_to_radec(l,b,degree=False,epoch=epoch) ra= radec[:,0] dec= radec[:,1] dec[dec == dec_ngp]+= 10.**-16 #deal w/ pole. sindec_ngp= numpy.sin(dec_ngp) cosdec_ngp= numpy.cos(dec_ngp) sindec= numpy.sin(dec) cosdec= numpy.cos(dec) sinrarangp= numpy.sin(ra-ra_ngp) cosrarangp= numpy.cos(ra-ra_ngp) #These were replaced by Poleski (2013)'s equivalent form that is better at the poles #cosphi= (sindec_ngp-sindec*sinb)/cosdec/cosb #sinphi= sinrarangp*cosdec_ngp/cosb cosphi= sindec_ngp*cosdec-cosdec_ngp*sindec*cosrarangp sinphi= sinrarangp*cosdec_ngp norm= numpy.sqrt(cosphi**2.+sinphi**2.) cosphi/= norm sinphi/= norm return (numpy.array([[cosphi,sinphi],[-sinphi,cosphi]]).T\ *numpy.array([[pmll,pmll],[pmbb,pmbb]]).T).sum(-1)
[docs]def cov_pmrapmdec_to_pmllpmbb(cov_pmradec,ra,dec,degree=False,epoch=2000.0): """ NAME: cov_pmrapmdec_to_pmllpmbb PURPOSE: propagate the proper motions errors through the rotation from (ra,dec) to (l,b) INPUT: covar_pmradec - uncertainty covariance matrix of the proper motion in ra (multplied with cos(dec)) and dec [2,2] or [:,2,2] ra - right ascension dec - declination degree - if True, ra and dec are given in degrees (default=False) epoch - epoch of ra,dec (right now only 2000.0 and 1950.0 are supported when not using astropy's transformations internally; when internally using astropy's coordinate transformations, epoch can be None for ICRS, 'JXXXX' for FK5, and 'BXXXX' for FK4) OUTPUT: covar_pmllbb [2,2] or [:,2,2] [pmll here is pmll x cos(b)] HISTORY: 2010-04-12 - Written - Bovy (NYU) """ if len(cov_pmradec.shape) == 3: out= numpy.zeros(cov_pmradec.shape) ndata= out.shape[0] lb = radec_to_lb(ra,dec,degree=degree,epoch=epoch) for ii in range(ndata): out[ii,:,:]= cov_pmradec_to_pmllbb_single(cov_pmradec[ii,:,:], ra[ii],dec[ii],lb[ii,1], degree,epoch) return out else: l,b = radec_to_lb(ra,dec,degree=degree,epoch=epoch) return cov_pmradec_to_pmllbb_single(cov_pmradec,ra,dec,b,degree,epoch)
def cov_pmradec_to_pmllbb_single(cov_pmradec,ra,dec,b,degree=False,epoch=2000.0): """ NAME: cov_pmradec_to_pmllbb_single PURPOSE: propagate the proper motions errors through the rotation from (ra,dec) to (l,b) for scalar inputs INPUT: covar_pmradec - uncertainty covariance matrix of the proper motion in ra (multplied with cos(dec)) and dec [2,2] or [:,2,2] ra - right ascension dec - declination degree - if True, ra and dec are given in degrees (default=False) epoch - epoch of ra,dec (right now only 2000.0 and 1950.0 are supported when not using astropy's transformations internally; when internally using astropy's coordinate transformations, epoch can be None for ICRS, 'JXXXX' for FK5, and 'BXXXX' for FK4) OUTPUT: cov_pmllbb HISTORY: 2010-04-12 - Written - Bovy (NYU) """ theta,dec_ngp,ra_ngp= get_epoch_angles(epoch) if degree: sindec_ngp= numpy.sin(dec_ngp) cosdec_ngp= numpy.cos(dec_ngp) sindec= numpy.sin(dec*_DEGTORAD) cosdec= numpy.cos(dec*_DEGTORAD) sinrarangp= numpy.sin(ra*_DEGTORAD-ra_ngp) cosrarangp= numpy.cos(ra*_DEGTORAD-ra_ngp) else: sindec_ngp= numpy.sin(dec_ngp) cosdec_ngp= numpy.cos(dec_ngp) sindec= numpy.sin(dec) cosdec= numpy.cos(dec) sinrarangp= numpy.sin(ra-ra_ngp) cosrarangp= numpy.cos(ra-ra_ngp) #These were replaced by Poleski (2013)'s equivalent form that is better at the poles #cosphi= (sindec_ngp-sindec*sinb)/cosdec/cosb #sinphi= sinrarangp*cosdec_ngp/cosb cosphi= sindec_ngp*cosdec-cosdec_ngp*sindec*cosrarangp sinphi= sinrarangp*cosdec_ngp norm= numpy.sqrt(cosphi**2.+sinphi**2.) cosphi/= norm sinphi/= norm P= numpy.array([[cosphi,sinphi],[-sinphi,cosphi]]) return numpy.dot(P,numpy.dot(cov_pmradec,P.T))
[docs]def cov_dvrpmllbb_to_vxyz(d,e_d,e_vr,pmll,pmbb,cov_pmllbb,l,b, plx=False,degree=False): """ NAME: cov_dvrpmllbb_to_vxyz PURPOSE: propagate distance, radial velocity, and proper motion uncertainties to Galactic coordinates INPUT: d - distance [kpc, as/mas for plx] e_d - distance uncertainty [kpc, [as/mas] for plx] e_vr - low velocity uncertainty [km/s] pmll - proper motion in l (*cos(b)) [ [as/mas]/yr ] pmbb - proper motion in b [ [as/mas]/yr ] cov_pmllbb - uncertainty covariance for proper motion [pmll is pmll x cos(b)] l - Galactic longitude b - Galactic lattitude KEYWORDS: plx - if True, d is a parallax, and e_d is a parallax uncertainty degree - if True, l and b are given in degree OUTPUT: cov(vx,vy,vz) [3,3] or [:,3,3] HISTORY: 2010-04-12 - Written - Bovy (NYU) """ if plx: d= 1./d e_d*= d**2. if degree: l*= _DEGTORAD b*= _DEGTORAD if numpy.array(d).shape == (): return cov_dvrpmllbb_to_vxyz_single(d,e_d,e_vr,pmll,pmbb,cov_pmllbb, l,b) else: ndata= len(d) out= numpy.zeros((ndata,3,3)) for ii in range(ndata): out[ii,:,:]= cov_dvrpmllbb_to_vxyz_single(d[ii],e_d[ii],e_vr[ii], pmll[ii],pmbb[ii], cov_pmllbb[ii,:,:], l[ii],b[ii]) return out
def cov_dvrpmllbb_to_vxyz_single(d,e_d,e_vr,pmll,pmbb,cov_pmllbb,l,b): """ NAME: cov_dvrpmllbb_to_vxyz PURPOSE: propagate distance, radial velocity, and proper motion uncertainties to Galactic coordinates for scalar inputs INPUT: d - distance [kpc, as/mas for plx] e_d - distance uncertainty [kpc, [as/mas] for plx] e_vr - low velocity uncertainty [km/s] pmll - proper motion in l (*cos(b)) [ [as/mas]/yr ] pmbb - proper motion in b [ [as/mas]/yr ] cov_pmllbb - uncertainty covariance for proper motion l - Galactic longitude [rad] b - Galactic lattitude [rad] OUTPUT: cov(vx,vy,vz) [3,3] HISTORY: 2010-04-12 - Written - Bovy (NYU) """ M= _K*numpy.array([[pmll,d,0.],[pmbb,0.,d]]) cov_dpmllbb= numpy.zeros((3,3)) cov_dpmllbb[0,0]= e_d**2. cov_dpmllbb[1:3,1:3]= cov_pmllbb cov_vlvb= numpy.dot(M,numpy.dot(cov_dpmllbb,M.T)) cov_vrvlvb= numpy.zeros((3,3)) cov_vrvlvb[0,0]= e_vr**2. cov_vrvlvb[1:3,1:3]= cov_vlvb R= numpy.array([[numpy.cos(l)*numpy.cos(b), numpy.sin(l)*numpy.cos(b), numpy.sin(b)], [-numpy.sin(l),numpy.cos(l),0.], [-numpy.cos(l)*numpy.sin(b),-numpy.sin(l)*numpy.sin(b), numpy.cos(b)]]) return numpy.dot(R.T,numpy.dot(cov_vrvlvb,R))
[docs]@scalarDecorator def XYZ_to_galcenrect(X,Y,Z,Xsun=1.,Zsun=0.,_extra_rot=True): """ NAME: XYZ_to_galcenrect PURPOSE: transform XYZ coordinates (wrt Sun) to rectangular Galactocentric coordinates INPUT: X - X Y - Y Z - Z Xsun - cylindrical distance to the GC Zsun - Sun's height above the midplane _extra_rot= (True) if True, perform an extra tiny rotation to align the Galactocentric coordinate frame with astropy's definition OUTPUT: (Xg, Yg, Zg) HISTORY: 2010-09-24 - Written - Bovy (NYU) 2016-05-12 - Edited to properly take into account the Sun's vertical position; dropped Ysun keyword - Bovy (UofT) 2018-04-18 - Tweaked to be consistent with astropy's Galactocentric frame - Bovy (UofT) """ if _extra_rot: X,Y,Z= numpy.dot(galcen_extra_rot,numpy.array([X,Y,Z])) dgc= numpy.sqrt(Xsun**2.+Zsun**2.) costheta, sintheta= Xsun/dgc, Zsun/dgc return numpy.dot(numpy.array([[costheta,0.,-sintheta], [0.,1.,0.], [sintheta,0.,costheta]]), numpy.array([-X+dgc,Y,numpy.sign(Xsun)*Z])).T
[docs]@scalarDecorator def galcenrect_to_XYZ(X,Y,Z,Xsun=1.,Zsun=0.,_extra_rot=True): """ NAME: galcenrect_to_XYZ PURPOSE: transform rectangular Galactocentric to XYZ coordinates (wrt Sun) coordinates INPUT: X, Y, Z - Galactocentric rectangular coordinates Xsun - cylindrical distance to the GC (can be array of same length as X) Zsun - Sun's height above the midplane (can be array of same length as X) _extra_rot= (True) if True, perform an extra tiny rotation to align the Galactocentric coordinate frame with astropy's definition OUTPUT: (X, Y, Z) HISTORY: 2011-02-23 - Written - Bovy (NYU) 2016-05-12 - Edited to properly take into account the Sun's vertical position; dropped Ysun keyword - Bovy (UofT) 2017-10-24 - Allowed Xsun/Zsun to be arrays - Bovy (UofT) 2018-04-18 - Tweaked to be consistent with astropy's Galactocentric frame - Bovy (UofT) """ dgc= numpy.sqrt(Xsun**2.+Zsun**2.) costheta, sintheta= Xsun/dgc, Zsun/dgc if isinstance(Xsun,numpy.ndarray): zero= numpy.zeros(len(Xsun)) one= numpy.ones(len(Xsun)) Carr= numpy.rollaxis(numpy.array([[-costheta,zero,-sintheta], [zero,one,zero], [-numpy.sign(Xsun)*sintheta,zero, numpy.sign(Xsun)*costheta]]),2) out= ((Carr*numpy.array([[X,X,X],[Y,Y,Y],[Z,Z,Z]]).T).sum(-1) +numpy.array([dgc,zero,zero]).T) else: out= numpy.dot(numpy.array([[-costheta,0.,-sintheta], [0.,1.,0.], [-numpy.sign(Xsun)*sintheta,0., numpy.sign(Xsun)*costheta]]), numpy.array([X,Y,Z])).T+numpy.array([dgc,0.,0.]) if _extra_rot: return numpy.dot(galcen_extra_rot.T,out.T).T else: return out
[docs]def rect_to_cyl(X,Y,Z): """ NAME: rect_to_cyl PURPOSE: convert from rectangular to cylindrical coordinates INPUT: X, Y, Z - rectangular coordinates OUTPUT: R,phi,z HISTORY: 2010-09-24 - Written - Bovy (NYU) 2019-06-21 - Changed such that phi in [-pi,pi] - Bovy (UofT) """ return (numpy.sqrt(X**2.+Y**2.),numpy.arctan2(Y,X),Z)
[docs]def cyl_to_rect(R,phi,Z): """ NAME: cyl_to_rect PURPOSE: convert from cylindrical to rectangular coordinates INPUT: R, phi, Z - cylindrical coordinates OUTPUT: X,Y,Z HISTORY: 2011-02-23 - Written - Bovy (NYU) """ return (R*numpy.cos(phi),R*numpy.sin(phi),Z)
def cyl_to_spher(R,Z, phi): """ NAME: cyl_to_spher PURPOSE: convert from cylindrical to spherical coordinates INPUT: R, Z, phi- cylindrical coordinates OUTPUT: R, theta, phi - spherical coordinates HISTORY: 2016-05-16 - Written - Aladdin """ theta = numpy.arctan2(R, Z) r = (R**2 + Z**2)**.5 return (r,theta, phi) def spher_to_cyl(r, theta, phi): """ NAME: spher_to_cyl PURPOSE: convert from spherical to cylindrical coordinates INPUT: r, theta, phi - spherical coordinates OUTPUT: R, z, phi - spherical coordinates HISTORY: 2016-05-20 - Written - Aladdin """ R = r*numpy.sin(theta) z = r*numpy.cos(theta) return (R,z, phi)
[docs]@scalarDecorator def XYZ_to_galcencyl(X,Y,Z,Xsun=1.,Zsun=0.,_extra_rot=True): """ NAME: XYZ_to_galcencyl PURPOSE: transform XYZ coordinates (wrt Sun) to cylindrical Galactocentric coordinates INPUT: X - X Y - Y Z - Z Xsun - cylindrical distance to the GC Zsun - Sun's height above the midplane _extra_rot= (True) if True, perform an extra tiny rotation to align the Galactocentric coordinate frame with astropy's definition OUTPUT: R,phi,z HISTORY: 2010-09-24 - Written - Bovy (NYU) """ XYZ= numpy.atleast_2d(XYZ_to_galcenrect(X,Y,Z,Xsun=Xsun,Zsun=Zsun, _extra_rot=_extra_rot)) return numpy.array(rect_to_cyl(XYZ[:,0],XYZ[:,1],XYZ[:,2])).T
[docs]@scalarDecorator def galcencyl_to_XYZ(R,phi,Z,Xsun=1.,Zsun=0.,_extra_rot=True): """ NAME: galcencyl_to_XYZ PURPOSE: transform cylindrical Galactocentric coordinates to XYZ coordinates (wrt Sun) INPUT: R, phi, Z - Galactocentric cylindrical coordinates Xsun - cylindrical distance to the GC (can be array of same length as R) Zsun - Sun's height above the midplane (can be array of same length as R) _extra_rot= (True) if True, perform an extra tiny rotation to align the Galactocentric coordinate frame with astropy's definition OUTPUT: X,Y,Z HISTORY: 2011-02-23 - Written - Bovy (NYU) 2017-10-24 - Allowed Xsun/Zsun to be arrays - Bovy (UofT) """ Xr,Yr,Zr= cyl_to_rect(R,phi,Z) return galcenrect_to_XYZ(Xr,Yr,Zr,Xsun=Xsun,Zsun=Zsun, _extra_rot=_extra_rot)
[docs]@scalarDecorator def vxvyvz_to_galcenrect(vx,vy,vz,vsun=[0.,1.,0.],Xsun=1.,Zsun=0., _extra_rot=True): """ NAME: vxvyvz_to_galcenrect PURPOSE: transform velocities in XYZ coordinates (wrt Sun) to rectangular Galactocentric coordinates for velocities INPUT: vx - U vy - V vz - W vsun - velocity of the sun in the GC frame ndarray[3] Xsun - cylindrical distance to the GC Zsun - Sun's height above the midplane _extra_rot= (True) if True, perform an extra tiny rotation to align the Galactocentric coordinate frame with astropy's definition OUTPUT: [:,3]= vXg, vYg, vZg HISTORY: 2010-09-24 - Written - Bovy (NYU) 2016-05-12 - Edited to properly take into account the Sun's vertical position; dropped Ysun keyword - Bovy (UofT) 2018-04-18 - Tweaked to be consistent with astropy's Galactocentric frame - Bovy (UofT) """ if _extra_rot: vx,vy,vz= numpy.dot(galcen_extra_rot,numpy.array([vx,vy,vz])) dgc= numpy.sqrt(Xsun**2.+Zsun**2.) costheta, sintheta= Xsun/dgc, Zsun/dgc return numpy.dot(numpy.array([[costheta,0.,-sintheta], [0.,1.,0.], [sintheta,0.,costheta]]), numpy.array([-vx,vy,numpy.sign(Xsun)*vz])).T+numpy.array(vsun)
[docs]@scalarDecorator def vxvyvz_to_galcencyl(vx,vy,vz,X,Y,Z,vsun=[0.,1.,0.],Xsun=1.,Zsun=0., galcen=False,_extra_rot=True): """ NAME: vxvyvz_to_galcencyl PURPOSE: transform velocities in XYZ coordinates (wrt Sun) to cylindrical Galactocentric coordinates for velocities INPUT: vx - U vy - V vz - W X - X in Galactocentric rectangular coordinates Y - Y in Galactocentric rectangular coordinates Z - Z in Galactocentric rectangular coordinates vsun - velocity of the sun in the GC frame ndarray[3] Xsun - cylindrical distance to the GC Zsun - Sun's height above the midplane galcen - if True, then X,Y,Z are in cylindrical Galactocentric coordinates rather than rectangular coordinates _extra_rot= (True) if True, perform an extra tiny rotation to align the Galactocentric coordinate frame with astropy's definition OUTPUT: vRg, vTg, vZg HISTORY: 2010-09-24 - Written - Bovy (NYU) """ vxyz= vxvyvz_to_galcenrect(vx,vy,vz,vsun=vsun,Xsun=Xsun,Zsun=Zsun, _extra_rot=_extra_rot) return numpy.array(\ rect_to_cyl_vec(vxyz[:,0],vxyz[:,1],vxyz[:,2],X,Y,Z,cyl=galcen)).T
[docs]@scalarDecorator def galcenrect_to_vxvyvz(vXg,vYg,vZg,vsun=[0.,1.,0.],Xsun=1.,Zsun=0., _extra_rot=True): """ NAME: galcenrect_to_vxvyvz PURPOSE: transform rectangular Galactocentric coordinates to XYZ coordinates (wrt Sun) for velocities INPUT: vXg - Galactocentric x-velocity vYg - Galactocentric y-velocity vZg - Galactocentric z-velocity vsun - velocity of the sun in the GC frame ndarray[3] (can be array of same length as vXg; shape [3,N]) Xsun - cylindrical distance to the GC (can be array of same length as vXg) Zsun - Sun's height above the midplane (can be array of same length as vXg) _extra_rot= (True) if True, perform an extra tiny rotation to align the Galactocentric coordinate frame with astropy's definition OUTPUT: [:,3]= vx, vy, vz HISTORY: 2011-02-24 - Written - Bovy (NYU) 2016-05-12 - Edited to properly take into account the Sun's vertical position; dropped Ysun keyword - Bovy (UofT) 2017-10-24 - Allowed vsun/Xsun/Zsun to be arrays - Bovy (UofT) 2018-04-18 - Tweaked to be consistent with astropy's Galactocentric frame - Bovy (UofT) """ dgc= numpy.sqrt(Xsun**2.+Zsun**2.) costheta, sintheta= Xsun/dgc, Zsun/dgc if isinstance(Xsun,numpy.ndarray): zero= numpy.zeros(len(Xsun)) one= numpy.ones(len(Xsun)) Carr= numpy.rollaxis(numpy.array([[-costheta,zero,-sintheta], [zero,one,zero], [-numpy.sign(Xsun)*sintheta,zero, numpy.sign(Xsun)*costheta]]),2) out= ((Carr *numpy.array([[vXg-vsun[0],vXg-vsun[0],vXg-vsun[0]], [vYg-vsun[1],vYg-vsun[1],vYg-vsun[1]], [vZg-vsun[2],vZg-vsun[2],vZg-vsun[2]]]).T).sum(-1)) else: out= numpy.dot(numpy.array([[-costheta,0.,-sintheta], [0.,1.,0.], [-numpy.sign(Xsun)*sintheta,0., numpy.sign(Xsun)*costheta]]), numpy.array([vXg-vsun[0],vYg-vsun[1],vZg-vsun[2]])).T if _extra_rot: return numpy.dot(galcen_extra_rot.T,out.T).T else: return out
[docs]@scalarDecorator def galcencyl_to_vxvyvz(vR,vT,vZ,phi,vsun=[0.,1.,0.],Xsun=1.,Zsun=0., _extra_rot=True): """ NAME: galcencyl_to_vxvyvz PURPOSE: transform cylindrical Galactocentric coordinates to XYZ (wrt Sun) coordinates for velocities INPUT: vR - Galactocentric radial velocity vT - Galactocentric tangential velocity vZ - Galactocentric vertical velocity phi - Galactocentric azimuth vsun - velocity of the sun in the GC frame ndarray[3] (can be array of same length as vRg; shape [3,N]) Xsun - cylindrical distance to the GC (can be array of same length as vRg) Zsun - Sun's height above the midplane (can be array of same length as vRg) _extra_rot= (True) if True, perform an extra tiny rotation to align the Galactocentric coordinate frame with astropy's definition OUTPUT: vx,vy,vz HISTORY: 2011-02-24 - Written - Bovy (NYU) 2017-10-24 - Allowed vsun/Xsun/Zsun to be arrays - Bovy (NYU) """ vXg, vYg, vZg= cyl_to_rect_vec(vR,vT,vZ,phi) return galcenrect_to_vxvyvz(vXg,vYg,vZg,vsun=vsun,Xsun=Xsun,Zsun=Zsun, _extra_rot=_extra_rot)
[docs]def rect_to_cyl_vec(vx,vy,vz,X,Y,Z,cyl=False): """ NAME: rect_to_cyl_vec PURPOSE: transform vectors from rectangular to cylindrical coordinates vectors INPUT: vx - vy - vz - X - X Y - Y Z - Z cyl - if True, X,Y,Z are already cylindrical OUTPUT: vR,vT,vz HISTORY: 2010-09-24 - Written - Bovy (NYU) """ if not cyl: R,phi,Z= rect_to_cyl(X,Y,Z) else: phi= Y vr=+vx*numpy.cos(phi)+vy*numpy.sin(phi) vt= -vx*numpy.sin(phi)+vy*numpy.cos(phi) return (vr,vt,vz)
[docs]def cyl_to_rect_vec(vr,vt,vz,phi): """ NAME: cyl_to_rect_vec PURPOSE: transform vectors from cylindrical to rectangular coordinate vectors INPUT: vr - radial velocity vt - tangential velocity vz - vertical velocity phi - azimuth OUTPUT: vx,vy,vz HISTORY: 2011-02-24 - Written - Bovy (NYU) """ vx= vr*numpy.cos(phi)-vt*numpy.sin(phi) vy= vr*numpy.sin(phi)+vt*numpy.cos(phi) return (vx,vy,vz)
def cyl_to_rect_jac(*args): """ NAME: cyl_to_rect_jac PURPOSE: calculate the Jacobian of the cylindrical to rectangular conversion INPUT: R, phi, Z- cylindrical coordinates vR, vT, vZ- cylindrical velocities if 6 inputs: R,vR,vT,z,vz,phi if 3: R, phi, Z OUTPUT: jacobian d(rect)/d(cyl) HISTORY: 2013-12-09 - Written - Bovy (IAS) """ out= numpy.zeros((6,6)) if len(args) == 3: R, phi, Z= args vR, vT, vZ= 0., 0., 0. outIndx= numpy.array([True,False,False,True,False,True],dtype='bool') elif len(args) == 6: R, vR, vT, Z, vZ, phi= args outIndx= numpy.ones(6,dtype='bool') cp= numpy.cos(phi) sp= numpy.sin(phi) out[0,0]= cp out[0,5]= -R*sp out[1,0]= sp out[1,5]= R*cp out[2,3]= 1. out[3,1]= cp out[3,2]= -sp out[3,5]= -vT*cp-vR*sp out[4,1]= sp out[4,2]= cp out[4,5]= -vT*sp+vR*cp out[5,4]= 1. if len(args) == 3: out= out[:3,outIndx] out[:,[1,2]]= out[:,[2,1]] return out def galcenrect_to_XYZ_jac(*args,**kwargs): """ NAME: galcenrect_to_XYZ_jac PURPOSE: calculate the Jacobian of the Galactocentric rectangular to Galactic coordinates INPUT: X,Y,Z- Galactocentric rectangular coordinates vX, vY, vZ- Galactocentric rectangular velocities if 6 inputs: X,Y,Z,vX,vY,vZ if 3: X,Y,Z Xsun - cylindrical distance to the GC Zsun - Sun's height above the midplane OUTPUT: jacobian d(galcen.)/d(Galactic) HISTORY: 2013-12-09 - Written - Bovy (IAS) """ Xsun= kwargs.get('Xsun',1.) dgc= numpy.sqrt(Xsun**2.+kwargs.get('Zsun',0.)**2.) costheta, sintheta= Xsun/dgc, kwargs.get('Zsun',0.)/dgc out= numpy.zeros((6,6)) out[0,0]= -costheta out[0,2]= -sintheta out[1,1]= 1. out[2,0]= -numpy.sign(Xsun)*sintheta out[2,2]= numpy.sign(Xsun)*costheta if len(args) == 3: return out[:3,:3] out[3,3]= -costheta out[3,5]= -sintheta out[4,4]= 1. out[5,3]= -numpy.sign(Xsun)*sintheta out[5,5]= numpy.sign(Xsun)*costheta return out def lbd_to_XYZ_jac(*args,**kwargs): """ NAME: lbd_to_XYZ_jac PURPOSE: calculate the Jacobian of the Galactic spherical coordinates to Galactic rectangular coordinates transformation INPUT: l,b,D- Galactic spherical coordinates vlos,pmll,pmbb- Galactic spherical velocities (some as proper motions) if 6 inputs: l,b,D,vlos,pmll x cos(b),pmbb if 3: l,b,D degree= (False) if True, l and b are in degrees OUTPUT: jacobian HISTORY: 2013-12-09 - Written - Bovy (IAS) """ out= numpy.zeros((6,6)) if len(args) == 3: l,b,D= args vlos, pmll, pmbb= 0., 0., 0. elif len(args) == 6: l,b,D,vlos,pmll,pmbb= args if kwargs.get('degree',False): l*= _DEGTORAD b*= _DEGTORAD cl= numpy.cos(l) sl= numpy.sin(l) cb= numpy.cos(b) sb= numpy.sin(b) out[0,0]= -D*cb*sl out[0,1]= -D*sb*cl out[0,2]= cb*cl out[1,0]= D*cb*cl out[1,1]= -D*sb*sl out[1,2]= cb*sl out[2,1]= D*cb out[2,2]= sb if len(args) == 3: if kwargs.get('degree',False): out[:,0]*= _DEGTORAD out[:,1]*= _DEGTORAD return out[:3,:3] out[3,0]= -sl*cb*vlos-cl*_K*D*pmll+sb*sl*_K*D*pmbb out[3,1]= -cl*sb*vlos-cb*cl*_K*D*pmbb out[3,2]= -sl*_K*pmll-sb*cl*_K*pmbb out[3,3]= cl*cb out[3,4]= -sl*_K*D out[3,5]= -cl*sb*_K*D out[4,0]= cl*cb*vlos-sl*_K*D*pmll-cl*sb*_K*D*pmbb out[4,1]= -sl*sb*vlos-sl*cb*_K*D*pmbb out[4,2]= cl*_K*pmll-sl*sb*_K*pmbb out[4,3]= sl*cb out[4,4]= cl*_K*D out[4,5]= -sl*sb*_K*D out[5,1]= cb*vlos-sb*_K*D*pmbb out[5,2]= cb*_K*pmbb out[5,3]= sb out[5,5]= cb*_K*D if kwargs.get('degree',False): out[:,0]*= _DEGTORAD out[:,1]*= _DEGTORAD return out
[docs]def dl_to_rphi_2d(d,l,degree=False,ro=1.,phio=0.): """ NAME: dl_to_rphi_2d PURPOSE: convert Galactic longitude and distance to Galactocentric radius and azimuth INPUT: d - distance l - Galactic longitude [rad/deg if degree] KEYWORDS: degree= (False): l is in degrees rather than rad ro= (1) Galactocentric radius of the observer phio= (0) Galactocentric azimuth of the observer [rad/deg if degree] OUTPUT: (R,phi); phi in degree if degree HISTORY: 2012-01-04 - Written - Bovy (IAS) """ scalarOut, listOut= False, False if isinstance(d,(int,float)): d= numpy.array([d]) scalarOut= True elif isinstance(d,list): d= numpy.array(d) listOut= True if isinstance(l,(int,float)): l= numpy.array([l]) elif isinstance(l,list): l= numpy.array(l) if degree: l*= _DEGTORAD R= numpy.sqrt(ro**2.+d**2.-2.*d*ro*numpy.cos(l)) phi= numpy.arcsin(d/R*numpy.sin(l)) indx= (ro/numpy.cos(l) < d)*(numpy.cos(l) > 0.) phi[indx]= numpy.pi-numpy.arcsin(d[indx]/R[indx]*numpy.sin(l[indx])) if degree: phi/= _DEGTORAD phi+= phio if scalarOut: return (R[0],phi[0]) elif listOut: return (list(R),list(phi)) else: return (R,phi)
[docs]def rphi_to_dl_2d(R,phi,degree=False,ro=1.,phio=0.): """ NAME: rphi_to_dl_2d PURPOSE: convert Galactocentric radius and azimuth to distance and Galactic longitude INPUT: R - Galactocentric radius phi - Galactocentric azimuth [rad/deg if degree] KEYWORDS: degree= (False): phi is in degrees rather than rad ro= (1) Galactocentric radius of the observer phio= (0) Galactocentric azimuth of the observer [rad/deg if degree] OUTPUT: (d,l); phi in degree if degree HISTORY: 2012-01-04 - Written - Bovy (IAS) """ scalarOut, listOut= False, False if isinstance(R,(int,float)): R= numpy.array([R]) scalarOut= True elif isinstance(R,list): R= numpy.array(R) listOut= True if isinstance(phi,(int,float)): phi= numpy.array([phi]) elif isinstance(phi,list): phi= numpy.array(phi) phi-= phio if degree: phi*= _DEGTORAD d= numpy.sqrt(R**2.+ro**2.-2.*R*ro*numpy.cos(phi)) l= numpy.arcsin(R/d*numpy.sin(phi)) indx= (ro/numpy.cos(phi) < R)*(numpy.cos(phi) > 0.) l[indx]= numpy.pi-numpy.arcsin(R[indx]/d[indx]*numpy.sin(phi[indx])) if degree: l/= _DEGTORAD if scalarOut: return (d[0],l[0]) elif listOut: return (list(d),list(l)) else: return (d,l)
[docs]def Rz_to_coshucosv(R,z,delta=1.,oblate=False): """ NAME: Rz_to_coshucosv PURPOSE: calculate prolate confocal cosh(u) and cos(v) coordinates from R,z, and delta INPUT: R - radius z - height delta= focus oblate= (False) if True, compute oblate confocal coordinates instead of prolate OUTPUT: (cosh(u),cos(v)) HISTORY: 2012-11-27 - Written - Bovy (IAS) 2017-10-11 - Added oblate coordinates - Bovy (UofT) """ if oblate: d12= (R+delta)**2.+z**2. d22= (R-delta)**2.+z**2. else: d12= (z+delta)**2.+R**2. d22= (z-delta)**2.+R**2. coshu= 0.5/delta*(numpy.sqrt(d12)+numpy.sqrt(d22)) cosv= 0.5/delta*(numpy.sqrt(d12)-numpy.sqrt(d22)) if oblate: # cosv is currently really sinv cosv= numpy.sqrt(1.-cosv**2.) return (coshu,cosv)
[docs]def Rz_to_uv(R,z,delta=1.,oblate=False): """ NAME: Rz_to_uv PURPOSE: calculate prolate or oblate confocal u and v coordinates from R,z, and delta INPUT: R - radius z - height delta= focus oblate= (False) if True, compute oblate confocal coordinates instead of prolate OUTPUT: (u,v) HISTORY: 2012-11-27 - Written - Bovy (IAS) 2017-10-11 - Added oblate coordinates - Bovy (UofT) """ coshu, cosv= Rz_to_coshucosv(R,z,delta,oblate=oblate) u= numpy.arccosh(coshu) v= numpy.arccos(cosv) return (u,v)
[docs]def uv_to_Rz(u,v,delta=1.,oblate=False): """ NAME: uv_to_Rz PURPOSE: calculate R and z from prolate confocal u and v coordinates INPUT: u - confocal u v - confocal v delta= focus oblate= (False) if True, compute oblate confocal coordinates instead of prolate OUTPUT: (R,z) HISTORY: 2012-11-27 - Written - Bovy (IAS) 2017-10-11 - Added oblate coordinates - Bovy (UofT) """ if oblate: R= delta*numpy.cosh(u)*numpy.sin(v) z= delta*numpy.sinh(u)*numpy.cos(v) else: R= delta*numpy.sinh(u)*numpy.sin(v) z= delta*numpy.cosh(u)*numpy.cos(v) return (R,z)
[docs]def vRvz_to_pupv(vR,vz,R,z,delta=1.,oblate=False,uv=False): """ NAME: vRvz_to_pupv PURPOSE: calculate momenta in prolate or oblate confocal u and v coordinates from cylindrical velocities vR,vz for a given focal length delta INPUT: vR - radial velocity in cylindrical coordinates vz - vertical velocity in cylindrical coordinates R - radius z - height delta= focus oblate= (False) if True, compute oblate confocal coordinates instead of prolate uv= (False) if True, the given R,z are actually u,v OUTPUT: (pu,pv) HISTORY: 2017-11-28 - Written - Bovy (UofT) """ if not uv: u,v= Rz_to_uv(R,z,delta,oblate=oblate) else: u,v= R,z if oblate: pu= delta*(vR*numpy.sinh(u)*numpy.sin(v)+vz*numpy.cosh(u)*numpy.cos(v)) pv= delta*(vR*numpy.cosh(u)*numpy.cos(v)-vz*numpy.sinh(u)*numpy.sin(v)) else: pu= delta*(vR*numpy.cosh(u)*numpy.sin(v)+vz*numpy.sinh(u)*numpy.cos(v)) pv= delta*(vR*numpy.sinh(u)*numpy.cos(v)-vz*numpy.cosh(u)*numpy.sin(v)) return (pu,pv)
[docs]def pupv_to_vRvz(pu,pv,u,v,delta=1.,oblate=False): """ NAME: pupv_to_vRvz PURPOSE: calculate cylindrical vR and vz from momenta in prolate or oblate confocal u and v coordinates for a given focal length delta INPUT: pu - u momentum pv - v momentum u - u coordinate v - v coordinate delta= focus oblate= (False) if True, compute oblate confocal coordinates instead of prolate OUTPUT: (vR,vz) HISTORY: 2017-12-04 - Written - Bovy (UofT) """ if oblate: denom= delta*(numpy.sinh(u)**2.+numpy.cos(v)**2.) vR= (pu*numpy.sinh(u)*numpy.sin(v)+pv*numpy.cosh(u)*numpy.cos(v))/denom vz= (pu*numpy.cosh(u)*numpy.cos(v)-pv*numpy.sinh(u)*numpy.sin(v))/denom else: denom= delta*(numpy.sinh(u)**2.+numpy.sin(v)**2.) vR= (pu*numpy.cosh(u)*numpy.sin(v)+pv*numpy.sinh(u)*numpy.cos(v))/denom vz= (pu*numpy.sinh(u)*numpy.cos(v)-pv*numpy.cosh(u)*numpy.sin(v))/denom return (vR,vz)
def Rz_to_lambdanu(R,z,ac=5.,Delta=1.): """ NAME: Rz_to_lambdanu PURPOSE: calculate the prolate spheroidal coordinates (lambda,nu) from galactocentric cylindrical coordinates (R,z) by solving eq. (2.2) in Dejonghe & de Zeeuw (1988a) for (lambda,nu): R^2 = (l+a) * (n+a) / (a-g) z^2 = (l+g) * (n+g) / (g-a) Delta^2 = g-a INPUT: R - Galactocentric cylindrical radius z - vertical height ac - axis ratio of the coordinate surfaces (a/c) = sqrt(-a) / sqrt(-g) (default: 5.) Delta - focal distance that defines the spheroidal coordinate system (default: 1.) Delta=sqrt(g-a) OUTPUT: (lambda,nu) HISTORY: 2015-02-13 - Written - Trick (MPIA) """ g = Delta**2 / (1.-ac**2) a = g - Delta**2 term = R**2 + z**2 - a - g discr = (R**2 + z**2 - Delta**2)**2 + (4. * Delta**2 * R**2) l = 0.5 * (term + numpy.sqrt(discr)) n = 0.5 * (term - numpy.sqrt(discr)) if isinstance(z,float) and z == 0.: l = R**2 - a n = -g elif isinstance(z,numpy.ndarray) and numpy.sum(z == 0.) > 0: if isinstance(R,float): l[z==0.] = R**2 - a if isinstance(R,numpy.ndarray): l[z==0.] = R[z==0.]**2 - a n[z==0.] = -g return (l,n) def Rz_to_lambdanu_jac(R,z,Delta=1.): """ NAME: Rz_to_lambdanu_jac PURPOSE: calculate the Jacobian of the cylindrical (R,z) to prolate spheroidal (lambda,nu) conversion INPUT: R - Galactocentric cylindrical radius z - vertical height Delta - focal distance that defines the spheroidal coordinate system (default: 1.) Delta=sqrt(g-a) OUTPUT: jacobian d((lambda,nu))/d((R,z)) HISTORY: 2015-02-13 - Written - Trick (MPIA) """ discr = (R**2 + z**2 - Delta**2)**2 + (4. * Delta**2 * R**2) dldR = R * (1. + (R**2 + z**2 + Delta**2) / numpy.sqrt(discr)) dndR = R * (1. - (R**2 + z**2 + Delta**2) / numpy.sqrt(discr)) dldz = z * (1. + (R**2 + z**2 - Delta**2) / numpy.sqrt(discr)) dndz = z * (1. - (R**2 + z**2 - Delta**2) / numpy.sqrt(discr)) dim = 1 if isinstance(R,numpy.ndarray): dim = len(R) elif isinstance(z,numpy.ndarray): dim = len(z) jac = numpy.zeros((2,2,dim)) jac[0,0,:] = dldR jac[0,1,:] = dldz jac[1,0,:] = dndR jac[1,1,:] = dndz if dim == 1: return jac[:,:,0] else: return jac def Rz_to_lambdanu_hess(R,z,Delta=1.): """ NAME: Rz_to_lambdanu_hess PURPOSE: calculate the Hessian of the cylindrical (R,z) to prolate spheroidal (lambda,nu) conversion INPUT: R - Galactocentric cylindrical radius z - vertical height Delta - focal distance that defines the spheroidal coordinate system (default: 1.) Delta=sqrt(g-a) OUTPUT: hessian [d^2(lamda)/d(R,z)^2 , d^2(nu)/d(R,z)^2] HISTORY: 2015-02-13 - Written - Trick (MPIA) """ D = Delta R2 = R**2 z2 = z**2 D2 = D**2 discr = (R2 + z2 - D2)**2 + (4. * D2 * R2) d2ldR2 = 1. + (3.*R2+ z2+D2)/discr**0.5 - (2.*R2*(R2+z2+D2)**2)/discr**1.5 d2ndR2 = 1. - (3.*R2+ z2+D2)/discr**0.5 + (2.*R2*(R2+z2+D2)**2)/discr**1.5 d2ldz2 = 1. + ( R2+3.*z2-D2)/discr**0.5 - (2.*z2*(R2+z2-D2)**2)/discr**1.5 d2ndz2 = 1. - ( R2+3.*z2-D2)/discr**0.5 + (2.*z2*(R2+z2-D2)**2)/discr**1.5 d2ldRdz = 2.*R*z/discr**0.5 * ( 1. - ((R2+z2)**2-D**4)/discr) d2ndRdz = 2.*R*z/discr**0.5 * (-1. + ((R2+z2)**2-D**4)/discr) dim = 1 if isinstance(R,numpy.ndarray): dim = len(R) elif isinstance(z,numpy.ndarray): dim = len(z) hess = numpy.zeros((2,2,2,dim)) #Hessian for lambda: hess[0,0,0,:] = d2ldR2 hess[0,0,1,:] = d2ldRdz hess[0,1,0,:] = d2ldRdz hess[0,1,1,:] = d2ldz2 #Hessian for nu: hess[1,0,0,:] = d2ndR2 hess[1,0,1,:] = d2ndRdz hess[1,1,0,:] = d2ndRdz hess[1,1,1,:] = d2ndz2 if dim == 1: return hess[:,:,:,0] else: return hess def lambdanu_to_Rz(l,n,ac=5.,Delta=1.): """ NAME: lambdanu_to_Rz PURPOSE: calculate galactocentric cylindrical coordinates (R,z) from prolate spheroidal coordinates (lambda,nu), cf. eq. (2.2) in Dejonghe & de Zeeuw (1988a) INPUT: l - prolate spheroidal coordinate lambda n - prolate spheroidal coordinate nu ac - axis ratio of the coordinate surfaces (a/c) = sqrt(-a) / sqrt(-g) (default: 5.) Delta - focal distance that defines the spheroidal coordinate system (default: 1.) Delta=sqrt(g-a) OUTPUT: (R,z) HISTORY: 2015-02-13 - Written - Trick (MPIA) """ g = Delta**2 / (1.-ac**2) a = g - Delta**2 r2 = (l + a) * (n + a) / (a - g) z2 = (l + g) * (n + g) / (g - a) index = (r2 < 0.) * ((n+a) > 0.) * ((n+a) < 1e-10) if numpy.any(index): if isinstance(r2,numpy.ndarray): r2[index] = 0. else: r2 = 0. index = (z2 < 0.) * ((n+g) < 0.) * ((n+g) > -1e-10) if numpy.any(index): if isinstance(z2,numpy.ndarray): z2[index] = 0. else: z2 = 0. return (numpy.sqrt(r2),numpy.sqrt(z2))
[docs]@scalarDecorator @degreeDecorator([0,1],[0,1]) def radec_to_custom(ra,dec,T=None,degree=False): """ NAME: radec_to_custom PURPOSE: transform from equatorial coordinates to a custom set of sky coordinates INPUT: ra - right ascension dec - declination T= matrix defining the transformation: new_rect= T dot old_rect, where old_rect = [cos(dec)cos(ra),cos(dec)sin(ra),sin(dec)] and similar for new_rect degree - (Bool) if True, ra and dec are given in degree and l and b will be as well OUTPUT: custom longitude, custom latitude (with longitude -180 to 180) For vector inputs [:,2] HISTORY: 2009-11-12 - Written - Bovy (NYU) 2014-06-14 - Re-written w/ numpy functions for speed and w/ decorators for beauty - Bovy (IAS) 2019-03-02 - adjusted angle ranges - Nathaniel (UofT) """ if T is None: raise ValueError("Must set T= for radec_to_custom") #Whether to use degrees and scalar input is handled by decorators XYZ= numpy.array([numpy.cos(dec)*numpy.cos(ra), numpy.cos(dec)*numpy.sin(ra), numpy.sin(dec)]) galXYZ= numpy.dot(T,XYZ) b= numpy.arcsin(galXYZ[2]) # [-pi/2, pi/2] l= numpy.arctan2(galXYZ[1], galXYZ[0]) l[l<0] += 2 * numpy.pi # fix range to [0, 2 pi] out= numpy.array([l,b]) return out.T
[docs]@scalarDecorator @degreeDecorator([2,3],[]) def pmrapmdec_to_custom(pmra,pmdec,ra,dec,T=None,degree=False): """ NAME: pmrapmdec_to_custom PURPOSE: rotate proper motions in (ra,dec) to proper motions in a custom set of sky coordinates (phi1,phi2) INPUT: pmra - proper motion in ra (multplied with cos(dec)) [mas/yr] pmdec - proper motion in dec [mas/yr] ra - right ascension dec - declination T= matrix defining the transformation: new_rect= T dot old_rect, where old_rect = [cos(dec)cos(ra),cos(dec)sin(ra),sin(dec)] and similar for new_rect degree= (False) if True, ra and dec are given in degrees (default=False) OUTPUT: (pmphi1 x cos(phi2),pmph2) for vector inputs [:,2] HISTORY: 2016-10-24 - Written - Bovy (UofT/CCA) 2019-03-09 - uses custom_to_radec - Nathaniel Starkman (UofT) """ if T is None: raise ValueError("Must set T= for pmrapmdec_to_custom") # Need to figure out ra_ngp and dec_ngp for this custom set of sky coords ra_ngp, dec_ngp= custom_to_radec(0., numpy.pi/2, T=T) #Whether to use degrees and scalar input is handled by decorators dec[dec == dec_ngp]+= 10.**-16 #deal w/ pole. sindec_ngp= numpy.sin(dec_ngp) cosdec_ngp= numpy.cos(dec_ngp) sindec= numpy.sin(dec) cosdec= numpy.cos(dec) sinrarangp= numpy.sin(ra-ra_ngp) cosrarangp= numpy.cos(ra-ra_ngp) cosphi= sindec_ngp*cosdec-cosdec_ngp*sindec*cosrarangp sinphi= sinrarangp*cosdec_ngp norm= numpy.sqrt(cosphi**2.+sinphi**2.) cosphi/= norm sinphi/= norm return (numpy.array([[cosphi,-sinphi],[sinphi,cosphi]]).T\ *numpy.array([[pmra,pmra],[pmdec,pmdec]]).T).sum(-1)
[docs]def custom_to_radec(phi1,phi2,T=None,degree=False): """ NAME: custom_to_radec PURPOSE: rotate a custom set of sky coordinates (phi1, phi2) to (ra, dec) given the rotation matrix T for (ra, dec) -> (phi1, phi2) INPUT: phi1 - custom sky coord phi2 - custom sky coord T - matrix defining the transformation (ra, dec) -> (phi1, phi2) degree - default: False. If True, phi1 and phi2 in degrees OUTPUT: (ra, dec) for vector inputs [:, 2] HISTORY: 2018-10-23 - Written - Nathaniel (UofT) """ if T is None: raise ValueError("Must set T= for custom_to_radec") return radec_to_custom(phi1, phi2, T=numpy.transpose(T), # T.T = inv(T) degree=degree)
[docs]def custom_to_pmrapmdec(pmphi1,pmphi2,phi1,phi2,T=None,degree=False): """ NAME: custom_to_pmrapmdec PURPOSE: rotate proper motions in a custom set of sky coordinates (phi1,phi2) to ICRS (ra,dec) INPUT: pmphi1 - proper motion in custom (multplied with cos(phi2)) [mas/yr] pmphi2 - proper motion in phi2 [mas/yr] phi1 - custom longitude phi2 - custom latitude T= matrix defining the transformation in cartesian coordinates: new_rect = T dot old_rect where old_rect = [cos(dec)cos(ra), cos(dec)sin(ra), sin(dec)] and similar for new_rect degree= (False) if True, phi1 and phi2 are given in degrees (default=False) OUTPUT: (pmra x cos(dec), dec) for vector inputs [:,2] HISTORY: 2019-03-02 - Written - Nathaniel Starkman (UofT) """ if T is None: raise ValueError("Must set T= for custom_to_pmrapmdec") return pmrapmdec_to_custom(pmphi1, pmphi2, phi1, phi2, T=numpy.transpose(T), # T.T = inv(T) degree=degree)
def get_epoch_angles(epoch=2000.0): """ NAME: get_epoch_angles PURPOSE: get the angles relevant for the transformation from ra, dec to l,b for the given epoch INPUT: epoch - epoch of ra,dec (right now only 2000.0 and 1950.0 are supported when not using astropy's transformations internally; when internally using astropy's coordinate transformations, epoch can be None for ICRS, 'JXXXX' for FK5, and 'BXXXX' for FK4 [but for B1950 FK4 with no E abberation terms is assumed... really, there's no reason to use B1950 in 2018 when using galpy...)) OUTPUT: set of angles HISTORY: 2010-04-07 - Written - Bovy (NYU) 2016-05-13 - Added support for using astropy's coordinate transformations and for non-standard epochs - Bovy (UofT) 2018-04-18 - Edited J2000 angles to be fully consistent with astropy - BOvy (UofT) """ if epoch == 2000.0: # Following astropy's definition here theta= 122.9319185680026/180.*numpy.pi dec_ngp= 27.12825118085622/180.*numpy.pi ra_ngp= 192.8594812065348/180.*numpy.pi elif epoch == 1950.0: theta= 123./180.*numpy.pi dec_ngp= 27.4/180.*numpy.pi ra_ngp= 192.25/180.*numpy.pi elif epoch == None: # obtained below theta= theta_icrs dec_ngp= dec_ngp_icrs ra_ngp= ra_ngp_icrs elif _APY_LOADED: # Use astropy to get the angles epoch, frame= _parse_epoch_frame_apy(epoch) c= apycoords.SkyCoord(180.*units.deg,90.*units.deg, frame=frame,equinox=epoch) c= c.transform_to(apycoords.Galactic) theta= c.l.to(units.rad).value c= apycoords.SkyCoord(180.*units.deg,90.*units.deg, frame='galactic') if not epoch is None and 'J' in epoch: c= c.transform_to(apycoords.FK5(equinox=epoch)) elif not epoch is None and 'B' in epoch: c= c.transform_to(apycoords.FK4(equinox=epoch)) else: # pragma: no cover raise ValueError('epoch input not understood; should be None for ICRS, JXXXX, or BXXXX') dec_ngp= c.dec.to(units.rad).value ra_ngp= c.ra.to(units.rad).value else: raise IOError("Only epochs 1950 and 2000 are supported if you don't have astropy") return (theta,dec_ngp,ra_ngp) # Get ICRS angles once when astropy is installed if _APY_LOADED: c= apycoords.SkyCoord(180.*units.deg,90.*units.deg,frame='icrs') c= c.transform_to(apycoords.Galactic) theta_icrs= c.l.to(units.rad).value c= apycoords.SkyCoord(180.*units.deg,90.*units.deg, frame='galactic') c= c.transform_to(apycoords.ICRS) dec_ngp_icrs= c.dec.to(units.rad).value ra_ngp_icrs= c.ra.to(units.rad).value else: theta_icrs= 2.1455668515225916 dec_ngp_icrs= 0.4734773249532947 ra_ngp_icrs= 3.366032882941063 def _parse_epoch_frame_apy(epoch): if epoch == 2000.0 or epoch == '2000': epoch= 'J2000' elif epoch == 1950.0 or epoch == '1950': epoch= 'B1950' if not epoch is None and 'J' in epoch: frame= 'fk5' elif not epoch is None and 'B' in epoch: frame= 'fk4' else: frame= 'icrs' return (epoch,frame) # Matrix to rotate to the astropy Galactocentric frame: astropy's # Galactocentric frame is slightly off from the one that we get by simply # taking Galactic coordinates and transforming them: transformation from # Bovy (2011) maps NGP --> (0,0,1), astropy to (v. small, v. small, 1-v. small) # so we rotate Bovy (2011) such that we agree; for that we compute what NGP # goes to using astropy's transformations theta,dec_ngp,ra_ngp= get_epoch_angles(None) # None = ICRS, basis for astropy dec_gc,ra_gc= -28.936175/180.*numpy.pi,266.4051/180.*numpy.pi # from apy def. eta= 58.5986320306/180.*numpy.pi # astropy 'roll' angle gc_vec= numpy.array(\ [numpy.cos(theta)*(-numpy.sin(dec_ngp)*numpy.cos(dec_gc)*numpy.cos(ra_gc-ra_ngp) +numpy.cos(dec_ngp)*numpy.sin(dec_gc)) +numpy.sin(theta)*numpy.cos(dec_gc)*numpy.sin(ra_gc-ra_ngp), numpy.sin(theta)*(-numpy.sin(dec_ngp)*numpy.cos(dec_gc)*numpy.cos(ra_gc-ra_ngp) +numpy.cos(dec_ngp)*numpy.sin(dec_gc)) -numpy.cos(theta)*numpy.cos(dec_gc)*numpy.sin(ra_gc-ra_ngp), numpy.cos(dec_ngp)*numpy.cos(dec_gc)*numpy.cos(ra_gc-ra_ngp) +numpy.sin(dec_ngp)*numpy.sin(dec_gc)]) galcen_extra_rot1= _rotate_to_arbitrary_vector(numpy.atleast_2d(gc_vec), numpy.array([1.,0.,0.]), inv=False,_dontcutsmall=True)[0] ngp_vec= numpy.dot(galcen_extra_rot1,numpy.array(\ [-numpy.cos(dec_gc)*numpy.cos(dec_ngp)*numpy.cos(ra_ngp-ra_gc) -numpy.sin(dec_gc)*numpy.sin(dec_ngp), numpy.cos(eta)*numpy.cos(dec_ngp)*numpy.sin(ra_ngp-ra_gc) +numpy.sin(eta)*(-numpy.sin(dec_gc)*numpy.cos(dec_ngp)*numpy.cos(ra_ngp-ra_gc) +numpy.cos(dec_gc)*numpy.sin(dec_ngp)), -numpy.sin(eta)*numpy.cos(dec_ngp)*numpy.sin(ra_ngp-ra_gc) +numpy.cos(eta)*(-numpy.sin(dec_gc)*numpy.cos(dec_ngp)*numpy.cos(ra_ngp-ra_gc) +numpy.cos(dec_gc)*numpy.sin(dec_ngp))])) galcen_extra_rot2= _rotate_to_arbitrary_vector(numpy.atleast_2d(ngp_vec), numpy.array([0.,0.,1.]), inv=True,_dontcutsmall=True)[0] # Leave x axis alone, because in place by rot1 galcen_extra_rot2[0,0]= 1. galcen_extra_rot2[0,1:]= 0. galcen_extra_rot2[1:,0]= 0. galcen_extra_rot= numpy.dot(galcen_extra_rot2,galcen_extra_rot1)