In [2]:

    

%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.pyplot as plt
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 [3]:

    

cat = Table.read('catalog_highres_z0.37.fits', format='fits')

mstar = 'stellar_mass_maps'
mtot = 'total_mass_maps'

print(cat.colnames)


    

    




['id', 'id_subhalo', 'central', 'logM500', 'r500', 'logM200', 'r200', 'logM30_3D', 'logM100_3D']

In [4]:

    

isophote = '/iraf/iraf/extern/stsdas/bin.macosx/x_isophote.e'
xttools = '/iraf/iraf/extern/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) / 7


    

    
Out[5]:





19.0

In [6]:

    

# Prepare the plot
fig = plt.figure(figsize=(14, 38))
fig.subplots_adjust(left=0.0, right=1.0, 
                    bottom=0.0, top=1.0,
                    wspace=0.00, hspace=0.00)
gs = gridspec.GridSpec(19, 7)
gs.update(wspace=0.0, hspace=0.00)

xArr, yArr, baArr, paArr, m10Arr, m30Arr, m100Arr, mtotArr = [], [], [], [], [], [], [], []
aperProf, aperMass = [], []
r50Arr = []

with ProgressBar(len(cat)) as bar:
    for item in enumerate(cat):
        ii, col = item
        galID = str(col['id'])
        bar.update()
        #print("## Deal with galaxy : %s" % galID)

        # Convert the mass density map to just mass maps
        imgTest = os.path.join(mstar, 'galaxy_' + galID + '_stars_no_sats.fits')
        imgSave = os.path.join(mstar, 'galaxy_' + galID + '_stars_no_sats_ms.fits')

        imgN = fits.open(imgTest)[0].data / 1.0E6 * 25.0 
        hdu = fits.PrimaryHDU(imgN)
        hdulist = fits.HDUList([hdu])
        hdulist.writeto(imgSave, clobber=True)

        # 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'] - 50.0) ** 2.0 + (objects['y'] - 50.0) ** 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 / 5.0)
        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, 250.0]) / 5.0
        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]))
        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 * 5.0)

        # 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=20 * 2.0,
                        height=(20 * 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.85, "3D: %5.2f" % col['logM100_3D'], 
                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[2]), 
                verticalalignment='center', horizontalalignment='center',
                fontsize=20.0, transform=ax.transAxes, weight='bold', 
                color='r', alpha=0.8)
    
fig.savefig('bhamas_map_aperture_hres_z0.37.pdf', dpi=140)


    

    




/usr/local/lib/python2.7/site-packages/ipykernel/__main__.py:60: RuntimeWarning: divide by zero encountered in log10






    

In [7]:

    

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.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 [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(cat['logM100_3D'], 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})\ (100 \mathrm{kpc},\ \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('m100_compare_hres_z0.37.pdf', dpi=90)


    

In [8]:

    

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)

#ax1.plot(np.asarray([11.4, 12.2]), np.asarray([11.4, 12.2]), linestyle='--', linewidth=8, c='k', 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('cog_fraction_hres_z0.37.pdf', dpi=90)


    

In [9]:

    

cat.add_column(Column(np.asarray(r50Arr), name='r50_prof'))


    

In [14]:

    

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['m100_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})\ (100 \mathrm{kpc},\ \mathrm{2D})$', size=50)
ax1.set_ylabel('$R_{50}\ (\mathrm{kpc})$', size=50)
# --------------------------------------------------------------------------------------- #
## X, Y limits
ax1.set_xlim(11.49, 12.19)
ax1.set_ylim(0.61, 1.79)

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


    

In [10]:

    

# Prepare the plot
fig = plt.figure(figsize=(14, 38))
fig.subplots_adjust(left=0.0, right=1.0, 
                    bottom=0.0, top=1.0,
                    wspace=0.00, hspace=0.00)
gs = gridspec.GridSpec(19, 7)
gs.update(wspace=0.0, hspace=0.00)

m10ArrC, m30ArrC, m100ArrC = [], [], []
aperProfC, aperMassC = [], []

with ProgressBar(len(cat)) as bar:
    for item in enumerate(cat):
        ii, col = item
        galID = str(col['id'])
        bar.update()
        #print("## Deal with galaxy : %s" % galID)

        # Convert the mass density map to just mass maps
        imgTest = os.path.join(mstar, 'galaxy_' + galID + '_stars_no_sats.fits')
        imgN = fits.open(imgTest)[0].data / 1.0E6 * 25.0 
        
        # 
        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 / 5.0)
        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]) / 5.0
        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]))

        # 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=20 * 2.0,
                        height=(20 * 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" % col['logM100_3D'], 
                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[2]), 
                verticalalignment='center', horizontalalignment='center',
                fontsize=20.0, transform=ax.transAxes, weight='bold', 
                color='r', alpha=0.8)
    
fig.savefig('bhamas_map_apertureC_hres_z0.37.pdf', dpi=140)


    

    




/usr/local/lib/python2.7/site-packages/ipykernel/__main__.py:36: RuntimeWarning: divide by zero encountered in log10






    

In [11]:

    

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


    

In [18]:

    

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['logM100_3D'], 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})\ (100 \mathrm{kpc},\ \mathrm{3D,\ Circ})$', size=46)
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('m100c_compare_hres_z0.37.pdf', dpi=90)


    

In [12]:

    

# Prepare the plot
fig = plt.figure(figsize=(14, 38))
fig.subplots_adjust(left=0.0, right=1.0, 
                    bottom=0.0, top=1.0,
                    wspace=0.00, hspace=0.00)
gs = gridspec.GridSpec(19, 7)
gs.update(wspace=0.0, hspace=0.00)

xArr, yArr, ba1Arr, pa1Arr, ba2Arr, pa2Arr, ba3Arr, pa3Arr = [], [], [], [], [], [], [], []
smaProf, ellProf, paProf, sbpProf = [], [], [], []

with ProgressBar(len(cat)) as bar:
    for item in enumerate(cat):
        ii, col = item
        galID = str(col['id'])
        bar.update()
        #print("## Deal with galaxy : %s" % galID)

        # Convert the mass density map to just mass maps
        imgTest = os.path.join(mstar, 'galaxy_' + galID + '_stars_no_sats.fits')
        imgSave = os.path.join(mstar, 'galaxy_' + galID + '_stars_no_sats_ms_conv.fits')

        imgN = fits.open(imgTest)[0].data / 1.0E6 * 25.0 
        imgC = ndimage.convolve(imgN, kernel2, mode='constant', cval=0.0) 
        
        hdu = fits.PrimaryHDU(imgC)
        hdulist = fits.HDUList([hdu])
        hdulist.writeto(imgSave, clobber=True)

        try:
            ellOut, binOut = galSBP.galSBP(imgSave, galX=50.0, galY=50.0, 
                                           maxSma=32, iniSma=5.0, 
                                           verbose=False, savePng=False, saveOut=False,
                                           useZscale=False, expTime=1.0, pix=5.0, zpPhoto=0.0,
                                           galQ=col['ba_aper'], galPA=col['pa_aper'], 
                                           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='mean', 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.0*s, vmax=m+10*s, origin='lower')
        
        if ellOut is not None:
            # Estimate the average ellipticity and position angle 
            sma = ellOut['sma']
            kpc = sma * 5.0
            xArr.append(ellOut[1]['x0'])
            yArr.append(ellOut[1]['y0'])

            avgEll1 = np.nanmedian(ellOut['ell'][(kpc > 12.0) & (kpc < 40.0)])
            avgPA1 = np.nanmedian(ellOut['pa'][(kpc > 12.0) & (kpc < 40.0)])
            ba1Arr.append(1.0 - avgEll1)
            pa1Arr.append(avgPA1)
            avgEll2 = np.nanmedian(ellOut['ell'][(kpc > 40.0) & (kpc < 80.0)])
            avgPA2 = np.nanmedian(ellOut['pa'][(kpc > 40.0) & (kpc < 80.0)])
            ba2Arr.append(1.0 - avgEll2)
            pa2Arr.append(avgPA2)
            avgEll3 = np.nanmedian(ellOut['ell'][(kpc > 80.0) & (kpc < 120.0)])
            avgPA3 = np.nanmedian(ellOut['pa'][(kpc > 80.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.15)
                    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)            
            
fig.savefig('bhamas_map_ellip_step2_hres_z0.37.pdf', dpi=160)


    

    




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
###  ELLIPSE RUN FAILED IN ATTEMPT:  0
###  Error Information :  XXX Can not find the outBin: stellar_mass_maps/galaxy_2_stars_no_sats_ms_conv_ellip_2.bin!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
----------------------------------------------------------------------------------------------------
###  Maxsma   32.0 -->   28.8
----------------------------------------------------------------------------------------------------
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
###  ELLIPSE RUN FAILED IN ATTEMPT:  1
###  Error Information :  XXX Can not find the outBin: stellar_mass_maps/galaxy_2_stars_no_sats_ms_conv_ellip_2.bin!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
----------------------------------------------------------------------------------------------------
###  Maxsma   28.8 -->   25.9
----------------------------------------------------------------------------------------------------
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
###  ELLIPSE RUN FAILED IN ATTEMPT:  2
###  Error Information :  XXX Can not find the outBin: stellar_mass_maps/galaxy_2_stars_no_sats_ms_conv_ellip_2.bin!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
----------------------------------------------------------------------------------------------------
###  Maxsma   25.9 -->   23.3
----------------------------------------------------------------------------------------------------
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
###  ELLIPSE RUN FAILED IN ATTEMPT:  3
###  Error Information :  XXX Can not find the outBin: stellar_mass_maps/galaxy_2_stars_no_sats_ms_conv_ellip_2.bin!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
----------------------------------------------------------------------------------------------------
###  Maxsma   23.3 -->   21.0
----------------------------------------------------------------------------------------------------
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
###  ELLIPSE RUN FAILED IN ATTEMPT:  4
###  Error Information :  XXX Can not find the outBin: stellar_mass_maps/galaxy_2_stars_no_sats_ms_conv_ellip_2.bin!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
----------------------------------------------------------------------------------------------------
###  Maxsma   21.0 -->   18.9
###  Step     0.08 -->   0.09
----------------------------------------------------------------------------------------------------
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
###  ELLIPSE RUN FAILED IN ATTEMPT:  5
###  Error Information :  XXX Can not find the outBin: stellar_mass_maps/galaxy_2_stars_no_sats_ms_conv_ellip_2.bin!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
----------------------------------------------------------------------------------------------------
###  Maxsma   18.9 -->   17.0
###  Step     0.09 -->   0.10
###  Flag     0.60 -->   0.63
----------------------------------------------------------------------------------------------------
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
###  ELLIPSE RUN FAILED IN ATTEMPT:  6
###  Error Information :  XXX Can not find the outBin: stellar_mass_maps/galaxy_2_stars_no_sats_ms_conv_ellip_2.bin!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
----------------------------------------------------------------------------------------------------
###  Maxsma   17.0 -->   15.3
###  Step     0.10 -->   0.10
###  Flag     0.63 -->   0.66
----------------------------------------------------------------------------------------------------
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
###  ELLIPSE RUN FAILED IN ATTEMPT:  7
###  Error Information :  XXX Can not find the outBin: stellar_mass_maps/galaxy_2_stars_no_sats_ms_conv_ellip_2.bin!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
----------------------------------------------------------------------------------------------------
###  Maxsma   15.3 -->   13.8
###  Step     0.10 -->   0.11
###  Flag     0.66 -->   0.69
----------------------------------------------------------------------------------------------------
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
###  ELLIPSE RUN FAILED IN ATTEMPT:  8
###  Error Information :  XXX Can not find the outBin: stellar_mass_maps/galaxy_2_stars_no_sats_ms_conv_ellip_2.bin!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
----------------------------------------------------------------------------------------------------
###  Maxsma   13.8 -->   12.4
###  Step     0.11 -->   0.12
###  Flag     0.69 -->   0.72
----------------------------------------------------------------------------------------------------
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
###  ELLIPSE RUN FAILED AFTER   9 ATTEMPTS!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!






    




/Users/songhuang/Dropbox/work/project/hs_hsc/py/galsbp/galSBP.py:577: RuntimeWarning: divide by zero encountered in log10
  sbpOri = zp - 2.5 * np.log10(intensOri / (pixArea * exptime))
/Users/songhuang/Dropbox/work/project/hs_hsc/py/galsbp/galSBP.py:578: RuntimeWarning: divide by zero encountered in log10
  sbpSub = zp - 2.5 * np.log10(intensSub / (pixArea * exptime))
/Users/songhuang/Dropbox/work/project/hs_hsc/py/galsbp/galSBP.py:588: RuntimeWarning: divide by zero encountered in log10
  (pixArea * exptime))
/Users/songhuang/Dropbox/work/project/hs_hsc/py/galsbp/galSBP.py:589: RuntimeWarning: invalid value encountered in subtract
  sbp_err = (sbpSub - sbp_low)
/Users/songhuang/Dropbox/work/project/hs_hsc/py/galsbp/galSBP.py:1617: RuntimeWarning: divide by zero encountered in log10
  expTime))






    

In [13]:

    

try:
    cat.remove_column('xCen_ellip')
    cat.remove_column('yCen_ellip')
    cat.remove_column('ba1_ellip')
    cat.remove_column('pa1_ellip')
    cat.remove_column('ba2_ellip')
    cat.remove_column('pa2_ellip')
    cat.remove_column('ba3_ellip')
    cat.remove_column('pa3_ellip')
    cat.remove_column('sma2')
    cat.remove_column('ell_prof')
    cat.remove_column('pa_prof')
    cat.remove_column('sbp2_prof')
except Exception:
    pass

cat.add_column(Column(xArr, name='xCen_ellip'))
cat.add_column(Column(yArr, name='yCen_ellip'))
cat.add_column(Column(ba1Arr, name='ba1_ellip'))
cat.add_column(Column(pa1Arr, name='pa1_ellip'))
cat.add_column(Column(ba2Arr, name='ba2_ellip'))
cat.add_column(Column(pa2Arr, name='pa2_ellip'))
cat.add_column(Column(ba3Arr, name='ba3_ellip'))
cat.add_column(Column(pa3Arr, name='pa3_ellip'))

cat.add_column(Column(smaProf, name='sma2'))
cat.add_column(Column(ellProf, name='ell_prof'))
cat.add_column(Column(paProf, name='pa_prof'))
cat.add_column(Column(sbpProf, name='sbp2_prof'))


    

In [24]:

    

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['ba_aper'], cat['ba1_ellip'], s=390, c=ORG(0.8), alpha=0.6, marker='s')
ax1.scatter(cat['ba_aper'], cat['ba2_ellip'], s=390, c=GRN(0.8), alpha=0.6, marker='o')
ax1.scatter(cat['ba_aper'], cat['ba3_ellip'], s=390, c=BLU(0.8), alpha=0.6, marker='h')

ax1.plot(np.asarray([0.2, 1.0]), np.asarray([0.2, 1.0]), linestyle='--', linewidth=8, c='k', zorder=0)

ax1.set_xlabel('$(b/a)\ \mathrm{SEP}$', size=50)
ax1.set_ylabel('$(b/a)\ \mathrm{ELLIPSE}$', size=50)
# --------------------------------------------------------------------------------------- #
## X, Y limits
ax1.set_xlim(0.41, 0.995)
ax1.set_ylim(0.41, 0.995)

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


    

In [14]:

    

# Prepare the plot
fig = plt.figure(figsize=(14, 38))
fig.subplots_adjust(left=0.0, right=1.0, 
                    bottom=0.0, top=1.0,
                    wspace=0.00, hspace=0.00)
gs = gridspec.GridSpec(19, 7)
gs.update(wspace=0.0, hspace=0.00)

smaProf1, sbpProf1, cogProf1 = [], [], []

with ProgressBar(len(cat)) as bar:
    for item in enumerate(cat):
        ii, col = item
        galID = str(col['id'])
        bar.update()
        #print("## Deal with galaxy : %s" % galID)

        # Convert the mass density map to just mass maps
        imgTest = os.path.join(mstar, 'galaxy_' + galID + '_stars_no_sats.fits')
        imgSave = os.path.join(mstar, 'galaxy_' + galID + '_stars_no_sats_ms.fits')
        imgN = fits.open(imgTest)[0].data / 1.0E6 * 25.0 
        
        try:
            baUse = col['ba1_ellip']
            paUse = hUtil.normAngle(col['pa1_ellip'], lower=-90, upper=90, b=True)
            ellOut, binOut = galSBP.galSBP(imgSave, galX=col['xCen_ellip'], galY=col['yCen_ellip'], 
                                           maxSma=50, iniSma=4.0, 
                                           verbose=False, savePng=False, saveOut=False,
                                           useZscale=False, expTime=1.0, pix=5.0, zpPhoto=0.0,
                                           galQ=baUse, galPA=paUse, 
                                           stage=3, minSma=0.0, ellipStep=0.10,
                                           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(imgN), np.std(imgN)
        im = ax.imshow(imgN, interpolation='nearest', cmap=CUBE,
                       vmin=m-1.0*s, vmax=m+10*s, origin='lower')
        
        if ellOut is not None:
            # Estimate the average ellipticity and position angle 
            sma = ellOut['sma']
            kpc = sma * 5.0

            # Ellipticity and position angle profiles
            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.15)
                    ax.add_artist(e)
            ax.set_aspect('equal')
        else: 
            # Ellipticity and position angle profiles
            smaProf1.append(None)
            cogProf1.append(None)
            sbpProf1.append(None)
        
fig.savefig('bhamas_map_ellip_step3_1_hres_z0.37.pdf', dpi=150)


    

In [15]:

    

#cat.remove_column('sma3a')
#cat.remove_column('sbp3a_prof')
#cat.remove_column('cog3a_prof')
cat.add_column(Column(smaProf1, name='sma3a'))
cat.add_column(Column(sbpProf1, name='sbp3a_prof'))
cat.add_column(Column(cogProf1, name='cog3a_prof'))


    

In [36]:

    

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)

m100e = []
for gal in cat: 
    try:
        intrp = interp1d(gal['sma3a'], gal['cog3a_prof'])
        m100e.append(np.log10(intrp(100)) + 0.05)
    except Exception:
        m100e.append(None)
m100e = np.asarray(m100e)
    
ax1.scatter(cat['m100_aper'], m100e, 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})\ (100 \mathrm{kpc},\ \mathrm{2D})$', size=50)
ax1.set_ylabel('$\log (M_{\star}/M_{\odot})\ (100 \mathrm{kpc},\ \mathrm{1D})$', size=50)
# --------------------------------------------------------------------------------------- #
## X, Y limits
ax1.set_xlim(11.45, 12.09)
ax1.set_ylim(11.45, 12.09)

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


    

In [16]:

    

# Prepare the plot
fig = plt.figure(figsize=(14, 38))
fig.subplots_adjust(left=0.0, right=1.0, 
                    bottom=0.0, top=1.0,
                    wspace=0.00, hspace=0.00)
gs = gridspec.GridSpec(19, 7)
gs.update(wspace=0.0, hspace=0.00)

smaProf2, sbpProf2, cogProf2 = [], [], []

with ProgressBar(len(cat)) as bar:
    for item in enumerate(cat):
        ii, col = item
        galID = str(col['id'])
        bar.update()
        #print("## Deal with galaxy : %s" % galID)

        # Convert the mass density map to just mass maps
        imgTest = os.path.join(mstar, 'galaxy_' + galID + '_stars_no_sats.fits')
        imgSave = os.path.join(mstar, 'galaxy_' + galID + '_stars_no_sats_ms.fits')
        imgN = fits.open(imgTest)[0].data / 1.0E6 * 25.0 
        
        try:
            baUse = col['ba2_ellip']
            paUse = hUtil.normAngle(col['pa2_ellip'], lower=-90, upper=90, b=True)
            ellOut, binOut = galSBP.galSBP(imgSave, galX=col['xCen_ellip'], galY=col['yCen_ellip'], 
                                           maxSma=50, iniSma=4.0, 
                                           verbose=False, savePng=False, saveOut=False,
                                           useZscale=False, expTime=1.0, pix=5.0, zpPhoto=0.0,
                                           galQ=baUse, galPA=paUse, 
                                           stage=3, minSma=0.0, ellipStep=0.10,
                                           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(imgN), np.std(imgN)
        im = ax.imshow(imgN, interpolation='nearest', cmap=CUBE,
                       vmin=m-1.0*s, vmax=m+10*s, origin='lower')
        
        if ellOut is not None:
            # Estimate the average ellipticity and position angle 
            sma = ellOut['sma']
            kpc = sma * 5.0

            # Ellipticity and position angle profiles
            smaProf2.append(np.asarray(kpc))
            cogProf2.append(np.asarray(ellOut['growth_cor']))
            sbpProf2.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.15)
                    ax.add_artist(e)
            ax.set_aspect('equal')
        else:
            # Ellipticity and position angle profiles
            smaProf2.append(None)
            cogProf2.append(None)
            sbpProf2.append(None)
        
fig.savefig('bhamas_map_ellip_step3_2_hres_z0.37.pdf', dpi=150)


    

In [17]:

    

#cat.remove_column('sma3b')
#cat.remove_column('sbp3b_prof')
#cat.remove_column('cog3b_prof')
cat.add_column(Column(smaProf2, name='sma3b'))
cat.add_column(Column(sbpProf2, name='sbp3b_prof'))
cat.add_column(Column(cogProf2, name='cog3b_prof'))


    

In [40]:

    

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)

m100e = []
for gal in cat: 
    try:
        intrp = interp1d(gal['sma3b'], gal['cog3b_prof'])
        m100e.append(np.log10(intrp(100)) + 0.05)
    except Exception:
        m100e.append(None)
m100e = np.asarray(m100e)
    
ax1.scatter(cat['m100_aper'], m100e, 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})\ (100 \mathrm{kpc},\ \mathrm{2D})$', size=50)
ax1.set_ylabel('$\log (M_{\star}/M_{\odot})\ (100 \mathrm{kpc},\ \mathrm{1D})$', size=50)
# --------------------------------------------------------------------------------------- #
## X, Y limits
ax1.set_xlim(11.45, 12.09)
ax1.set_ylim(11.45, 12.09)

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


    

In [42]:

    

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)

m30e = []
for gal in cat: 
    try:
        intrp = interp1d(gal['sma3b'], gal['cog3b_prof'])
        m30e.append(np.log10(intrp(30)) + 0.05)
    except Exception:
        m30e.append(None)
m30e = np.asarray(m30e)
    
ax1.scatter(cat['m30_aper'], m30e, s=390, c=ORG(0.8), alpha=0.8)

ax1.plot(np.asarray([11.0, 12.2]), np.asarray([11.0, 12.2]), linestyle='--', linewidth=8, c='k', zorder=0)

ax1.set_xlabel('$\log (M_{\star}/M_{\odot})\ (30 \mathrm{kpc},\ \mathrm{2D})$', size=50)
ax1.set_ylabel('$\log (M_{\star}/M_{\odot})\ (30 \mathrm{kpc},\ \mathrm{1D})$', size=50)
# --------------------------------------------------------------------------------------- #
## X, Y limits
ax1.set_xlim(11.05, 12.09)
ax1.set_ylim(11.05, 12.09)

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


    

In [43]:

    

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)

m10e = []
for gal in cat: 
    try:
        intrp = interp1d(gal['sma3b'], gal['cog3b_prof'])
        m10e.append(np.log10(intrp(10)) + 0.05)
    except Exception:
        m10e.append(None)
m10e = np.asarray(m10e)
    
ax1.scatter(cat['m10_aper'], m10e, s=390, c=ORG(0.8), alpha=0.8)

ax1.plot(np.asarray([10.4, 12.2]), np.asarray([10.4, 12.2]), linestyle='--', linewidth=8, c='k', zorder=0)

ax1.set_xlabel('$\log (M_{\star}/M_{\odot})\ (10 \mathrm{kpc},\ \mathrm{2D})$', size=50)
ax1.set_ylabel('$\log (M_{\star}/M_{\odot})\ (10 \mathrm{kpc},\ \mathrm{1D})$', size=50)
# --------------------------------------------------------------------------------------- #
## X, Y limits
ax1.set_xlim(10.45, 11.69)
ax1.set_ylim(10.45, 11.69)

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


    

In [62]:

    

# Prepare the plot
fig = plt.figure(figsize=(14, 50))
fig.subplots_adjust(left=0.0, right=1.0, 
                    bottom=0.0, top=1.0,
                    wspace=0.00, hspace=0.00)
gs = gridspec.GridSpec(25, 7)
gs.update(wspace=0.0, hspace=0.00)

smaDist, sbpDist = [], []

with ProgressBar(len(cat)) as bar:
    for item in enumerate(cat):
        ii, col = item
        galID = str(col['id'])
        bar.update()
        #print("## Deal with galaxy : %s" % galID)

        # Convert the mass density map to just mass maps
        imgTest = os.path.join(mstar, 'galaxy_' + galID + '_stars_no_sats.fits')
        imgS = fits.open(imgTest)[0].data / 1.0E6
        
        try:
            baUse = col['ba2_ellip']
            paUse = hUtil.normAngle(col['pa2_ellip'], lower=-90, upper=90, b=True)
            xUse = col['xCen_ellip']
            yUse = col['yCen_ellip']
        except Exception:
            baUse = 1.0
            paUse = 0.0
            xUse = col['xCen_aper']
            yUse = col['yCen_aper']
        
        xShape, yShape = imgS.shape
        
        rArr = np.asarray([(hUtil.ellipDist(xx+0.5, yy+0.5, xUse, yUse, pa=paUse, q=baUse) * 5.0) for xx, yy in itertools.product(np.arange(xShape), np.arange(yShape))])
        mArr = np.asarray([imgS[xx, yy] for xx, yy in itertools.product(np.arange(xShape), np.arange(yShape))])
        
        mArrU = np.log10(mArr[mArr >= 9900.0])
        rArrU = rArr[mArr >= 9900.0] ** 0.25
        
        # Show the aperture
        ax = plt.subplot(gs[ii])
        ax.yaxis.set_major_formatter(NullFormatter())
        ax.xaxis.set_major_formatter(NullFormatter())

        # Ellipticity and position angle profiles
        smaDist.append(rArrU)
        sbpDist.append(mArrU)

        # plot an ellipse for each object
        ax.scatter(rArrU, mArrU, marker='+', s=10, c=BLU(0.5), alpha=0.4, rasterized=True)
        #ax.set_aspect('equal')
        
        try:
            ax.plot(RSMA_COMMON, rm0_aml[2], linestyle='--', linewidth=1.0, 
                    c=cmap1(0.9), alpha=0.9, zorder=8, label=label1, 
                    dashes=(10,6), rasterized=True)
            ax.plot(RSMA_COMMON, rm1_aml[2], linestyle='-', linewidth=1.0, 
                    c=cmap1(0.9), alpha=0.9, zorder=8, label=label2,
                    rasterized=True)
            ax.plot(RSMA_COMMON, rm2_aml[2], linestyle='--', linewidth=1.0, 
                    c=cmap1(0.9), alpha=0.8, zorder=8, label=label3, 
                    dashes=(24,8), rasterized=True)
        except Exception:
            pass
        
        ax.set_xlim(1.0, 3.5)
        ax.set_ylim(4.0, 11.0)
        
fig.savefig('bhamas_map_dist_hres_z0.37.pdf', dpi=110)