plot.py

import gala.potential as gp
import numpy as np
from idlplotInd import tvhist2d
import astropy.table as atpy
import astropy.io.fits as pyfits
import astropy.units as auni
import astropy.coordinates as acoo
import matplotlib.pyplot as plt
import sys

def doplot():
  
    rmax = 40
    rmin = -20
 
    T = atpy.Table().read('mws_figure_data1.fits', mask_invalid=True)
    DD = pyfits.getdata('mws_figure_data2.fits')
    plt.figure(1, figsize=(3.37, 3.37))
    tvhist2d(
        -T['x'],  
        T['z'],
        rmin,
        rmax,
        -10,
        10,
        weights=T['feh'],
        statistic='median',
        cmap='turbo',
        vmin=-1.5,
        bar=True,
        vmax=0.2,
        bins=[100, 100])
    if True:
        plt.contour(np.log10(DD),
                    extent=[-50, 50, -50, 50],
                    levels=[3, 4, 5],
                    colors='orange',
                    alpha=0.5)
    plt.axhline(0, color='orange')
    aind = (np.abs(T['z']) > 10)

    plt.scatter(-T['x'][aind],
                T['z'][aind],
                c=T['feh'][aind],
                s=3,
                alpha=.3,
                cmap='turbo',
                norm=plt.Normalize(vmin=-1.5, vmax=.2))
    plt.xlim(-50, 50)
    plt.ylim(-50, 50)
    plt.gcf().set_size_inches(3.37 * 2, 3.37 * 2.)
    plt.xlabel('X [kpc]')
    plt.ylabel('Z [kpc]')
    plt.tight_layout()
    plt.savefig('fig_mws.pdf')