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
In [3]:
cat = Table.read('catalog_highres_z0.37.fits', format='fits')
mstar = 'stellar_mass_maps'
mtot = 'total_mass_maps'
print(cat.colnames)
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]:
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)
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)
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)