In [3]:
%load_ext autoreload
%reload_ext autoreload
%autoreload 2
% matplotlib inline
from __future__ import (division,
print_function)
import os
import sys
import copy
import fnmatch
import warnings
import collections
import numpy as np
import scipy
try:
from scipy.stats import scoreatpercentile
except:
scoreatpercentile = False
from scipy.interpolate import interp1d
import cPickle as pickle
from scipy import ndimage
# Astropy
from astropy.io import fits
from astropy import units as u
from astropy.stats import sigma_clip
from astropy.table import Table, Column
from astropy.utils.console import ProgressBar
from astropy.convolution import convolve, Box1DKernel
# AstroML
from astroML.plotting import hist
from astroML.density_estimation import KNeighborsDensity
try:
from sklearn.neighbors import KernelDensity
use_sklearn_KDE = True
except:
import warnings
warnings.warn("KDE will be removed in astroML version 0.3. Please "
"upgrade to scikit-learn 0.14+ and use "
"sklearn.neighbors.KernelDensity.", DeprecationWarning)
from astroML.density_estimation import KDE
use_sklearn_KDE = False
from sklearn.neighbors import KDTree
from sklearn.neighbors import BallTree
# Matplotlib related
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.style.use('classic')
import matplotlib.patches as mpatches
import matplotlib.gridspec as gridspec
from matplotlib.patches import Ellipse
from matplotlib.ticker import NullFormatter, MaxNLocator, FormatStrFormatter
from matplotlib.collections import PatchCollection
tickFormat = FormatStrFormatter('$\mathbf{%g}$')
# Personal
import hscUtils as hUtil
import galSBP
import coaddCutoutGalfitSimple as gSimple
# Cosmology
import cosmology
c=cosmology.Cosmo(H0=70.0, omega_m=0.3, omega_l=0.7, flat=1)
# SEP
import sep
# Color map
from palettable.colorbrewer.sequential import Greys_9, OrRd_9, Blues_9, Purples_9, YlGn_9
BLK = Greys_9.mpl_colormap
ORG = OrRd_9.mpl_colormap
BLU = Blues_9.mpl_colormap
GRN = YlGn_9.mpl_colormap
PUR = Purples_9.mpl_colormap
# Personal tools
import galSBP as sbp
from hscUtils import songPlotSetup, removeIsNullCol
from hscUtils import confidence_interval, ma_confidence_interval_1d, confidence_interval_1d
## Constants
# SDSS pivot wavelength
sdss_u_pivot = 3551.0
sdss_g_pivot = 4686.0
sdss_r_pivot = 6165.0
sdss_i_pivot = 7481.0
sdss_z_pivot = 8931.0
# HSC pivot wavelength
hsc_g_pivot = 4782.2
hsc_r_pivot = 6101.7
hsc_i_pivot = 7648.0
hsc_z_pivot = 8883.0
hsc_y_pivot = 9750.8
hscFiltWave = np.asarray([hsc_g_pivot, hsc_r_pivot, hsc_i_pivot, hsc_z_pivot, hsc_y_pivot])
"""
Absolute magnitude of the Sun in HSC filters
Right now, just use the DES filters
"""
SUN_G = 5.08
SUN_R = 4.62
SUN_I = 4.52
SUN_Z = 4.52
SUN_Y = 4.51
# Solar stellar metallicity
Z_SUN = 0.02
# definitions for the axes
left, width = 0.15, 0.64
right = left + width
bottom, height = 0.13, 0.86
bottom_h = left_h = left + width + 0.02
recScat = [left, bottom, width, height]
recHist = [right, bottom, 0.20, height]
SBP1 = [0.124, 0.085, 0.865, 0.33]
SBP2 = [0.124, 0.41, 0.865, 0.55]
EC1 = [0.135, 0.072, 0.862, 0.295]
EC2 = [0.135, 0.366, 0.862, 0.295]
EC3 = [0.135, 0.666, 0.862, 0.295]
REC = [0.12, 0.11, 0.87, 0.87]
COG1 = [0.143, 0.10, 0.850, 0.43]
COG2 = [0.143, 0.53, 0.850, 0.43]
# Universal RSMA array
RSMA_COMMON = np.arange(0.4, 4.2, 0.01)
RR50_COMMON = np.arange(0.0, 9.0, 0.01)
EMPTY = (RSMA_COMMON * np.nan)
# Color
BLUE0 = "#92c5de"
BLUE1 = "#0571b0"
RED0 = "#f4a582"
RED1 = "#ca0020"
PURPLE0 = '#af8dc3'
PURPLE1 = '#762a83'
BROWN0 = '#bf812d'
BROWN1 = '#543005'
GREEN0 = '#7fbf7b'
GREEN1 = '#1b7837'
# 3-sigma
SIGMA1 = 0.3173
SIGMA2 = 0.0455
SIGMA3 = 0.0027
# Cubehelix
from palettable.cubehelix import Cubehelix
CUBE = Cubehelix.make(start=0.3, rotation=-0.5, n=16).mpl_colormap
import itertools
In [48]:
def loadPkl(filename):
try:
import cPickle as pickle
except:
warnings.warn("## cPickle is not available!!")
import pickle
if os.path.isfile(filename):
pklFile = open(filename, 'rb')
data = pickle.load(pklFile)
pklFile.close()
return data
else:
warnings.warn("## Can not find %s, return None" % filename)
return None
In [4]:
imgDir = 'example'
In [5]:
ls example
In [6]:
#isophote = '/iraf/iraf/extern/stsdas/bin.macosx/x_isophote.e'
#xttools = '/iraf/iraf/extern/tables/bin.macosx/x_ttools.e'
isophote = '/Users/song/anaconda/pkgs/iraf-2.16.1-0/variants/common/iraf/stsci_iraf/stsdas/bin.macosx/x_isophote.e'
xttools = '/Users/song/anaconda/pkgs/iraf-2.16.1-0/variants/common/iraf/stsci_iraf/tables/bin.macosx/x_ttools.e'
# The Gaussian kernal used for convolution
kernel1 = np.asarray([[0.560000, 0.980000, 0.560000],
[0.980000, 1.000000, 0.980000],
[0.560000, 0.980000, 0.560000]])
kernel2 = np.asarray([[0.000000, 0.220000, 0.480000, 0.220000, 0.000000],
[0.220000, 0.990000, 1.000000, 0.990000, 0.220000],
[0.480000, 1.000000, 1.000000, 1.000000, 0.480000],
[0.220000, 0.990000, 1.000000, 0.990000, 0.220000],
[0.000000, 0.220000, 0.480000, 0.220000, 0.000000]])
kernel3 = np.asarray([[0.092163, 0.221178, 0.296069, 0.221178, 0.092163],
[0.221178, 0.530797, 0.710525, 0.530797, 0.221178],
[0.296069, 0.710525, 0.951108, 0.710525, 0.296069],
[0.221178, 0.530797, 0.710525, 0.530797, 0.221178],
[0.092163, 0.221178, 0.296069, 0.221178, 0.092163]])
In [5]:
len(cat)
Out[5]:
In [8]:
imgMs = fits.open('example/map_stellar_central.fits')[0].data
headMs = fits.open('example/map_stellar_central.fits')[0].header
imgMd = fits.open('example/map_stellar_halo.fits')[0].data
imgMg = fits.open('example/map_gas_halo.fits')[0].data
imgMh = fits.open('example/map_darkmatter_halo.fits')[0].data
In [26]:
print(imgMs.shape, imgMd.shape, imgMg.shape, imgMh.shape)
example = [imgMs, imgMd, imgMg, imgMh]
labels = [r'$\mathrm{Stellar\ Mass\ (Central)}$',
r'$\mathrm{Stellar\ Mass\ (Halo)}$',
r'$\mathrm{Gas\ Map}$',
r'$\mathrm{Dark\ Matter\ Map}$']
exampleFile = ['example/map_stellar_central.fits',
'example/map_stellar_halo.fits',
'example/map_gas_halo.fits',
'example/map_darkmatter_halo.fits']
print(len(example))
In [51]:
h0 = 0.7
scale = 1.86 / h0 # kpc
cenX, cenY = 81.0, 81.0
In [52]:
# Prepare the
fig = plt.figure(figsize=(14, 14))
fig.subplots_adjust(left=0.0, right=1.0,
bottom=0.0, top=1.0,
wspace=0.00, hspace=0.00)
gs = gridspec.GridSpec(2, 2)
gs.update(wspace=0.0, hspace=0.00)
xArr, yArr, baArr, paArr, m10Arr, m30Arr, m100Arr, mtotArr = [], [], [], [], [], [], [], []
m150Arr, m200Arr = [], []
aperProf, aperMass = [], []
r50Arr = []
with ProgressBar(len(example)) as bar:
for item in enumerate(example):
ii, img = item
bar.update()
imgN = img.byteswap().newbyteorder()
# Extract objects
bkg = sep.Background(imgN, bw=20, bh=20, fw=4, fh=4)
imgSub = imgN - bkg
objects = sep.extract(imgSub, 15.0, filter_kernel=kernel3)
# Find the object at the center
if len(objects) == 1:
index = 0
else:
index = np.argmin(np.sqrt((objects['x'] - cenX) ** 2.0 + (objects['y'] - cenY) ** 2.0))
# Get the naive ba, theta, xcen, ycen
ba = objects[index]['b'] / objects[index]['a']
theta = objects[index]['theta']
xcen, ycen = objects[index]['x'], objects[index]['y']
baArr.append(ba)
paArr.append(theta * 180.0 / np.pi)
xArr.append(xcen)
yArr.append(ycen)
# Aperture photometry and naive "profiles"
rad = np.linspace(1.2, 4.0, 50)
kpc = (rad ** 4.0)
pix = (kpc / scale)
area = (np.pi * (kpc ** 2.0) * ba)
maper = sep.sum_ellipse(imgN, xcen, ycen, pix, pix * ba, theta, 1.0,
bkgann=None, subpix=11)[0]
mring = maper[1:] - maper[0:-1]
aring = area[1:] - area[0:-1]
rho = np.log10(mring / aring)
aperMass.append(maper)
aperProf.append(rho)
# Mass within different apertures
rad = np.asarray([10.0, 30.0, 100.0, 150.0, 200.0]) / scale
maper = sep.sum_ellipse(imgN, xcen, ycen, rad, rad * ba, theta, 1.0,
bkgann=None, subpix=11)[0]
m10Arr.append(np.log10(maper[0]))
m30Arr.append(np.log10(maper[1]))
m100Arr.append(np.log10(maper[2]))
m150Arr.append(np.log10(maper[3]))
m200Arr.append(np.log10(maper[4]))
mtotArr.append(np.log10(np.sum(imgN)))
# Get r50
r50, flag = sep.flux_radius(imgN, xcen, ycen, 30.0, 0.5,
subpix=11)
r50Arr.append(r50 * scale)
# Show the aperture
ax = plt.subplot(gs[ii])
ax.yaxis.set_major_formatter(NullFormatter())
ax.xaxis.set_major_formatter(NullFormatter())
m, s = np.mean(imgN), np.std(imgN)
im = ax.imshow(imgN, interpolation='nearest', cmap=CUBE,
vmin=m-1.5*s, vmax=m+10*s, origin='lower')
# plot an ellipse for each object
for i in range(len(objects)):
e = Ellipse(xy=(xcen, ycen),
width=40 * 2.0,
height=(40 * 2.0 * ba),
angle=(theta * 180. / np.pi))
e.set_facecolor('none')
e.set_edgecolor('r')
e.set_alpha(0.5)
e.set_linewidth(1.0)
ax.add_artist(e)
ax.set_aspect('equal')
# Show the 2D and 3D mass within 100 kpc
ax.text(0.5, 0.90, labels[ii],
verticalalignment='center', horizontalalignment='center',
fontsize=40.0, transform=ax.transAxes, weight='bold',
color='r', alpha=0.8)
ax.text(0.5, 0.15, "2D: %5.2f" % np.log10(maper[3]),
verticalalignment='center', horizontalalignment='center',
fontsize=20.0, transform=ax.transAxes, weight='bold',
color='r', alpha=0.8)
fig.savefig('mblack2_map_aperture_example.pdf', dpi=90)
In [53]:
baArr[0], paArr[0]
Out[53]:
In [54]:
# Prepare the plot
fig = plt.figure(figsize=(14, 14))
fig.subplots_adjust(left=0.0, right=1.0,
bottom=0.0, top=1.0,
wspace=0.00, hspace=0.00)
gs = gridspec.GridSpec(2, 2)
gs.update(wspace=0.0, hspace=0.00)
xArr, yArr, ba1Arr, pa1Arr, ba2Arr, pa2Arr, ba3Arr, pa3Arr = [], [], [], [], [], [], [], []
smaProf, ellProf, paProf, sbpProf = [], [], [], []
with ProgressBar(len(exampleFile)) as bar:
for item in enumerate(exampleFile):
ii, name = item
img = example[ii]
bar.update()
# Convert the mass density map to just mass maps
imgSave = name.replace('.fits', '_conv.fits')
imgN = img.byteswap().newbyteorder()
imgC = ndimage.convolve(imgN, kernel2, mode='constant', cval=0.0)
hdu = fits.PrimaryHDU(imgC)
hdulist = fits.HDUList([hdu])
hdulist.writeto(imgSave, overwrite=True)
try:
ellOut, binOut = galSBP.galSBP(imgSave, galX=cenX, galY=cenY,
maxSma=120.0, iniSma=15.0,
verbose=False, savePng=False, saveOut=False,
useZscale=False, expTime=1.0,
pix=scale, zpPhoto=0.0,
galQ=baArr[0], galPA=paArr[0],
stage=2, minSma=0.5, ellipStep=0.08,
isophote=isophote, xttools=xttools,
uppClip=9.0, lowClip=9.0, fracBad=0.6,
harmonics="none", maxTry=9, olthresh=0.2,
nClip=0, intMode='median', updateIntens=False)
except Exception:
ellOut, binOut = None, None
# Show the aperture
ax = plt.subplot(gs[ii])
ax.yaxis.set_major_formatter(NullFormatter())
ax.xaxis.set_major_formatter(NullFormatter())
m, s = np.mean(imgC), np.std(imgC)
im = ax.imshow(imgC, interpolation='nearest', cmap=CUBE,
vmin=m-1.5*s, vmax=m+10*s, origin='lower')
if ellOut is not None:
# Estimate the average ellipticity and position angle
sma = ellOut['sma']
kpc = sma * scale
xArr.append(ellOut[1]['x0'])
yArr.append(ellOut[1]['y0'])
avgEll1 = np.nanmedian(ellOut['ell'][(kpc > 6.0) & (kpc < 30.0)])
avgPA1 = np.nanmedian(ellOut['pa'][(kpc > 6.0) & (kpc < 30.0)])
ba1Arr.append(1.0 - avgEll1)
pa1Arr.append(avgPA1)
avgEll2 = np.nanmedian(ellOut['ell'][(kpc > 30.0) & (kpc < 50.0)])
avgPA2 = np.nanmedian(ellOut['pa'][(kpc > 30.0) & (kpc < 50.0)])
ba2Arr.append(1.0 - avgEll2)
pa2Arr.append(avgPA2)
avgEll3 = np.nanmedian(ellOut['ell'][(kpc > 50.0) & (kpc < 120.0)])
avgPA3 = np.nanmedian(ellOut['pa'][(kpc > 50.0) & (kpc < 120.0)])
ba3Arr.append(1.0 - avgEll3)
pa3Arr.append(avgPA3)
# Ellipticity and position angle profiles
smaProf.append(np.asarray(kpc))
ellProf.append(np.asarray(ellOut['ell']))
paProf.append(np.asarray(ellOut['pa']))
sbpProf.append(np.asarray(ellOut['sbp']))
# plot an ellipse for each object
for i in range(len(objects)):
for k, iso in enumerate(ellOut):
if k % 5 == 0:
e = Ellipse(xy=(iso['x0'], iso['y0']),
height=iso['sma'] * 2.0,
width=iso['sma'] * 2.0 * (1.0 - iso['ell']),
angle=iso['pa'])
e.set_facecolor('none')
e.set_edgecolor('r')
e.set_alpha(0.20)
ax.add_artist(e)
ax.set_aspect('equal')
else:
xArr.append(None)
yArr.append(None)
ba1Arr.append(None)
pa1Arr.append(None)
ba2Arr.append(None)
pa2Arr.append(None)
ba3Arr.append(None)
pa3Arr.append(None)
smaProf.append(None)
ellProf.append(None)
paProf.append(None)
sbpProf.append(None)
ax.text(0.5, 0.9, labels[ii],
verticalalignment='center', horizontalalignment='center',
fontsize=40.0, transform=ax.transAxes, weight='bold',
color='r', alpha=0.8)
fig.savefig('mblack2_map_ell_2_example.pdf', dpi=90)
In [55]:
print(ba1Arr[0], pa1Arr[0])
print(ba2Arr[0], pa2Arr[0])
print(ba3Arr[0], pa3Arr[0])
In [56]:
print(xArr[0], yArr[0])
In [68]:
print(scale)
In [67]:
# Prepare the plot
fig = plt.figure(figsize=(14, 14))
fig.subplots_adjust(left=0.0, right=1.0,
bottom=0.0, top=1.0,
wspace=0.00, hspace=0.00)
gs = gridspec.GridSpec(2, 2)
gs.update(wspace=0.0, hspace=0.00)
smaProf1, sbpProf1, cogProf1 = [], [], []
with ProgressBar(len(exampleFile)) as bar:
for item in enumerate(exampleFile):
ii, name = item
img = example[ii]
bar.update()
# Convert the mass density map to just mass maps
imgN = img.byteswap().newbyteorder()
imgC = copy.copy(imgN)
try:
baUse = ba1Arr[0]
paUse = hUtil.normAngle(pa1Arr[0], lower=-90, upper=90, b=True)
ellOut, binOut = galSBP.galSBP(name,
galX=xArr[0], galY=yArr[0],
maxSma=108, iniSma=15.0,
verbose=False, savePng=False, saveOut=False,
useZscale=False, expTime=1.0, pix=scale, zpPhoto=0.0,
galQ=baUse, galPA=paUse,
stage=3, minSma=0.0, ellipStep=0.08,
isophote=isophote, xttools=xttools,
fracBad=0.6, uppClip=6.0, lowClip=6.0, maxTry=9,
nClip=0, intMode='median', updateIntens=False)
except Exception:
ellOut, binOut = None, None
# Show the aperture
ax = plt.subplot(gs[ii])
ax.yaxis.set_major_formatter(NullFormatter())
ax.xaxis.set_major_formatter(NullFormatter())
m, s = np.mean(imgC), np.std(imgC)
im = ax.imshow(imgC, interpolation='nearest', cmap=CUBE,
vmin=m-1.5*s, vmax=m+10*s, origin='lower')
if ellOut is not None:
sma = ellOut['sma']
kpc = sma * scale
smaProf1.append(np.asarray(kpc))
cogProf1.append(np.asarray(ellOut['growth_cor']))
sbpProf1.append(np.asarray(ellOut['sbp']) / -2.5)
# plot an ellipse for each object
for i in range(len(objects)):
for k, iso in enumerate(ellOut):
if k % 5 == 0:
e = Ellipse(xy=(iso['x0'], iso['y0']),
height=iso['sma'] * 2.0,
width=iso['sma'] * 2.0 * (1.0 - iso['ell']),
angle=iso['pa'])
e.set_facecolor('none')
e.set_edgecolor('r')
e.set_alpha(0.20)
ax.add_artist(e)
ax.set_aspect('equal')
else:
smaProf1.append(None)
cogProf1.append(None)
sbpProf1.append(None)
ax.text(0.5, 0.9, labels[ii],
verticalalignment='center', horizontalalignment='center',
fontsize=40.0, transform=ax.transAxes, weight='bold',
color='r', alpha=0.8)
fig.savefig('mblack2_map_ell_3_example.pdf', dpi=90)
In [49]:
hscAvgProf0 = loadPkl("hscAvgProf0.pkl")
hscAvgProf1 = loadPkl("hscAvgProf1.pkl")
hscAvgProf2 = loadPkl("hscAvgProf2.pkl")
rm0_sl, rm0_ml, rm0_aml = hscAvgProf0['all'], hscAvgProf0['med'], hscAvgProf0['avg']
rm1_sl, rm1_ml, rm1_aml = hscAvgProf1['all'], hscAvgProf1['med'], hscAvgProf1['avg']
rm2_sl, rm2_ml, rm2_aml = hscAvgProf2['all'], hscAvgProf2['med'], hscAvgProf2['avg']
In [66]:
vline1, vline2 = 10.0, 100.0
matchR, highlight1, highlight2 = 100.0, False, True
xmin, xmax = 1.02, 3.99
#-------------------------------------------------------------------------------#
ymin, ymax = 4.01, 9.89
norm, integrate, normR1 = False, False, 10.0
dmin, dmax = -0.199, 0.399
showLegend = True
#-------------------------------------------------------------------------------#
label1="$11.4 < \log{M_{\star,\ 100\mathrm{kpc}}} \leq 11.6$"
label2="$11.6 < \log{M_{\star,\ 100\mathrm{kpc}}} \leq 11.8$"
label3="$11.8 < \log{M_{\star,\ 100\mathrm{kpc}}} \leq 12.0$"
showInfo1=True
showInfo2=True
showLegend=True
rPsfKpc=5.5
kpcArr=[2.0, 5.0, 10.0, 20.0, 50.0, 100.0, 150.0]
color1a=BLK(0.5)
color1b=BLK(0.7)
cmap1=BLK
color2a=ORG(0.5)
color2b=ORG(0.7)
cmap2=ORG
xtickFormat='$\mathbf{%4.1f}$'
ytickFormat='$\mathbf{%g}\ $'
ytickFormat2='$\mathbf{%g}$'
# --------------------------------------------------------------------------------------- #
fig = plt.figure(figsize=(14, 14))
fig.subplots_adjust(left=0.13, right=0.995,
bottom=0.13, top=0.94,
wspace=0.00, hspace=0.00)
# --------------------------------------------------------------------------------------- #
ax1 = fig.add_subplot(111)
ax1 = songPlotSetup(ax1, ylabel=50, xlabel=50, border=6.5,
xtickFormat=xtickFormat, ytickFormat=ytickFormat)
# --------------------------------------------------------------------------------------- #
## Mark the two interesting radius
if highlight1:
ax1.axvline(vline1 ** 0.25, linewidth=5.5, c='k', linestyle='--',
zorder=0, alpha=0.5, dashes=(30, 6))
else:
ax1.axvline(vline1 ** 0.25, linewidth=4.0, c='k', linestyle='--',
zorder=0, alpha=0.2)
if highlight2:
ax1.axvline(vline2 ** 0.25, linewidth=5.5, c='k', linestyle='--',
zorder=0, alpha=0.5, dashes=(30, 6))
else:
ax1.axvline(vline2 ** 0.25, linewidth=4.0, c='k', linestyle='--',
zorder=0, alpha=0.2)
# --------------------------------------------------------------------------------------- #
## Individual profiles
for ii, aa in enumerate(rm0_sl):
if (ii % 7) == 0:
ax1.plot(RSMA_COMMON, aa, c=color1a, alpha=0.10, linewidth=1.1)
for ii, bb in enumerate(rm1_sl):
if (ii % 7) == 0:
ax1.plot(RSMA_COMMON, bb, c=color1a, alpha=0.15, linewidth=1.1)
for ii, cc in enumerate(rm2_sl):
if (ii % 7) == 0:
ax1.plot(RSMA_COMMON, cc, c=color1a, alpha=0.15, linewidth=1.1)
# --------------------------------------------------------------------------------------- #
## Median profiles
ax1.plot(RSMA_COMMON, rm0_aml[2], linestyle='--', linewidth=10.0,
c=cmap1(0.9), alpha=0.9, zorder=8, label=label1, dashes=(10,6))
ax1.plot(RSMA_COMMON, rm1_aml[2], linestyle='-', linewidth=10.0,
c=cmap1(0.9), alpha=0.9, zorder=8, label=label2)
ax1.plot(RSMA_COMMON, rm2_aml[2], linestyle='--', linewidth=10.0,
c=cmap1(0.9), alpha=0.8, zorder=8, label=label3, dashes=(24,8))
# ========= #
ax1.plot((smaProf1[0] * 1.0) ** 0.25, sbpProf1[0],
linestyle='-', linewidth=5.0,
c=BLU(0.7), alpha=0.9, zorder=0)
#ax1.plot((smaProf1[1] * 0.7) ** 0.25, sbpProf1[1],
# linestyle='-', linewidth=3.0,
# c=BLU(0.9), alpha=0.4, zorder=0)
ax1.plot((smaProf1[3] * 1.0) ** 0.25, sbpProf1[3],
linestyle='--', linewidth=5.0,
c=GRN(0.9), alpha=0.5, zorder=0, dashes=(40,8))
# ========= #
# --------------------------------------------------------------------------------------- #
## Y Lables
if norm:
ax1.set_ylabel('$\mathrm{Normalized}\ \ \log ({\mu}_{\star}/[M_{\odot}\ \mathrm{kpc}^{-2}])$',
size=41)
else:
ax1.set_ylabel('$\log ({\mu}_{\star}/[M_{\odot}\ \mathrm{kpc}^{-2}])$', size=65)
ax1.set_xlabel('$R^{1/4}\ (\mathrm{kpc})$', size=60)
# --------------------------------------------------------------------------------------- #
## X, Y limits
ax1.set_xlim(xmin, xmax)
ax1.set_ylim(ymin, ymax)
# --------------------------------------------------------------------------------------- #
## Region affected by PSF
# z = 0.4 : 1"=5.4 Kpc
# z = 0.5 : 1"=6.1 Kpc
ySep = (ymax - ymin) / 5.0
ax1.fill_between([0.0, rPsfKpc ** 0.25],
[ymin - ySep, ymin - ySep],
[ymax + ySep, ymax + ySep],
facecolor='k', edgecolor='k', alpha=0.15, zorder=0)
## Label the PSF region
ax1.text(0.05, 0.20, '$\mathrm{PSF}$', rotation='vertical',
verticalalignment='bottom', horizontalalignment='left',
fontsize=50.0, transform=ax1.transAxes, weight='bold',
color='k', alpha=0.4)
# --------------------------------------------------------------------------------------- #
## Legend
if showLegend:
ax1.legend(loc=(0.325, 0.72), shadow=True, fancybox=True,
numpoints=1, fontsize=40, scatterpoints=1,
markerscale=1.2, borderpad=0.4, handletextpad=0.6)
# --------------------------------------------------------------------------------------- #
## Secondary Axis
ax2 = ax1.twiny()
kpcs = np.asarray(kpcArr)
kpcTicks= (kpcs ** 0.25)
ax2.set_xlim(xmin, xmax)
ax2.set_xticks(kpcTicks)
ax2.set_xticklabels(['$\mathbf{%g}$' % kpc for kpc in kpcs], fontsize=50)
for tick in ax2.xaxis.get_major_ticks():
tick.label.set_fontsize(40)
for tick in ax2.yaxis.get_major_ticks():
tick.label.set_fontsize(40)
ax2.text(0.92, 1.0040, '$\mathrm{kpc}$',
verticalalignment='bottom', horizontalalignment='left',
fontsize=50.0, transform=ax2.transAxes)
# --------------------------------------------------------------------------------------- #
ax1.text(0.41, 0.13, '$\mathrm{Average\ Profiles\ of}$',
verticalalignment='bottom', horizontalalignment='center',
fontsize=50.0, transform=ax1.transAxes)
ax1.text(0.41, 0.04, '$\mathrm{Different\ }M_{\star}\ \mathrm{Bins}$',
verticalalignment='bottom', horizontalalignment='center',
fontsize=50.0, transform=ax1.transAxes)
# --------------------------------------------------------------------------------------- #
vline1, vline2 = 10.0, 100.0
matchR, highlight1, highlight2 = 100.0, False, True
xmin, xmax = 1.02, 3.99
#-------------------------------------------------------------------------------#
ymin, ymax = 4.01, 9.89
norm, integrate, normR1 = False, False, 10.0
dmin, dmax = -0.199, 0.399
showLegend = True
showInfo1=True
showInfo2=True
showLegend=True
rPsfKpc=5.5
kpcArr=[2.0, 5.0, 10.0, 20.0, 50.0, 100.0, 150.0]
color1a=BLK(0.5)
color1b=BLK(0.7)
cmap1=BLK
color2a=ORG(0.5)
color2b=ORG(0.7)
cmap2=ORG
xtickFormat='$\mathbf{%4.1f}$'
ytickFormat='$\mathbf{%g}$'
ytickFormat2='$\mathbf{%g}$'
fig.savefig('mblack2_average_mass_profiles_example.pdf', dpi=90)
plt.show()
In [6]:
cenX, cenY = 200.0, 200.0
scale = 1.0
In [17]:
# Prepare the
fig = plt.figure(figsize=(14, 76))
fig.subplots_adjust(left=0.0, right=1.0,
bottom=0.0, top=1.0,
wspace=0.00, hspace=0.00)
gs = gridspec.GridSpec(38, 7)
gs.update(wspace=0.0, hspace=0.00)
xArr, yArr, baArr, paArr, m10Arr, m30Arr, m100Arr, mtotArr = [], [], [], [], [], [], [], []
m150Arr, m200Arr = [], []
aperProf, aperMass = [], []
r50Arr = []
with ProgressBar(len(e)) as bar:
for item in enumerate(cat):
ii, col = item
bar.update()
imgN = img.byteswap().newbyteorder()
# Extract objects
bkg = sep.Background(imgN, bw=20, bh=20, fw=4, fh=4)
imgSub = imgN - bkg
objects = sep.extract(imgSub, 15.0, filter_kernel=kernel3)
# Find the object at the center
if len(objects) == 1:
index = 0
else:
index = np.argmin(np.sqrt((objects['x'] - cenX) ** 2.0 + (objects['y'] - cenY) ** 2.0))
#print(len(objects), index)
# Get the naive ba, theta, xcen, ycen
ba = objects[index]['b'] / objects[index]['a']
theta = objects[index]['theta']
xcen, ycen = objects[index]['x'], objects[index]['y']
baArr.append(ba)
paArr.append(theta * 180.0 / np.pi)
xArr.append(xcen)
yArr.append(ycen)
# Aperture photometry and naive "profiles"
rad = np.linspace(1.2, 4.0, 50)
kpc = (rad ** 4.0)
pix = (kpc / scale)
area = (np.pi * (kpc ** 2.0) * ba)
maper = sep.sum_ellipse(imgN, xcen, ycen, pix, pix * ba, theta, 1.0,
bkgann=None, subpix=11)[0]
mring = maper[1:] - maper[0:-1]
aring = area[1:] - area[0:-1]
rho = np.log10(mring / aring)
aperMass.append(maper)
aperProf.append(rho)
# Mass within different apertures
rad = np.asarray([10.0, 30.0, 100.0, 150.0, 200.0]) / scale
maper = sep.sum_ellipse(imgN, xcen, ycen, rad, rad * ba, theta, 1.0,
bkgann=None, subpix=11)[0]
m10Arr.append(np.log10(maper[0]))
m30Arr.append(np.log10(maper[1]))
m100Arr.append(np.log10(maper[2]))
m150Arr.append(np.log10(maper[3]))
m200Arr.append(np.log10(maper[4]))
mtotArr.append(np.log10(np.sum(imgN)))
# Get r50
r50, flag = sep.flux_radius(imgN, xcen, ycen, 30.0, 0.5,
subpix=11)
r50Arr.append(r50 * scale)
# Show the aperture
ax = plt.subplot(gs[ii])
ax.yaxis.set_major_formatter(NullFormatter())
ax.xaxis.set_major_formatter(NullFormatter())
m, s = np.mean(imgN), np.std(imgN)
im = ax.imshow(imgN, interpolation='nearest', cmap=CUBE,
vmin=m-1.5*s, vmax=m+10*s, origin='lower')
# plot an ellipse for each object
for i in range(len(objects)):
e = Ellipse(xy=(xcen, ycen),
width=40 * 2.0,
height=(40 * 2.0 * ba),
angle=(theta * 180. / np.pi))
e.set_facecolor('none')
e.set_edgecolor('r')
e.set_alpha(0.5)
e.set_linewidth(1.0)
ax.add_artist(e)
ax.set_aspect('equal')
ax.text(0.5, 0.15, "2D: %5.2f" % np.log10(maper[3]),
verticalalignment='center', horizontalalignment='center',
fontsize=20.0, transform=ax.transAxes, weight='bold',
color='r', alpha=0.8)
In [11]:
cat.add_column(Column(xArr, name='xCen_aper'))
cat.add_column(Column(yArr, name='yCen_aper'))
cat.add_column(Column(baArr, name='ba_aper'))
cat.add_column(Column(paArr, name='pa_aper'))
cat.add_column(Column(m10Arr, name='m10_aper'))
cat.add_column(Column(m30Arr, name='m30_aper'))
cat.add_column(Column(m100Arr, name='m100_aper'))
cat.add_column(Column(m150Arr, name='m150_aper'))
cat.add_column(Column(m200Arr, name='m200_aper'))
#cat.remove_column('r50_aper')
cat.add_column(Column(r50Arr, name='r50_aper'))
#cat.remove_column('sbp_aper')
cat.add_column(Column(aperProf, name='sbp_aper'))
cat.add_column(Column(aperMass, name='mass_aper'))
In [9]:
xtickFormat='$\mathbf{%4.1f}$'
ytickFormat='$\mathbf{%4.1f}$'
# --------------------------------------------------------------------------------------- #
fig = plt.figure(figsize=(12, 12))
fig.subplots_adjust(left=0.19, right=0.995,
bottom=0.13, top=0.995,
wspace=0.00, hspace=0.00)
# --------------------------------------------------------------------------------------- #
ax1 = fig.add_subplot(111)
ax1 = songPlotSetup(ax1, ylabel=45, xlabel=45, border=6.5,
xtickFormat=xtickFormat, ytickFormat=ytickFormat)
ax1.scatter(np.log10(cat['col6']), cat['m100_aper'],
s=390, c=ORG(0.8), alpha=0.8)
ax1.plot(np.asarray([11.4, 12.2]), np.asarray([11.4, 12.2]), linestyle='--', linewidth=8, c='k', zorder=0)
ax1.set_xlabel('$\log (M_{\star}/M_{\odot})\ \mathrm{3D}$', size=50)
ax1.set_ylabel('$\log (M_{\star}/M_{\odot})\ (100 \mathrm{kpc},\ \mathrm{2D})$', size=50)
# --------------------------------------------------------------------------------------- #
## X, Y limits
ax1.set_xlim(11.45, 12.09)
ax1.set_ylim(11.45, 12.09)
fig.savefig('horizon_mass_compare_m100_z0.37.pdf', dpi=90)
In [10]:
xtickFormat='$\mathbf{%4.1f}$'
ytickFormat='$\mathbf{%4.1f}$'
# --------------------------------------------------------------------------------------- #
fig = plt.figure(figsize=(12, 12))
fig.subplots_adjust(left=0.19, right=0.995,
bottom=0.13, top=0.995,
wspace=0.00, hspace=0.00)
# --------------------------------------------------------------------------------------- #
ax1 = fig.add_subplot(111)
ax1 = songPlotSetup(ax1, ylabel=45, xlabel=45, border=6.5,
xtickFormat=xtickFormat, ytickFormat=ytickFormat)
ax1.scatter(np.log10(cat['col6']), cat['m150_aper'],
s=390, c=ORG(0.8), alpha=0.8)
ax1.plot(np.asarray([11.4, 12.2]), np.asarray([11.4, 12.2]), linestyle='--', linewidth=8, c='k', zorder=0)
ax1.set_xlabel('$\log (M_{\star}/M_{\odot})\ \mathrm{3D}$', size=50)
ax1.set_ylabel('$\log (M_{\star}/M_{\odot})\ (150 \mathrm{kpc},\ \mathrm{2D})$', size=50)
# --------------------------------------------------------------------------------------- #
## X, Y limits
ax1.set_xlim(11.45, 12.09)
ax1.set_ylim(11.45, 12.09)
fig.savefig('horizon_mass_compare_m150_z0.37.pdf', dpi=90)
In [12]:
xtickFormat='$\mathbf{%4.1f}$'
ytickFormat='$\mathbf{%4.1f}$'
rad = np.linspace(1.2, 4.0, 50)
kpc = (rad ** 4.0)
pix = (kpc / 5.0)
r50Arr = []
# --------------------------------------------------------------------------------------- #
fig = plt.figure(figsize=(12, 12))
fig.subplots_adjust(left=0.19, right=0.995,
bottom=0.13, top=0.940,
wspace=0.00, hspace=0.00)
# --------------------------------------------------------------------------------------- #
ax1 = fig.add_subplot(111)
ax1 = songPlotSetup(ax1, ylabel=45, xlabel=45, border=6.5,
xtickFormat=xtickFormat, ytickFormat=ytickFormat)
for gal in cat:
mProf = gal['mass_aper']
mMax = np.nanmax(mProf[(kpc <= 150.0)])
fracMax = (mProf / mMax)
fracInterp = interp1d(fracMax, kpc)
r50 = fracInterp(0.50)
r50Arr.append(r50)
ax1.plot(kpc ** 0.25, mProf / mMax, linestyle='--', linewidth=2, c=BLU(0.8), alpha=0.7)
ax1.axhline(0.5, linestyle='-', linewidth=5.0, c=BLK(0.6), alpha=0.5, zorder=0)
## Label
ax1.set_xlabel('$R^{1/4}\ (\mathrm{kpc})$', size=46)
ax1.set_ylabel('$\mathrm{Fraction}$', size=45)
# --------------------------------------------------------------------------------------- #
## X, Y limits
ax1.set_xlim(1.05, 3.75)
ax1.set_ylim(0.02, 1.19)
# --------------------------------------------------------------------------------------- #
## Secondary Axis
kpcArr = [2.0, 5.0, 10.0, 20.0, 40.0, 60.0, 100.0]
ax1b = ax1.twiny()
kpcs = np.asarray(kpcArr)
kpcTicks= (kpcs ** 0.25)
ax1b.set_xlim(1.05, 3.75)
ax1b.set_xticks(kpcTicks)
ax1b.set_xticklabels(['$\mathbf{%g}$' % kpc for kpc in kpcs], fontsize=35)
for tick in ax1b.xaxis.get_major_ticks():
tick.label.set_fontsize(30)
for tick in ax1b.yaxis.get_major_ticks():
tick.label.set_fontsize(30)
ax1b.text(0.91, 1.0035, '$\mathrm{kpc}$',
verticalalignment='bottom', horizontalalignment='left',
fontsize=35.0, transform=ax1b.transAxes)
fig.savefig('horizon_cog_fraction_hres_z0.37.pdf', dpi=90)
In [13]:
cat.add_column(Column(np.asarray(r50Arr), name='r50_prof'))
In [13]:
xtickFormat='$\mathbf{%4.1f}$'
ytickFormat='$\mathbf{%4.1f}$'
# --------------------------------------------------------------------------------------- #
fig = plt.figure(figsize=(12, 12))
fig.subplots_adjust(left=0.19, right=0.995,
bottom=0.13, top=0.995,
wspace=0.00, hspace=0.00)
# --------------------------------------------------------------------------------------- #
ax1 = fig.add_subplot(111)
ax1 = songPlotSetup(ax1, ylabel=45, xlabel=45, border=6.5,
xtickFormat=xtickFormat, ytickFormat=ytickFormat)
ax1.scatter(cat['m150_aper'], np.log10(cat['r50_prof']),
s=390, c=ORG(0.8), alpha=0.8)
#ax1.plot(np.asarray([11.4, 12.2]), np.asarray([11.4, 12.2]), linestyle='--', linewidth=8, c='k', zorder=0)
ax1.set_xlabel('$\log (M_{\star}/M_{\odot})\ (150 \mathrm{kpc},\ \mathrm{2D})$', size=50)
ax1.set_ylabel('$R_{50}\ (\mathrm{kpc})$', size=50)
# --------------------------------------------------------------------------------------- #
## X, Y limits
ax1.set_xlim(11.41, 12.49)
ax1.set_ylim(0.81, 1.79)
fig.savefig('horizon_m100_r50_hres_z0.37.pdf', dpi=90)
In [14]:
# Prepare the plot
fig = plt.figure(figsize=(14, 76))
fig.subplots_adjust(left=0.0, right=1.0,
bottom=0.0, top=1.0,
wspace=0.00, hspace=0.00)
gs = gridspec.GridSpec(38, 7)
gs.update(wspace=0.0, hspace=0.00)
m10ArrC, m30ArrC, m100ArrC = [], [], []
m150ArrC, m200ArrC = [], []
aperProfC, aperMassC = [], []
with ProgressBar(len(cat)) as bar:
for item in enumerate(cat):
ii, col = item
galID = col['col1']
bar.update()
#print("## Deal with galaxy : %s" % galID)
# Convert the mass density map to just mass maps
imgTest = os.path.join(imgDir, 'galaxy_{:03d}'.format(galID) + '.fits')
imgN = fits.open(imgTest)[0].data
imgN = imgN.byteswap().newbyteorder()
xcen = col['xCen_aper']
ycen = col['yCen_aper']
# Aperture photometry and naive "profiles"
rad = np.linspace(1.2, 4.0, 50)
kpc = (rad ** 4.0)
pix = (kpc / scale)
area = (np.pi * (kpc ** 2.0) * ba)
maper = sep.sum_ellipse(imgN, xcen, ycen, pix, pix, 0.0, 1.0,
bkgann=None, subpix=11)[0]
mring = maper[1:] - maper[0:-1]
aring = area[1:] - area[0:-1]
rho = np.log10(mring / aring)
aperMassC.append(maper)
aperProfC.append(rho)
# Mass within different apertures
rad = np.asarray([10.0, 30.0, 100.0, 150.0, 200.0]) / scale
maper = sep.sum_ellipse(imgN, xcen, ycen, rad, rad, 0.0, 1.0,
bkgann=None, subpix=11)[0]
m10ArrC.append(np.log10(maper[0]))
m30ArrC.append(np.log10(maper[1]))
m100ArrC.append(np.log10(maper[2]))
m150ArrC.append(np.log10(maper[3]))
m200ArrC.append(np.log10(maper[4]))
# Show the aperture
ax = plt.subplot(gs[ii])
ax.yaxis.set_major_formatter(NullFormatter())
ax.xaxis.set_major_formatter(NullFormatter())
m, s = np.mean(imgN), np.std(imgN)
im = ax.imshow(imgN, interpolation='nearest', cmap=CUBE,
vmin=m-1.0*s, vmax=m+10*s, origin='lower')
# plot an ellipse for each object
for i in range(len(objects)):
e = Ellipse(xy=(xcen, ycen),
width=40 * 2.0,
height=(40 * 2.0),
angle=0.0)
e.set_facecolor('none')
e.set_edgecolor('r')
e.set_alpha(0.5)
e.set_linewidth(1.0)
ax.add_artist(e)
ax.set_aspect('equal')
# Show the 2D and 3D mass within 100 kpc
ax.text(0.5, 0.85, "3D: %5.2f" % np.log10(col['col6']),
verticalalignment='center', horizontalalignment='center',
fontsize=20.0, transform=ax.transAxes, weight='bold',
color='r', alpha=0.8)
ax.text(0.5, 0.15, "2D: %5.2f" % np.log10(maper[3]),
verticalalignment='center', horizontalalignment='center',
fontsize=20.0, transform=ax.transAxes, weight='bold',
color='r', alpha=0.8)
fig.savefig('horizon_map_apertureC_z0.37.pdf', dpi=120)