In [36]:
%matplotlib inline
%load_ext autoreload
%autoreload 2
from __future__ import (division,
print_function)
## ---------------------------------------------------------------------------------------------##
## Imports
## ---------------------------------------------------------------------------------------------##
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 scipy.ndimage as ndimage
import cPickle as pickle
## ---------------------------------------------------------------------------------------------##
# Astropy
from astropy.io import fits
from astropy import wcs
from astropy import units as u
from astropy.table import Table, Column, vstack
from astropy.stats import sigma_clip
from astropy.nddata import Cutout2D
from astropy.utils.console import ProgressBar
from astropy import coordinates as coords
from astropy.stats import sigma_clip
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
from matplotlib.patches import Ellipse
from matplotlib.ticker import NullFormatter, MaxNLocator, FormatStrFormatter
from matplotlib.collections import PatchCollection
tickFormat = FormatStrFormatter('$\mathbf{%g}$')
mpl.style.use('classic')
## ---------------------------------------------------------------------------------------------##
# Cosmology
import cosmology
c=cosmology.Cosmo(H0=70.0, omega_m=0.3, omega_l=0.7, flat=1)
## ---------------------------------------------------------------------------------------------##
# 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
## ---------------------------------------------------------------------------------------------##
# Define the region of interests:
from shapely.geometry import Polygon, Point
from descartes import PolygonPatch
## ---------------------------------------------------------------------------------------------##
# Emcee
import emcee
import corner
## ---------------------------------------------------------------------------------------------##
# Personal
import hscUtils as hUtil
import coaddCutoutGalfitSimple as gSimple
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
# Extinction correction factor for HSC
## A\_lambda = Coeff * E(B-V)
A_HSC_G = 3.233
A_HSC_R = 2.291
A_HSC_I = 1.635
A_HSC_Z = 1.261
A_HSC_Y = 1.076
A_G = 0.528
A_I = 0.260
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
RECSCAT = [0.16, 0.11, 0.59, 0.88]
RECHIST = [0.75, 0.11, 0.24, 0.88]
SBP1 = [0.13, 0.12, 0.865, 0.30]
SBP2 = [0.13, 0.42, 0.865, 0.54]
EC1 = [0.135, 0.066, 0.862, 0.30]
EC2 = [0.135, 0.366, 0.862, 0.30]
EC3 = [0.135, 0.666, 0.862, 0.30]
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]
## ---------------------------------------------------------------------------------------------##
# 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
## ---------------------------------------------------------------------------------------------##
In [4]:
def getExtinction(ra, dec, a_lambda=None):
"""
Estimate the Galactic extinction for HSC filters.
Parameters:
ra, dec : The input coordinates can be arrays
"""
# First try mwdust from Jo Bovy
try:
import mwdust
sfd = mwdust.SFD(sf10=True)
from astropy.coordinates import SkyCoord
coords = SkyCoord(ra, dec, frame='icrs', unit='deg')
galactic = coords.galactic
l, b = galactic.l, galactic.b
ebv = sfd(l, b, 0)
except ImportError:
try:
# Then try sncosmo
from sncosmo import SFD98Map
dustDir = os.environ.get('DUST_DIR')
if (not os.path.isfile(os.path.join(dustDir,
'SFD_dust_4096_ngp.fits'))) or (
not os.path.isfile(os.path.join(dustDir,
'SFD_dust_4096_sgp.fits'))):
print('# DUST_DIR : %s' % dustDir)
raise Exception("# Can not find the SFD dust map!")
else:
sfd = SFD98Map(dustDir)
ebv = sfd.get_ebv((ra, dec))
except ImportError:
raise Exception("# Both mwdust and sncosmo are not available")
if a_lambda is not None:
return (ebv * a_lambda)
else:
return ebv
def zscale(img, contrast=0.25, samples=500):
"""
Imaging scaling function.
Parameters:
contrast :
samples :
"""
ravel = img.ravel()
if len(ravel) > samples:
imsort = np.sort(np.random.choice(ravel, size=samples))
else:
imsort = np.sort(ravel)
n = len(imsort)
idx = np.arange(n)
med = imsort[n/2]
w = 0.25
i_lo, i_hi = int((0.5 - w) * n), int((0.5 + w) * n)
p = np.polyfit(idx[i_lo:i_hi], imsort[i_lo:i_hi], 1)
slope, intercept = p
z1 = med - (slope/contrast)*(n/2-n*w)
z2 = med + (slope/contrast)*(n/2-n*w)
return z1, z2
def toColorArr(data, bottom=None, top=None):
"""
Convert a data array to "color array" (between 0 and 1).
Parameters:
bottom, top :
"""
if top is not None:
data[data >= top] = top
if bottom is not None:
data[data <= bottom] = bottom
return ((data - np.nanmin(data)) /
(np.nanmax(data) - np.nanmin(data))) * 255.0
In [5]:
# redshift of Perseus cluster
z_per = 0.017900
# Pixel scale
pix = 0.168
asecScale = 0.366
scale = (asecScale * pix) * u.kpc
# R200 x 0.1
rCluster = (180.0 / asecScale) / 3600.0
# Number of pixels for 1 Kpc at the distance of Perseus
kpcPix = (1.0 / scale)
# Boundary of the useful data
raLow, raUpp = 48.9, 51.0
decLow, decUpp = 40.7, 42.3
# Distance module (error is about 0.15 mag)
distMod = 34.36
cosDim = 0.07332
print("At z=%6.4f, 1 HSC pixel (%5.3f arcsec) equals to %5.3f Kpc" % (z_per, pix, scale.value))
print(" 1 Kpc == %5.3f pixels" % kpcPix.value)
print(" 0.1 x R200 ~ 180 Kpc ~ %6.3f" % rCluster)
print(
"""
# NGC 1272
Cosmology-Corrected Quantities [Ho = 70.50 km/sec/Mpc, Ωmatter = 0.27, Ωvacuum = 0.73]
[Redshift 0.012577 as corrected to the Reference Frame defined by the (Virgo + GA + Shapley)]
Luminosity Distance : 54 Mpc (m-M) = 33.66 mag
Angular-Size Distance : 52.7 Mpc (m-M) = 33.61 mag
Co-Moving Radial Distance : 53.3 Mpc (m-M) = 33.64 mag
Co-Moving Tangential Dist. : 53.3 Mpc (m-M) = 33.64 mag
Co-Moving Volume : 0.000636 Gpc^3
Light Travel-Time : 0.173 Gyr
Age at Redshift 0.012577 : 13.598 Gyr
Age of Universe : 13.770 Gyr
Scale (Cosmology Corrected): 255 pc/arcsec = 0.255 kpc/arcsec = 15.32 kpc/arcmin = 0.92 Mpc/degree
Surface Brightness Dimming : Flux Density per Unit Area = 0.95124; Magnitude per Unit Area = 0.05428 mag""")
print(
"""
# NGC 1275
Luminosity Distance : 74.4 Mpc (m-M) = 34.36 mag
Angular-Size Distance : 71.9 Mpc (m-M) = 34.28 mag
Co-Moving Radial Distance : 73.2 Mpc (m-M) = 34.32 mag
Co-Moving Tangential Dist. : 73.2 Mpc (m-M) = 34.32 mag
Co-Moving Volume : 0.00164 Gpc^3
Light Travel-Time : 0.237 Gyr
Age at Redshift 0.017264 : 13.534 Gyr
Age of Universe : 13.770 Gyr
Scale (Cosmology Corrected): 349 pc/arcsec = 0.349 kpc/arcsec = 20.92 kpc/arcmin = 1.26 Mpc/degree
Surface Brightness Dimming : Flux Density per Unit Area = 0.93382; Magnitude per Unit Area = 0.07434 mag"""
)
print(
"""
# NGC 1277
[Redshift 0.016622 as corrected to the Reference Frame defined by the (Virgo + GA + Shapley)]
Luminosity Distance : 71.6 Mpc (m-M) = 34.27 mag
Angular-Size Distance : 69.3 Mpc (m-M) = 34.20 mag
Co-Moving Radial Distance : 70.4 Mpc (m-M) = 34.24 mag
Co-Moving Tangential Dist. : 70.4 Mpc (m-M) = 34.24 mag
Co-Moving Volume : 0.00146 Gpc^3
Light Travel-Time : 0.228 Gyr
Age at Redshift 0.016622 : 13.543 Gyr
Age of Universe : 13.770 Gyr
Scale (Cosmology Corrected): 336 pc/arcsec = 0.336 kpc/arcsec = 20.16 kpc/arcmin = 1.21 Mpc/degree
Surface Brightness Dimming : Flux Density per Unit Area = 0.93619; Magnitude per Unit Area = 0.07159 mag
"""
)
print(
"""
# NGC 1278
[Redshift 0.019954 as corrected to the Reference Frame defined by the (Virgo + GA + Shapley)]
Luminosity Distance : 86.2 Mpc (m-M) = 34.68 mag
Angular-Size Distance : 82.9 Mpc (m-M) = 34.59 mag
Co-Moving Radial Distance : 84.5 Mpc (m-M) = 34.63 mag
Co-Moving Tangential Dist. : 84.5 Mpc (m-M) = 34.63 mag
Co-Moving Volume : 0.00253 Gpc^3
Light Travel-Time : 0.273 Gyr
Age at Redshift 0.019954 : 13.498 Gyr
Age of Universe : 13.770 Gyr
Scale (Cosmology Corrected): 402 pc/arcsec = 0.402 kpc/arcsec = 24.10 kpc/arcmin = 1.45 Mpc/degree
Surface Brightness Dimming : Flux Density per Unit Area = 0.92401; Magnitude per Unit Area = 0.08581 mag
"""
)
distMod1, asec1 = 34.68, 0.402
distMod2, asec2 = 33.66, 0.255
scale1 = (asec1 * pix) * u.kpc
scale2 = (asec2 * pix) * u.kpc
print(scale1, scale2)
In [6]:
loc = '/Users/songhuang/work/project/perseus/final/i/'
expMod = glob.glob(loc + '*_exp.fits')
serMod = glob.glob(loc + '*_ser.fits')
print(len(expMod))
print(len(serMod))
In [40]:
expMag, expRe, expAr, expNser = [], [], [], []
expChi2, expNdof, expRchi2 = [], [], []
for modFile in expMod:
#try:
mod = gPar.GalfitResults(modFile, verbose=False)
#except Exception:
# continue
nComp = mod.num_components
for ii in range(nComp):
icomp = (ii + 1)
strComp = '_'.join(['component', str(icomp).strip()])
comp = getattr(mod, strComp)
if comp.component_type == 'sersic':
expChi2.append(mod.chisq)
expNdof.append(mod.ndof)
expRchi2.append(mod.reduced_chisq)
expMag.append(comp.mag)
expRe.append(comp.re)
expAr.append(comp.ar)
expNser.append(comp.n)
else:
continue
In [41]:
expRchi2
Out[41]: