In [1]:
import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import fitsio
import astropy.units as u
from astropy.io import ascii
from astropy.table import Table
from astropy.cosmology import FlatLambdaCDM
In [2]:
%pylab inline
In [3]:
mpl.rcParams.update({'font.size': 18})
cosmo = FlatLambdaCDM(H0=70, Om0=0.3)
In [4]:
massdir = os.path.join( os.getenv('HIZEA_PROJECT'), 'massprofiles', 'isedfit' )
etcdir = os.path.join( os.getenv('HIZEA_DIR'), 'etc' )
In [5]:
photfile = os.path.join(etcdir, 'sg_fluxtable_nm.txt')
isedfile = os.path.join(massdir, 'massprofiles_fsps_v2.4_miles_chab_charlot_sfhgrid01.fits.gz')
kcorrfile = os.path.join(massdir, 'massprofiles_fsps_v2.4_miles_chab_charlot_sfhgrid01_kcorr.z0.0.fits.gz')
In [6]:
print('Reading {}'.format(photfile))
phot = ascii.read(photfile)
phot[:2]
Out[6]:
In [7]:
print('Reading {}'.format(isedfile))
ised = Table(fitsio.read(isedfile, ext=1, upper=True))
ised[:2]
Out[7]:
In [8]:
print('Reading {}'.format(kcorrfile))
kcorr = Table(fitsio.read(kcorrfile, ext=1, upper=True))
kcorr[:2]
Out[8]:
In [9]:
galaxy = [gg[:5] for gg in phot['ID'].data]
galaxy = np.unique(galaxy)
ngal = len(galaxy)
In [10]:
nrad = 40
radpix = np.linspace(1.0, 40.0, nrad) # [pixels]
radarcsec = radpix * 0.05 # [arcsec]
In [11]:
mstar = ised['MSTAR_AVG'].data.reshape(ngal, nrad)
mstar_err = ised['MSTAR_ERR'].data.reshape(ngal, nrad)
redshift = phot['z'].data.reshape(ngal, nrad)[:, 0]
In [12]:
area = np.pi * np.insert(np.diff(radarcsec**2), 0, radarcsec[0]**2) # aperture annulus [arcsec2]
In [13]:
sigma = np.zeros_like(mstar) # surface mass density [Mstar/kpc2]
radkpc = np.zeros_like(mstar) # radius [comoving kpc]
In [14]:
for igal in range(ngal):
arcsec2kpc = cosmo.arcsec_per_kpc_comoving(redshift[igal]).value
radkpc[igal, :] = radarcsec / arcsec2kpc
areakpc2 = area / arcsec2kpc**2
sigma[igal, :] = np.log10( 10**mstar[igal, :] / areakpc2 )
In [15]:
massrange = (8, 10.2)
sigmarange = (6, 9.6)
In [16]:
fig, ax = plt.subplots(3, 4, figsize=(14, 8), sharey=True, sharex=True)
for ii, thisax in enumerate(ax.flat):
thisax.errorbar(radarcsec, mstar[ii, :], yerr=mstar_err[ii, :],
label=galaxy[ii])
thisax.set_ylim(massrange)
#thisax.legend(loc='upper right', frameon=False)
thisax.annotate(galaxy[ii], xy=(0.9, 0.9), xycoords='axes fraction',
size=16, ha='right', va='top')
fig.text(0.0, 0.5, r'$\log_{10}\, (M / M_{\odot})$', ha='center',
va='center', rotation='vertical')
fig.text(0.5, 0.0, 'Radius (arcsec)', ha='center',
va='center')
fig.subplots_adjust(wspace=0.05, hspace=0.05)
fig.tight_layout()
In [17]:
fig, ax = plt.subplots(figsize=(10, 7))
for igal in range(ngal):
ax.plot(radkpc[igal, :], np.log10(np.cumsum(10**mstar[igal, :])), label=galaxy[igal])
ax.legend(loc='lower right', fontsize=16, ncol=3, frameon=False)
ax.set_xlabel(r'Galactocentric Radius $r_{kpc}$ (Comoving kpc)')
ax.set_ylabel(r'$\log_{10}\, M(<r_{kpc})\ (M_{\odot})$')
Out[17]:
In [18]:
fig, ax = plt.subplots(figsize=(10, 7))
for igal in range(ngal):
ax.plot(radkpc[igal, :], sigma[igal, :], label=galaxy[igal])
ax.legend(loc='upper right', fontsize=16, ncol=3, frameon=False)
ax.set_xlabel(r'Galactocentric Radius $r_{kpc}$ (Comoving kpc)')
ax.set_ylabel(r'$\log_{10}\, \Sigma\ (M_{\odot}\ /\ {\rm kpc}^2)$')
ax.set_ylim(sigmarange)
Out[18]:
In [ ]: