dame_CO_with_midplanes.png

iau_midplane_but_Sun_is_in_plane.png

Model2.png

sgraframe.py

from astropy import coordinates
from astropy import units as u
from astropy.coordinates import (BaseCoordinateFrame, frame_transform_graph,
                                 RepresentationMapping,
                                 SphericalRepresentation,
                                 CartesianRepresentation)
from astropy.coordinates.baseframe import FrameAttribute
from astropy.coordinates.angles import rotation_matrix
from astropy import log
import numpy as np
import itertools

class GalacticPlane(BaseCoordinateFrame):
    """
    Galactocentric Coordinates.  Only Cartesian representation is directly supported.

    Parameters
    ----------
    representation : `BaseRepresentation` or None
        A representation object or None to have no data (or use the other keywords)
    distance : `~astropy.units.Quantity`, optional, must be keyword
        The Distance for this object along the line-of-sight.

    Example
    -------
    >>> x = coordinates.SkyCoord(x=-1*u.kpc, y=1*u.kpc, z=0.2*u.kpc,
    ...                          representation=coordinates.CartesianRepresentation)
    >>> x.transform_to(coordinates.Galactic)
    <SkyCoord (Galactic): l=220.780472997 deg, b=32.1851280336 deg, distance=1.42828568571 kpc>
    """

    default_representation = CartesianRepresentation

    gc_coord = FrameAttribute(coordinates.SkyCoord(l=359.944304966*u.deg,
                                                   b=-0.0461861959744*u.deg,
                                                   distance=8.34*u.kpc,
                                                   frame='galactic'))

    zsun =FrameAttribute(25*u.pc)
    



@frame_transform_graph.transform(coordinates.FunctionTransform,
                                 coordinates.Galactic, GalacticPlane)
def gal_to_gc(gal_coord, gc_frame):

    gal_recen = coordinates.SkyCoord(gal_coord.l-gc_frame.gc_coord.galactic.l,
                                     gal_coord.b-gc_frame.gc_coord.galactic.b,
                                     distance=gal_coord.distance)
    galcar = gal_recen.represent_as(coordinates.CartesianRepresentation)
    x,y,z = galcar.x, galcar.y, galcar.z

    # Determine x,y,z in galactic coordinates transformed to  gc coordinates (rotate coords around Sgr A*)
    gc_x, gc_y, gc_z = (1,0,0)*gc_frame.gc_coord.galactic.distance

    # rotate around y axis, centered on Sgr A*, to change x, z
    # x axis defined to go through sgr a and sun-Z_sun
    # z axis defined to be 'up', parallel to z_sun in Galactocentric and
    # slightly offset in Galactic coords
    gc_angle = np.arcsin(gc_frame.zsun/gc_frame.gc_coord.distance)
    rot = rotation_matrix(-gc_angle, 'y')
    rotated = np.dot(rot,
                     np.array([a.value for a in [x-gc_x,
                                                 y-gc_y,
                                                 z-gc_z]]))
    newx,newy,newz = np.array(rotated).squeeze()*x.unit

    newframe = GalacticPlane(newx, newy, newz, representation=coordinates.CartesianRepresentation,
                             gc_coord=gc_frame.gc_coord,
                             zsun=gc_frame.zsun)
    return newframe

@frame_transform_graph.transform(coordinates.FunctionTransform,
                                 GalacticPlane, coordinates.Galactic)
def gc_to_gal(gc_coord, gal_frame):

    gc_angle = np.arcsin(gc_coord.zsun/gc_coord.gc_coord.distance)
    rot = rotation_matrix(gc_angle, 'y')
    gc_x, gc_y, gc_z = (1,0,0)*gc_coord.gc_coord.galactic.distance

    xs,ys,zs = gc_coord.x, gc_coord.y, gc_coord.z

    xo,yo,zo = np.array(np.dot(rot, np.array([a.value
                                              for a in [xs,ys,zs]]))).squeeze()*xs.unit

    xo += gc_x
    yo += gc_y
    zo += gc_z

    newgal = coordinates.Galactic(xo, yo, zo, representation=coordinates.CartesianRepresentation)
    newskygal = coordinates.Galactic(newgal.spherical.lon+gc_coord.gc_coord.galactic.l,
                                     newgal.spherical.lat+gc_coord.gc_coord.galactic.b,
                                     distance=newgal.spherical.distance)

    return newskygal

def test_roundtrip():
    x = coordinates.SkyCoord(0.0*u.deg, 0.0*u.deg, distance=8.5*u.kpc, frame='galactic')
    y = x.transform_to(GalacticPlane).transform_to(coordinates.Galactic)
    np.testing.assert_almost_equal(x.l, y.l)
    np.testing.assert_almost_equal(x.b, y.b)

    x = coordinates.SkyCoord(5*u.deg, 12*u.deg, distance=3*u.kpc, frame='galactic')
    y = x.transform_to(GalacticPlane).transform_to(coordinates.Galactic)
    np.testing.assert_almost_equal(x.l, y.l)
    np.testing.assert_almost_equal(x.b, y.b)


def test_4kpc():
    c = coordinates.SkyCoord(-4*u.kpc, 0*u.kpc, 0*u.kpc, frame=GalacticPlane)
    log.info("Model 2: -4,0,0 in gc Galactic midplane: {0}".format(c.transform_to('galactic')))
    c = coordinates.SkyCoord(-4*u.kpc, 0*u.kpc, 0*u.kpc, frame=GalacticPlane,
                             zsun=0*u.pc)
    cgal = coordinates.SkyCoord(gc_to_gal(c.frame, coordinates.Galactic), frame='galactic')
    log.info("Model X: -4,0,0 in Sun-on-gc plane, correcting for IAU GC: {0}".format(cgal))
    c = coordinates.SkyCoord(-4*u.kpc, 0*u.kpc, 0*u.kpc, frame=GalacticPlane,
                             gc_coord=coordinates.SkyCoord(0*u.deg, 0*u.deg,
                                                           distance=8.5*u.kpc,
                                                           frame='galactic'))
    cgal = coordinates.SkyCoord(gc_to_gal(c.frame, coordinates.Galactic), frame='galactic')
    log.info("Model 3(?): -4,0,0 in IAU Galactic midplane: {0}".format(cgal))

def test_nessie():
    """ Try to reproduce https://www.authorea.com/users/23/articles/249/master/file/figures/sideview_both/sideview_both.jpg """
    log.info("Goodman paper: 'no tilt' = -0.48, 'tilt' = -0.36")
    xc = -(8.34-2.96)*u.kpc
    yc = -0.906*u.kpc
    c = coordinates.SkyCoord(xc, yc, 0*u.kpc, frame=GalacticPlane)
    log.info("Model 2: nessie in gc Galactic midplane: {0}".format(c.transform_to('galactic')))
    c = coordinates.SkyCoord(xc, yc, 0*u.kpc, frame=GalacticPlane, zsun=0*u.pc)
    cgal = coordinates.SkyCoord(gc_to_gal(c.frame, coordinates.Galactic), frame='galactic')
    log.info("Model X: nessie in Sun-on-gc plane, correcting for IAU GC: {0}".format(cgal))
    c = coordinates.SkyCoord(xc, yc, 0*u.kpc, frame=GalacticPlane,
                             gc_coord=coordinates.SkyCoord(0*u.deg, 0*u.deg,
                                                           distance=8.5*u.kpc,
                                                           frame='galactic'))
    cgal = coordinates.SkyCoord(gc_to_gal(c.frame, coordinates.Galactic), frame='galactic')
    log.info("Model 3(?): nessie in IAU Galactic midplane: {0}".format(cgal))
    c = coordinates.SkyCoord(xc, yc, 0*u.kpc, frame=GalacticPlane,
                             gc_coord=coordinates.SkyCoord(0*u.deg,
                                                           np.arcsin(25*u.pc/(8.5*u.kpc)),
                                                           distance=8.5*u.kpc,
                                                           frame='galactic'))
    cgal = coordinates.SkyCoord(gc_to_gal(c.frame, coordinates.Galactic), frame='galactic')
    log.info("Model 1(?): nessie in IAU Galactic midplane: {0}".format(cgal))


def w51_test():
    # Find the x,y coordinates approximately corresponding to W51
    w51fil = coordinates.SkyCoord(np.linspace(48,52)*u.deg, np.ones(50)*-0.3*u.deg, distance=5.1*u.kpc, frame='galactic')
    w51fil_plane = w51fil.transform_to(GalacticPlane)

    # Use x,y but replace z with zero to get midplane
    midplane = coordinates.SkyCoord(w51fil_plane.x, w51fil_plane.y, np.zeros(50)*u.kpc, frame=GalacticPlane)

    # Convert back to glon
    midplane_g = midplane.transform_to('galactic')
    return midplane_g

def full_midplane(distance=3*u.kpc, npix=500, bounds=[10,60], transform_frame=GalacticPlane):
    # Find the x,y coordinates approximately corresponding to the specified distance
    posplane = coordinates.SkyCoord(np.linspace(bounds[0], bounds[1],
                                                npix)*u.deg,
                                    np.ones(npix)*0.0*u.deg, distance=distance,
                                    frame='galactic')
    #posplane_plane = gal_to_gc(posplane, transform_frame)
    posplane_plane = posplane.transform_to(transform_frame)

    # Use x,y but replace z with zero to get midplane
    midplane = coordinates.SkyCoord(posplane_plane.x, posplane_plane.y,
                                    np.zeros(npix)*u.kpc,
                                    frame=GalacticPlane,
                                    gc_coord=transform_frame.gc_coord,
                                    zsun=transform_frame.zsun)

    # Convert back to glon
    midplane_g = gc_to_gal(midplane.frame, coordinates.Galactic)
    return coordinates.SkyCoord(midplane_g)

def show_midplanes_on_dame(bounds=[10,60], distances=[1,3,5,10]*u.kpc,
                           color_cycle=itertools.cycle('rgbcmy'), height=16):
    from astropy.utils.data import download_file
    import aplpy
    filename = download_file('http://www.cfa.harvard.edu/mmw/Wco_DHT2001.fits.gz', cache=True)
    F3 = aplpy.FITSFigure(filename)
    F3.show_grayscale()

    for ii,distance in enumerate(distances):
        mp = full_midplane(distance=distance, bounds=bounds)

        color = color_cycle.next()
        F3.show_lines(np.array([[mp.l.deg, mp.b.deg]]), color=color)
        F3.add_label(min(bounds)+5, (height/2-1)-(ii*16./height), str(distance), color=color)

    F3.recenter(np.mean(bounds), 0, width=np.abs(np.diff(bounds)), height=height)
    F3.save('dame_CO_with_midplanes.png')

    return F3

def compare_tilted_untilted_midplanes(bounds=[333, 342],
                                      distances=[3.1]*u.kpc,
                                      color_cycle=itertools.cycle('rgbmk'),
                                      height=2):
    from astropy.utils.data import download_file
    import aplpy
    filename = download_file('http://www.cfa.harvard.edu/mmw/Wco_DHT2001.fits.gz', cache=True)
    F3 = aplpy.FITSFigure(filename)
    F3.show_grayscale()

    tilted = GalacticPlane()
    b0l0 = GalacticPlane(zsun=25*u.pc,
                         gc_coord=coordinates.SkyCoord(0*u.deg, 0*u.deg,
                                                              frame='galactic',
                                                              distance=8.5*u.kpc))
    z0pc = GalacticPlane(zsun=0*u.pc)
    model1 = GalacticPlane(gc_coord=coordinates.SkyCoord(0*u.deg,
                                                                -10.305011748157108*u.arcmin,
                                                                distance=8.5*u.kpc,
                                                                frame='galactic'))

    for ii,distance in enumerate(distances):
        mpt = full_midplane(distance=distance, bounds=bounds, transform_frame=tilted)
        mp1 = full_midplane(distance=distance, bounds=bounds, transform_frame=b0l0)
        mp2 = full_midplane(distance=distance, bounds=bounds, transform_frame=z0pc)
        mp3 = full_midplane(distance=distance, bounds=bounds, transform_frame=model1)
        log.debug("z=0pc dl: {0}".format((mpt[0].l-mp2[0].l).to(u.deg)))
        log.debug("z=0pc db: {0}".format((mpt[0].b-mp2[0].b).to(u.deg)))
        log.debug("l=0,b=0 dl: {0}".format((mpt[0].l-mp1[0].l).to(u.deg)))
        log.debug("l=0,b=0 db: {0}".format((mpt[0].b-mp1[0].b).to(u.deg)))

        color = color_cycle.next()
        F3.show_lines(np.array([[mpt.l.deg, mpt.b.deg]]), color=color)
        F3.add_label(min(bounds)+5, -0.1, "gc centric, 25pc rotation "+str(distance), color=color)

        color = color_cycle.next()
        F3.show_lines(np.array([[mp2.l.deg, mp2.b.deg]]), color=color)
        F3.add_label(min(bounds)+5, +0.3, "gc centric, no rotation "+str(distance), color=color)

        color = color_cycle.next()
        F3.show_lines(np.array([[mp1.l.deg, mp1.b.deg]]), color=color)
        F3.add_label(min(bounds)+5, +0.1, "l=0,b=0 centric, 25 pc rotation "+str(distance), color=color)

        color = color_cycle.next()
        F3.show_lines(np.array([[mp3.l.deg, mp3.b.deg]]), color=color)
        F3.add_label(min(bounds)+5, -0.7, "l=0,b=0 centric, 25 pc rotation *twice* (Model 1)"+str(distance), color=color, weight='bold')

    F3.recenter(np.mean(bounds), 0, width=np.abs(np.diff(bounds)), height=height)
    F3.recenter(337.5, 0, width=9, height=2)
    F3.save('tilted_vs_untilted.png')


    return F3

def model_diagram():
    import pylab as pl
    fig1 = pl.figure(1)
    fig1.clf()
    ax = fig1.gca()
    ax.plot([-1,1], [0,0], 'k--', label='True Midplane')
    ax.plot([-1],[0.25], 'yo', label='Sun')
    ax.plot([-1,0.5], [0.25, 0], 'b-', label='IAU midplane')
    ax.plot([0,0],[0,0], 'r*', label='Sgr A*')
    ax.plot([-0.5], [0], 'gs', label='Object in midplane')
    ax.plot([-1,-0.5], [0.25,0], 'g-', label='Latitude of Midplane d=4kpc')
    ax.text(-0.8,0.15, '$\\theta_1$')
    ax.axis([-1.25,1.25,-0.25,0.5])
    ax.legend(loc='best')
    fig1.savefig('theta1.png')

    fig2 = pl.figure(2)
    fig2.clf()
    ax = fig2.gca()
    ax.plot([-1,1], [0,0], 'k--', label='True Midplane')
    ax.plot([-1],[0.25], 'yo', label='Sun')
    ax.plot([-1,-1],[0.25,0], 'y-', linewidth=2, zorder=-1, alpha=0.8)
    ax.plot([-1,0.5], [0.25, 0], 'b-', label='IAU midplane')
    ax.plot([0],[1/12.],'bs',label='IAU b=0')
    ax.plot([-1,0.0], [0.25, 0], 'b--', label='Redefined, nonphysical midplane through Sgr A*')
    ax.plot([0,0],[0,0], 'r*', label='Sgr A*')
    ax.plot([-0.5], [0], 'gs', label='Object in midplane')
    ax.plot([-0.5,0], [0,0], 'g-', linewidth=2, zorder=-1, alpha=0.5)
    ax.plot([-1,-0.5], [0.25,0], 'g-', label='Latitude of Midplane d=4kpc')
    ax.axis([-1.25,1.25,-0.25,0.5])
    ax.legend(loc='best')
    fig2.savefig('Model2.png')

    ax.text(-1.1, 0.125, '$z_{\\odot}$')
    ax.text(-0.25,-0.02,'$x_s$', color='g')
    ax.text(-0.79,0.1,'$d_s$', color='g')
    ax.text(-0.8,0.15, '$\\theta_2$', label='$\\theta_2 < \\theta_1$, ')
    ax.text(-0.2, 0.01, '$\phi_{mp}$')
    ax.text(-0.5, -0.12, r'$d_s = \sqrt{z_\odot^2 + (d_{GC}-x_s)^2}$')
    ax.text(-0.5, -0.15, r'$\theta_2 = \sin^{-1}(d_{GC}/z_\odot) - \sin^{-1}((d_{GC}-x_s)/z_\odot)$')
    ax.text(-0.5, -0.18, r"$\phi=\phi'$ (longitude is unaffected)")
    fig2.savefig('theta2.png')

    fig3 = pl.figure(3)
    fig3.clf()
    ax = fig3.gca()
    ax.plot([-1,1], [0,0], 'k--', label='True Midplane')
    ax.plot([-1],[0.0], 'yo', label='Sun (as if it were in the midplane)')
    ax.plot([-1,0.0], [0.0, 1/12.], 'b-', label='IAU midplane')
    ax.plot([0.0], [1/12.], 'bs')
    ax.plot([0,0],[0,0], 'r*', label='Sgr A*')
    ax.plot([-0.5], [0], 'gs', label='Object in midplane')
    #ax.text(-0.8,0.15, '$\\theta_1$')
    ax.axis([-1.25,1.25,-0.25,0.5])
    ax.set_title("Everything in the plane has the same latitude")
    ax.legend(loc='best')
    fig3.savefig('iau_midplane_but_Sun_is_in_plane.png')

    fig4 = pl.figure(4)
    fig4.clf()
    ax = fig4.gca()

    ax.plot(0,0,'ks')
    ax.plot([-1,1],[0,0],'k--')
    ax.plot([0,0],[-1,1],'k--')
    ax.text(0.01, 0.01, '$0,0$')
    ax.plot(0.25,-0.5, 'bo')
    ax.plot([0.25,0.25], [-1,1], 'b--')
    ax.plot([-1,1], [-0.5,-0.5], 'b--')
    ax.text(0.26, -0.51, '$a,b$', color='b')
    # y=mx+b
    # y=2x-1
    # y=-1/2 x - 0.375
    ax.plot([-1,1],[-3,1],'r--')
    ax.plot([-1,1],[0.125,-0.875],'r--')
    ax.text(0.4,-0.55,r'$\theta$', color='r')

    ax.plot(0.2,0.3,'g*')
    ax.text(0.26,0.31,r'$x,y \rightarrow$',color='k')
    ax.text(0.26,0.25,r'$x+a,y+b \rightarrow$',color='b')
    ax.text(0.26,0.19,r'$(x+a) \cos(\theta) + (y+b) \sin(\theta),$',color='r')
    ax.text(0.26,0.13,r'$ (x+a)(-\sin(\theta)) + (y+b)\cos(\theta)$',color='r')
    ax.axis([-1,1,-1,1])
    ax.set_aspect('equal')

    """
    Notes to self, for discarding...
    zsun = gc_coord._sun_above_midplane_distance
    dgc = gc_coord.gc_coord.galactic.distance
    phi_mp = gc_coord._gc_angle

    distance = (zsun**2 + (dgc - xs)**2)**0.5
    angle_from_below = np.arcsin((dgc-xs)/zsun)
    latitude = phi_mp - angle_from_below
    longitude = ()
    """

tilted_vs_untilted.png

Compare to https://www.authorea.com/users/23/articles/249/master/file/figures/sideview_both/sideview_both.jpg and
https://www.authorea.com/users/23/articles/249/_show_article

@keflavich - One side item - it looks to me like you're using a generator function to make a new frame. The way I (and @adrn) were thinking, instead there's just one Galactocentric frame, but it takes frame attributes that specify all the parameters like coordinates of Sgr A* and such. They would default to sensible values that more-or-less reproduce Galactic (maybe with the fixed Sgr A* location?), but in principal all of them are user-selectable via the FrameAttribute machinery. Did you implement it as a generator function for some reason I'm missing, or just for convinience/you're not familiar enough with frame attributes?

The latter: I've never heard of FrameAttributes before now.

@eteq - this is to ping you in case you get notifications