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


    

    




The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

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


    

    




all_xray_imag_lt_23.fits* map_gas_halo.fits*        map_stellar_halo.fits*
map_darkmatter_halo.fits* map_stellar_central.fits*

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]:





37.857142857142854

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))


    

    




(162, 162) (162, 162) (162, 162) (162, 162)
4

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]:





(0.61672132801398216, 19.187912123986127)

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)


    

    




!!! DANGEROUS : Outer boundary is not safe !!!
!!! DANGEROUS : Outer boundary is not safe !!!
!!! DANGEROUS : Outer boundary is not safe !!!
!!! DANGEROUS : Outer boundary is not safe !!!






    

In [55]:

    

print(ba1Arr[0], pa1Arr[0])
print(ba2Arr[0], pa2Arr[0])
print(ba3Arr[0], pa3Arr[0])


    

    




0.6492024 -66.13204
0.68869925 -64.757375
0.6090984 -67.93168

In [56]:

    

print(xArr[0], yArr[0])


    

    




81.0 82.0

In [68]:

    

print(scale)


    

    




2.65714285714

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)


    

    




!!! DANGEROUS : Outer boundary is not safe !!!
!!! DANGEROUS : Outer boundary is not safe !!!
!!! DANGEROUS : Outer boundary is not safe !!!
!!! DANGEROUS : Outer boundary is not safe !!!






    

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)


    

    




---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-17-41cef06dce52> in <module>()
     12 r50Arr = []
     13 
---> 14 with ProgressBar(len(cat)) as bar:
     15     for item in enumerate(cat):
     16         ii, col = item

NameError: name 'cat' is not defined





    






<matplotlib.figure.Figure at 0x113598650>

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)


    

    




/Users/song/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:38: RuntimeWarning: divide by zero encountered in log10






    

In [15]:

    

cat.add_column(Column(m10ArrC, name='m10c_aper'))
cat.add_column(Column(m30ArrC, name='m30c_aper'))
cat.add_column(Column(m100ArrC, name='m100c_aper'))
cat.add_column(Column(m150ArrC, name='m150c_aper'))
cat.add_column(Column(m200ArrC, name='m200c_aper'))


    

In [16]:

    

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['m100c_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,\ Circ})$', size=46)
# --------------------------------------------------------------------------------------- #
## X, Y limits
ax1.set_xlim(11.45, 12.09)
ax1.set_ylim(11.45, 12.09)

fig.savefig('horizon_m100c_compare_z0.37.pdf', dpi=90)