In [263]:
%pylab inline
# Matplotlib default settings
rcdef = plt.rcParams.copy()
pylab.rcParams['figure.figsize'] = 12, 10
pylab.rcParams['xtick.major.size'] = 8.0
pylab.rcParams['xtick.major.width'] = 1.5
pylab.rcParams['xtick.minor.size'] = 4.0
pylab.rcParams['xtick.minor.width'] = 1.5
pylab.rcParams['ytick.major.size'] = 8.0
pylab.rcParams['ytick.major.width'] = 1.5
pylab.rcParams['ytick.minor.size'] = 4.0
pylab.rcParams['ytick.minor.width'] = 1.5
rc('axes', linewidth=2)

import numpy as np
from astropy.io import fits 
from __future__ import division 
from astropy import units as u

import cubehelix  # Cubehelix color scheme
import copy

import os.path
from pyraf import iraf


Populating the interactive namespace from numpy and matplotlib
WARNING: pylab import has clobbered these variables: ['copy', 'e']
`%matplotlib` prevents importing * from pylab and numpy

Parameters that are need to know:

  1. Basic parameters for the galaxy: galX, galY, galR, galQ, galPA
  2. The pixel scale, gain, exptime (expTime), and photometric zeropoint (zpPhoto) of the image
  3. For three separate ELLIPSE run, parameters like:
    • iniSMA, minSMA, maxSMA : initial, minimum, and maximum SMA of the ELLIPSE run
    • ellipStep: non-linear step-size of the ELLIPSE run
    • uppClip, lowClip, nClip, fracBadPix:
    • intMode: Mode for integration, options are bi-linear, mean, median

Optional parameters

1.

Get the 1-D Surface Brightness Profile of the Galaxy


In [264]:
# The input image, object masks, and the PSF image

imgFile = 'red_21572_Icut_ori.fits'
mskFile = 'red_21572_Icut_msk.fits'
psfFile = 'red_21572_Ipsf.fits'

In [265]:
# Read in the image 
img = (fits.open(imgFile))[0].data
msk = (fits.open(mskFile))[0].data

# Size of the image 
imgSizeX, imgSizeY = img.shape
imgSizeMax = imgSizeX if (imgSizeX >= imgSizeY) else imgSizeY

In [266]:
# Basic information about the object 
galX  = 400.5 
galY  = 400.5 
galR  = 10.0 
galQ  = 0.88
galPA = 85.0

# Basic parameters of the images 
pixelScale = 0.168  # arcsec / pixel 
gain       = 3.0 
expTime    = 1.0 
zpPhoto    = 27.0

In [537]:
# Parameters for the ELLIPSE runs 
galInnerFact  = 0.4 
galOuterFact1 = 20.0
galOuterFact2 = 25.0 
minSmaSys     = 6.0 # Minimum SMA0 = 6.0 pixel
maxSmaSysFac  = 1.2

## 1. All Free 

iniSMA1  = (galR * galInnerFact) if ((galR * galInnerFact) > minSmaSys) else minSmaSys  
minSMA1  = 0.00 
maxSMA1  = (galR * galOuterFact1) if ((galR * galOuterFact1) < (imgSizeMax * maxSmaSysFac)) else (imgSizeMax * maxSmaSysFac)
ellStep1 = 0.1
uppClip1 = 2.5
lowClip1 = 2.5 
nClip1   = 3 
fracBad1 = 0.5
intMode1 = "median"

## 2. Fixed Center, Free Geometry

iniSMA2  = (galR * galInnerFact) if ((galR * galInnerFact) > minSmaSys) else minSmaSys  
minSMA2  = 0.00 
maxSMA2  = (galR * galOuterFact1) if ((galR * galOuterFact1) < (imgSizeMax * maxSmaSysFac)) else (imgSizeMax * maxSmaSysFac)
ellStep2 = 0.08
uppClip2 = 2.0
lowClip2 = 2.0 
nClip2   = 3 
fracBad2 = 0.5
intMode2 = "median"

## 3. Fixed Center, Fixed Geometry

iniSMA3  = (galR * galInnerFact) if ((galR * galInnerFact) > minSmaSys) else minSmaSys  
minSMA3  = 0.00 
maxSMA3  = (galR * galOuterFact2) if ((galR * galOuterFact2) < (imgSizeMax * maxSmaSysFac)) else (imgSizeMax * maxSmaSysFac)
ellStep3 = 0.08
uppClip3 = 2.0
lowClip3 = 2.0 
nClip3   = 3 
fracBad3 = 0.5
intMode3 = "median"

Make a .pl version of the mask image


In [268]:
# Name of the .pl mask file for IRAF 
mskIraf = imgFile.replace('.fits', '.pl')

if os.path.isfile(mskIraf): 
    os.remove(mskIraf)
    
# Convert the fits format mask into pl format.
iraf.imcopy.input  = "'" + mskFile + "'"
iraf.imcopy.output = "'" + mskIraf + "'"
iraf.imcopy()


red_21572_Icut_msk.fits -> red_21572_Icut_ori.pl

In [269]:
# Call the STSDAS.ANALYSIS.ISOPHOTE package
iraf.stsdas()
iraf.analysis()
iraf.isophote()

Some Useful Functions


In [ ]:
def ellipseConfig(galX, galY, maxSMA, galEll=0.05, galPA=0.0, 
                  iniSMA=6.0, minSMA=0.0, ellStep=0.1, 
                  recenter=True, 

#Define parameters for the ellipse run
# 1. Initial guess of the central X, Y 
iraf.ellipse.geompar.x0         = galX
iraf.ellipse.geompar.y0         = galY
# 2. Initial guess of the ellipticity and PA of the first ISOPHOTE
iraf.ellipse.geompar.ellip0     = galEll
iraf.ellipse.geompar.pa0        = galPA
# 3. Initial radius for ellipse fitting
iraf.ellipse.geompar.sma0       = iniSMA
# 4. The minimum and maximum radius for the ellipse fitting
iraf.ellipse.geompar.minsma     = minSMA
iraf.ellipse.geompar.maxsma     = maxSMA
# 5. Parameters about the stepsize during the fitting.
iraf.ellipse.geompar.linear     = "no"
iraf.ellipse.geompar.step       = ellStep
# 6. Do you want to allow the ellipse to decide the galaxy center during the
iraf.ellipse.geompar.recenter   = "yes"
# 7. The next three parameters control the behavior of the fit
iraf.ellipse.controlpar.conver  = 2
iraf.ellipse.controlpar.hcenter = "no"
iraf.ellipse.controlpar.hellip  = "no"
iraf.ellipse.controlpar.hpa     = "no"
# 8. Parameters about the iterations
# minit/maxit: minimun and maximum number of the iterations
iraf.ellipse.controlpar.minit   = 10
iraf.ellipse.controlpar.maxit   = 200
# 9. Threshold for the object locator algorithm
iraf.ellipse.controlpar.olthresh = 1.00000
# 10. Make sure the Interactive Mode is turned off
iraf.ellipse.interactive         = "no"
# 11. Magnitude Zeropoint 
iraf.ellipse.magpar.mag0         = zpPhoto
# 12. Sampler
iraf.ellipse.samplepar.integrmode  = intMode1  # "bi-linear" "mean" "median"
iraf.ellipse.samplepar.usclip      = uppClip1 
iraf.ellipse.samplepar.lsclip      = lowClip1
iraf.ellipse.samplepar.nclip       = nClip1
iraf.ellipse.samplepar.fflag       = fracBad1
# 13. Optional Harmonics
# iraf.ellipse.samplepar.harmonics   = '1,2'

In [281]:
def ellipRemoveIndef(outTabName): 
    
    import subprocess 
    
    if os.path.exists(outTabName): 
        sub = subprocess.call(['sed', '-i_back', 's/INDEF/NaN/g', outTabName])
    else:
        raise Exception('Can not find the input catalog!')
        
    return outTabName

In [282]:
def readEllipseOut(outTabName, pix=1.0, zp=27.0, exptime=1.0, bkg=0.0,
                   harmonic=False): 
    
    name = ellipRemoveIndef(outTabName)
    from astropy.table import Table, Column
    
    ellipseOut = Table.read(outTabName, format='ascii.no_header')
    
    # Rename all the columns 
    ellipseOut.rename_column('col1',  'sma')
    ellipseOut.rename_column('col2',  'intens')
    ellipseOut.rename_column('col3',  'int_err')
    ellipseOut.rename_column('col4',  'pix_var')
    ellipseOut.rename_column('col5',  'rms')
    ellipseOut.rename_column('col6',  'ell')
    ellipseOut.rename_column('col7',  'ell_err')
    ellipseOut.rename_column('col8',  'pa')
    ellipseOut.rename_column('col9',  'pa_err')
    ellipseOut.rename_column('col10', 'x0')
    ellipseOut.rename_column('col11', 'x0_err')
    ellipseOut.rename_column('col12', 'y0')
    ellipseOut.rename_column('col13', 'y0_err')
    ellipseOut.rename_column('col14', 'grad')
    ellipseOut.rename_column('col15', 'grad_err')
    ellipseOut.rename_column('col16', 'grad_r_err')
    ellipseOut.rename_column('col17', 'rsma')
    ellipseOut.rename_column('col18', 'mag')
    ellipseOut.rename_column('col19', 'mag_lerr')
    ellipseOut.rename_column('col20', 'mag_uerr')
    ellipseOut.rename_column('col21', 'tflux_e')
    ellipseOut.rename_column('col22', 'tflux_c')
    ellipseOut.rename_column('col23', 'tmag_e')
    ellipseOut.rename_column('col24', 'tmag_c')
    ellipseOut.rename_column('col25', 'npix_e')
    ellipseOut.rename_column('col26', 'npix_c')
    ellipseOut.rename_column('col27', 'a3')
    ellipseOut.rename_column('col28', 'a3_err')
    ellipseOut.rename_column('col29', 'b3')
    ellipseOut.rename_column('col30', 'b3_err')
    ellipseOut.rename_column('col31', 'a4')
    ellipseOut.rename_column('col32', 'a4_err')
    ellipseOut.rename_column('col33', 'b4')
    ellipseOut.rename_column('col34', 'b4_err')
    ellipseOut.rename_column('col35', 'ndata')
    ellipseOut.rename_column('col36', 'nflag')
    ellipseOut.rename_column('col37', 'niter')
    ellipseOut.rename_column('col38', 'stop')
    ellipseOut.rename_column('col39', 'a_big')
    ellipseOut.rename_column('col40', 'sarea')
    if harmonic: 
        ellipseOut.rename_column('col41', 'a1')
        ellipseOut.rename_column('col42', 'a1_err')
        ellipseOut.rename_column('col43', 'a2')
        ellipseOut.rename_column('col44', 'a2_err')
    
    # Normalize the PA 
    import angles 
    ellipseOut.add_column(Column(name='pa_norm', 
                          data=np.array([angles.normalize(pa, 0.0, 180.0) for pa in ellipseOut['pa']])))
    
    # Convert the intensity into surface brightness 
    parea = (pix ** 2.0)
    ellipseOut.add_column(Column(name='sbp', 
                          data=(zp-2.5*np.log10((ellipseOut['intens']-bkg)/(parea*exptime)))))
    ellipseOut.add_column(Column(name='sbp_lerr', 
                          data=(zp-2.5*np.log10((ellipseOut['intens']-ellipseOut['int_err']-bkg)/(parea*exptime)))))  
    ellipseOut.add_column(Column(name='sbp_uerr', 
                          data=(zp-2.5*np.log10((ellipseOut['intens']+ellipseOut['int_err']-bkg)/(parea*exptime))))) 
    
    # Convert the unit of radius into arcsecs
    ellipseOut.add_column(Column(name='sma_asec', 
                                 data=(ellipseOut['sma'] * pix)))
    ellipseOut.add_column(Column(name='rsma_asec', 
                                 data=(ellipseOut['sma'] * pix) ** 0.25))

    return ellipseOut

In [96]:
def convIso2Ell(ellTab, xpad=0.0, ypad=0.0): 
    
    from matplotlib.patches import Ellipse
    
    x  = ellTab['x0'] - xpad 
    y  = ellTab['y0'] - ypad
    pa = ellTab['pa']
    a  = ellTab['sma'] * 2.0
    b  = ellTab['sma'] * 2.0 * (1.0 - ellTab['ell'])
    
    ells = [Ellipse(xy=np.array([x[i], y[i]]), 
                    width=np.array(b[i]), 
                    height=np.array(a[i]), 
                    angle=np.array(pa[i])) 
            for i in range(x.shape[0])]
    
    return ells

In [97]:
def zscale(img, contrast=0.25, samples=500):

    # Image scaling function form http://hsca.ipmu.jp/hscsphinx/scripts/psfMosaic.html
    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

In [478]:
def ellipseGetGrowthCurve(ellipOut, relThreshold=(5E-3)): 
    
    # The area in unit of pixels covered by an elliptical isophote 
    ellArea = np.pi * (ellipOut['sma'] ** 2.0 * (1.0-ellipOut['ell']))
    # The area in unit covered by the "ring"
    isoArea = np.append(ellArea[0], [ellArea[1:] - ellArea[:-1]])
    # The total flux inside the "ring"
    isoFlux = np.append(ellArea[0], [ellArea[1:] - ellArea[:-1]]) * ellipOut['intens']
    # 
    isoTFlux = np.asarray(map(lambda x: np.nansum(isoFlux[0:x+1]), range(isoFlux.shape[0])))
    # Flux difference between each "ring"
    diffIsoTFlux = np.append(isoTFlux[0], [isoTFlux[1:] - isoTFlux[:-1]])
    # Relative change of flux in each "ring"
    relDiffIso = (diffIsoTFlux / isoTFlux)
    
    # TODO: Using more examples to show whether this really works well
    
    indFluxDecrease = np.where(relDiffIso < relThreshold)
    if (indFluxDecrease[0]).size is 0:
        import warnings
        warnings.warn("WARNING!! the flux increasement is never smaller than the threshold!")
        indMaxIso = (relDiffIso.shape[0] - 1)
    else: 
        indMaxIso = (np.where(diffIsoTFlux < 0))[0][0]
    
    maxIsoSma  = ellipOut['sma'][indMaxIso]
    
    maxIsoFlux = np.nanmax(isoTFlux[0:indMaxIso])
    
    # Get the growth curve
    isoGrowthCurve = np.asarray((isoTFlux / maxIsoFlux) * 100.0)
 
    return isoGrowthCurve, maxIsoSma

In [495]:
def ellipseGetR50(ellipseRsma, isoGrowthCurve, simple=True):
    
    if len(ellipseRsma) != len(isoGrowthCurve): 
        raise "The x and y should have the same size!", len(ellipseRsma), len(isoGrowthCurve)
    else:
        if simple:
            isoRsma50 = ellipseRsma[np.nanargmin(np.abs(isoGrowthCurve - 50.0))]
        else:
            isoRsma50 = (numpy.interp([50.0], isoGrowthCurve, ellipseRsma))[0]
            
    return isoRsma50

Run 1: All Parameters Free


In [538]:
outBin1 = imgFile.replace(".fits", "_ellipse_1.bin")
outTab1 = imgFile.replace(".fits", "_ellipse_1.tab")
outCdf1 = imgFile.replace(".fits", "_ellipse_1.cdf")
  • TODO: Need better defaults of these parameters

In [539]:
#Define parameters for the ellipse run
# 1. Initial guess of the central X, Y 
iraf.ellipse.geompar.x0         = galX
iraf.ellipse.geompar.y0         = galY
# 2. Initial guess of the ellipticity and PA of the first ISOPHOTE
iraf.ellipse.geompar.ellip0     = (1.0 - galQ)
iraf.ellipse.geompar.pa0        = galPA
# 3. Initial radius for ellipse fitting
iraf.ellipse.geompar.sma0       = iniSMA1
# 4. The minimum and maximum radius for the ellipse fitting
iraf.ellipse.geompar.minsma     = minSMA1
iraf.ellipse.geompar.maxsma     = maxSMA1
# 5. Parameters about the stepsize during the fitting.
iraf.ellipse.geompar.linear     = "no"
iraf.ellipse.geompar.step       = ellStep1
# 6. Do you want to allow the ellipse to decide the galaxy center during the
iraf.ellipse.geompar.recenter   = "yes"
# 7. The next three parameters control the behavior of the fit
iraf.ellipse.controlpar.conver  = 2
iraf.ellipse.controlpar.hcenter = "no"
iraf.ellipse.controlpar.hellip  = "no"
iraf.ellipse.controlpar.hpa     = "no"
# 8. Parameters about the iterations
# minit/maxit: minimun and maximum number of the iterations
iraf.ellipse.controlpar.minit   = 10
iraf.ellipse.controlpar.maxit   = 200
# 9. Threshold for the object locator algorithm
iraf.ellipse.controlpar.olthresh = 1.00000
# 10. Make sure the Interactive Mode is turned off
iraf.ellipse.interactive         = "no"
# 11. Magnitude Zeropoint 
iraf.ellipse.magpar.mag0         = zpPhoto
# 12. Sampler
iraf.ellipse.samplepar.integrmode  = intMode1  # "bi-linear" "mean" "median"
iraf.ellipse.samplepar.usclip      = uppClip1 
iraf.ellipse.samplepar.lsclip      = lowClip1
iraf.ellipse.samplepar.nclip       = nClip1
iraf.ellipse.samplepar.fflag       = fracBad1
# 13. Optional Harmonics
# iraf.ellipse.samplepar.harmonics   = '1,2'

In [540]:
# Check and remove outputs from the previous Ellipse run
if os.path.exists(outBin1):
    os.remove(outBin1)
    
# Start the fitting of the 1st round
iraf.ellipse(input=imgFile, output=outBin1)


Running object locator... Done.
#
# Semi-    Isophote      Ellipticity     Position   Grad.  Data Flag Iter. Stop
# major      mean                         Angle      rel.                  code
# axis     intensity                                error
#(pixel)                                 (degree)
#
   6.00     3.47(  0.09) 0.118(0.007)  86.00( 1.94) 0.038   36   0    20    0
   6.60     2.95(  0.08) 0.122(0.008)  83.84( 2.02) 0.046   38   1    10    0
   7.26     2.50(  0.10) 0.112(0.011)  86.49( 2.97) 0.058   43   0    10    0
   7.99     2.14(  0.09) 0.112(0.011) -88.63( 2.95) 0.062   47   0    10    0
   8.78     1.83(  0.10) 0.121(0.012) -80.51( 3.14) 0.069   52   0    10    0
   9.66     1.53(  0.10) 0.114(0.013) -80.51( 3.51) 0.079   56   1    10    0
  10.63     1.24(  0.10) 0.091(0.014) -77.59( 4.71) 0.085   58   6    10    0
  11.69     1.00(  0.12) 0.065(0.022) -65.88( 9.57) 0.130   64   7    10    0
  12.86     0.87(  0.13) 0.098(0.026) -65.88( 7.12) 0.159   68   9    10    0
  14.15     0.73(  0.12) 0.098(0.028) -65.88( 7.59) 0.159   75   9    10    0
  15.56     0.65(  0.14) 0.153(0.032) -67.77( 5.89) 0.204   81   9    10    0
  17.12     0.53(  0.11) 0.128(0.027) -55.88( 6.03) 0.148   91   9    10    0
  18.83     0.51(  0.10) 0.203(0.019) -58.33( 2.95) 0.103  104   1    10    0
  20.71     0.39(  0.09) 0.203(0.045) -71.82( 7.13) 0.285   50   0    10    0
  22.78     0.26(  0.08) 0.032(0.099) -21.46(88.34) 0.553   58   3    10    0
  25.06     0.24(  0.08) 0.032(INDEF) -21.46(INDEF) 0.228   57   4     1    4
  27.57     0.18(  0.07) 0.032(INDEF) -21.46(INDEF) 0.289   58   3     1    4
  30.33     0.14(  0.06) 0.032(INDEF) -21.46(INDEF) 0.254   61   0     1    4
  33.36     0.10(  0.06) 0.032(INDEF) -21.46(INDEF) 0.307   63   2     1    4
  36.70     0.07(  0.05) 0.032(INDEF) -21.46(INDEF) INDEF   66   5     1    4
  40.37     0.10(  0.08) 0.032(INDEF) -21.46(INDEF) 0.428   71   7     1    4
  44.40     0.07(  0.06) 0.032(INDEF) -21.46(INDEF) 0.322   75  11     1    4
  48.84     0.06(  0.11) 0.032(INDEF) -21.46(INDEF) 0.581   86   9     1    4
  53.73     0.04(  0.06) 0.032(INDEF) -21.46(INDEF) 0.475   95   9     1    4
  59.10     0.02(  0.05) 0.032(INDEF) -21.46(INDEF) 0.479  106   8     1    4
  65.01     0.02(  0.04) 0.032(INDEF) -21.46(INDEF) 0.612  121   3     1    4
  71.51     0.01(  0.04) 0.032(INDEF) -21.46(INDEF) INDEF  116  10     1    4
  78.66     0.02(  0.06) 0.032(INDEF) -21.46(INDEF) 0.916  113  13     1    4
  86.53     0.01(  0.04) 0.032(INDEF) -21.46(INDEF) 0.546  115  11     1    4
  95.18     0.00(  0.04) 0.032(INDEF) -21.46(INDEF) INDEF  117   9     1    4
 104.70     0.00(  0.03) 0.032(INDEF) -21.46(INDEF) INDEF  114  12     1    4
 115.17     0.01(  0.03) 0.032(INDEF) -21.46(INDEF) INDEF  114  12     1    4
 126.68     0.01(  0.03) 0.032(INDEF) -21.46(INDEF) 0.916  115  11     1    4
 139.35     0.00(  0.03) 0.032(INDEF) -21.46(INDEF) INDEF  114  12     1    4
 153.29     0.00(  0.02) 0.032(INDEF) -21.46(INDEF) INDEF  117   9     1    4
 168.61    -0.00(  0.03) 0.032(INDEF) -21.46(INDEF) 2.364  115  11     1    4
 185.48    -0.00(  0.03) 0.032(INDEF) -21.46(INDEF) INDEF  116  10     1    4
   5.45     4.00(  0.07) 0.105(0.005)  86.28( 1.62) 0.033   30   3    10    0
   4.96     4.61(  0.08) 0.097(0.007)  85.83( 2.07) 0.034   30   0    10    0
   4.51     5.37(  0.07) 0.108(0.005)  81.42( 1.50) 0.032   26   1    10    0
   4.10     6.09(  0.07) 0.101(0.005)  78.77( 1.62) 0.028   25   0    10    0
   3.73     6.99(  0.06) 0.110(0.004)  79.01( 1.19) 0.025   22   0    10    0
   3.39     7.89(  0.07) 0.114(0.005)  79.65( 1.38) 0.026   20   0    10    0
   3.08     8.77(  0.05) 0.090(0.004)  77.70( 1.22) 0.029   18   1    10    0
   2.80     9.83(  0.09) 0.088(0.006)  77.51( 2.07) 0.028   17   0    10    0
   2.54    10.85(  0.09) 0.092(0.007)  79.12( 2.33) 0.031   16   0    10    0
   2.31    11.71(  0.07) 0.069(0.006)  76.14( 2.56) 0.035   14   0    10    0
   2.10    12.68(  0.08) 0.076(0.007)  78.55( 2.92) 0.042   13   0    10    0
   1.91    13.54(  0.11) 0.068(0.011)  81.50( 4.78) 0.043   13   0    10    0
   1.74    14.41(  0.15) 0.063(0.014)  82.60( 6.98) 0.059   13   0    10    0
   1.58    15.21(  0.19) 0.057(0.020)  82.75(10.67) 0.079   13   0    10    0
   1.44    15.90(  0.18) 0.050(0.022)  78.24(13.51) 0.096   13   0    10    0
   1.31    16.51(  0.16) 0.040(0.023)  80.29(17.46) 0.113   13   0    10    0
   1.19    17.07(  0.13) 0.038(0.020)  80.49(15.99) 0.102   13   0    10    0
   1.08    17.52(  0.15) 0.022(0.028)  78.13(37.10) 0.113   13   0    10    0
   0.98    18.04(  0.20) 0.079(0.047) -87.04(18.40) 0.225   13   0    10    0
   0.89    18.42(  0.27) 0.108(0.071)  89.71(20.67) 0.330   13   0    10    0
   0.81    18.72(  0.26) 0.132(0.083)  89.52(20.00) 0.404   13   0    10    0
   0.74    18.94(  0.24) 0.132(0.088)  89.52(20.99) 0.440   13   0    10    0
   0.67    19.15(  0.23) 0.133(0.089)  89.54(21.19) 0.447   13   0    10    0
   0.61    19.34(  0.21) 0.133(0.090)  89.56(21.38) 0.452   13   0    10    0
   0.55    19.52(  0.19) 0.134(0.091)  89.58(21.53) 0.458   13   0    10    0
   0.50    19.68(  0.18) 0.134(0.092)  89.58(21.69) 0.462   13   0    10    0
   0.00    21.32
   2.48  CPU seconds.
   0.03  minutes elapsed.

In [541]:
# Dump the useful information into a ASCII table
iraf.tables()
iraf.ttools()

if os.path.exists(outTab1):
    os.remove(outTab1)
if os.path.exists(outCdf1):
    os.remove(outCdf1)

iraf.tdump.columns=''
iraf.tdump(outBin1, datafil=outTab1, cdfile=outCdf1)


IMAGE    t red_21572_Icut_ori.fits

In [542]:
ellipseOut1 = readEllipseOut(outTab1, zp=zpPhoto, pix=pixelScale, exptime=expTime)

nIso1 = len(ellipseOut1)
print "%d elliptical isophotes have been extracted" % nIso1

maxRad1 = np.nanmax(ellipseOut1['sma'])
print "The maximum radius is %7.2f pixels" % maxRad1


64 elliptical isophotes have been extracted
The maximum radius is  185.48 pixels
/Users/songhuang/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:58: RuntimeWarning: invalid value encountered in log10
/Users/songhuang/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:60: RuntimeWarning: invalid value encountered in log10

In [543]:
isoGrowthCurve1, maxIsoSma1 = ellipseGetGrowthCurve(ellipseOut1)
isoRsma50_1 = ellipseGetR50(ellipseOut1['rsma'], isoGrowthCurve1)

print "The maximum iso SMA is %10.3f pixels" % maxIsoSma1 
print "The isophotal R50 is %10.3f pixels"   % (isoRsma50_1 ** 4.0)


The maximum iso SMA is    168.615 pixels
The isophotal R50 is     14.148 pixels
/Users/songhuang/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:14: RuntimeWarning: invalid value encountered in true_divide
/Users/songhuang/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:18: RuntimeWarning: invalid value encountered in less

In [544]:
fig = plt.figure(figsize=(8, 5))
fig.subplots_adjust(hspace=0.1, wspace=0.1, 
                    top=0.95, right=0.95)

ax = gca()
fontsize = 14
ax.minorticks_on()

for tick in ax.xaxis.get_major_ticks():
    tick.label1.set_fontsize(fontsize)
for tick in ax.yaxis.get_major_ticks():
    tick.label1.set_fontsize(fontsize)

ax.set_xlabel('RSMA (pixel)', fontsize=23)
ax.set_ylabel('Growth Curve', fontsize=23)

ax.axhline(100.0, linestyle='-', color='k', alpha=0.5, linewidth=2, 
           label='$f_{100}$')
ax.axhline(50.0,  linestyle='--', color='k', alpha=0.5, linewidth=2, 
           label='$f_{50}$')
ax.plot(ellipseOut1['rsma'], isoGrowthCurve1, '-', color='r', linewidth=2.0, 
        label='$growth$')
ax.axvline(maxIsoSma1**0.25, linestyle='-.', color='g', alpha=0.6, 
           linewidth=2, label='max $R$')

ax.axvline(isoRsma50_1, linestyle=':', color='c', alpha=0.6, linewidth=2, 
           label='$R_{50}$')
ax.legend(loc=[0.8, 0.05])

ax.set_title("Ellipse Run #1", fontsize=22)
ax.title.set_position((0.5,1.01))


Overplot the Ellipse on the Image


In [545]:
fig = plt.figure(figsize=(10, 10))
fig.subplots_adjust(hspace=0.1, wspace=0.1,
                    top=0.95, right=0.95)
ax = gca()
fontsize = 14
ax.minorticks_on()

for tick in ax.xaxis.get_major_ticks():
    tick.label1.set_fontsize(fontsize)
for tick in ax.yaxis.get_major_ticks():
    tick.label1.set_fontsize(fontsize)
    
ax.set_title('Elliptical Isophotes (1)', fontsize=25, fontweight='bold')
ax.title.set_position((0.5,1.01))

# NaN-Mask the image
imgMsk = copy.deepcopy(img)
imin, imax = zscale(imgMsk, contrast=0.6, samples=500)
imgMsk[msk > 0] = np.nan 

# Color map
cmap = cubehelix.cmap(start=0.5, rot=-1.0, minSat=1.2, maxSat=1.2, 
                      minLight=0., maxLight=1., gamma=0.5)
cmap.set_bad('k',1.)

# Crop the image
if (galX > maxRad1) and (galY > maxRad1): 
    zoomReg = imgMsk[np.int(galX-maxRad1):np.int(galX+maxRad1), 
                     np.int(galY-maxRad1):np.int(galY+maxRad1)]
    # Define the new center of the cropped images
    xPad = (imgSizeX / 2.0 - maxRad1)
    yPad = (imgSizeY / 2.0 - maxRad1)
else: 
    zoomReg = imgMsk
    xPad = 0 
    yPad = 0

# Show the image
ax.imshow(np.arcsinh(zoomReg), interpolation="none", 
          vmin=imin, vmax=imax, cmap=cmap)

# Get the Shapes
ellipIso1 = convIso2Ell(ellipseOut1, xpad=xPad, ypad=yPad)

# Overlay the ellipses on the image
for e in ellipIso1:
    ax.add_artist(e)
    e.set_clip_box(ax.bbox)
    e.set_alpha(0.9)
    e.set_edgecolor('r')
    e.set_facecolor('none')
    e.set_linewidth(2.0)



In [546]:
def ellipseGetAvgCen(ellipseOut, isoRsma50, minSma=0.5, ratio=1.2):

    # Get the Average X0/Y0 
    radOuter = (isoRsma50 ** 4.0) * ratio 
    
    avgCenX = np.nanmedian(ellipseOut['x0'][np.logical_and((ellipseOut['sma'] <= radOuter), 
                                                           (ellipseOut['sma'] >= minSma))])
    print ellipseOut['x0'][np.logical_and((ellipseOut['sma'] <= radOuter), 
                                                           (ellipseOut['sma'] >= minSma))]
    avgCenY = np.nanmedian(ellipseOut['y0'][np.logical_and((ellipseOut['sma'] <= radOuter), 
                                                           (ellipseOut['sma'] >= minSma))])
    
    return avgCenX, avgCenY

print "The average X0/Y0 is %7.2f/%7.2f" % ellipseGetAvgCen(ellipseOut1, isoRsma50_1)


   x0   
--------
400.9116
400.9025
400.8924
400.8813
400.8689
400.8552
400.8476
400.8465
400.8696
400.8788
     ...
401.0165
401.0188
401.0078
401.0687
401.1534
401.2739
401.4755
401.3889
401.3066
401.7171
401.8196
Length = 37 rows
The average X0/Y0 is  400.94/ 400.94

In [547]:
fig = plt.figure(figsize=(8, 5))
fig.subplots_adjust(hspace=0.1, wspace=0.1, 
                    top=0.95, right=0.95)

ax = gca()
fontsize = 14
ax.minorticks_on()
#ax.invert_yaxis()

for tick in ax.xaxis.get_major_ticks():
    tick.label1.set_fontsize(fontsize)
for tick in ax.yaxis.get_major_ticks():
    tick.label1.set_fontsize(fontsize)

ax.set_xlabel('RSMA (pixel)',      fontsize=23)
ax.set_ylabel('X0 or Y0 (pixel)', fontsize=23)

ax.errorbar(ellipseOut1['rsma'], ellipseOut1['x0'], 
            yerr=ellipseOut1['x0_err'], fmt='o', markersize=3, ecolor='r')
ax.plot(ellipseOut1['rsma'], ellipseOut1['x0'], '-', 
        color='r', linewidth=2.0, label='X0')
ax.axhline(avgCenX, linestyle='--', color='r', alpha=0.6, linewidth=3.0)

ax.errorbar(ellipseOut1['rsma'], ellipseOut1['y0'], 
            yerr=ellipseOut1['y0_err'], fmt='o', markersize=3, ecolor='b')
ax.plot(ellipseOut1['rsma'], ellipseOut1['y0'], '-', 
        color='b', linewidth=2.0, label='Y0')
ax.axhline(avgCenY, linestyle='--', color='b', alpha=0.6, linewidth=3.0)

# Put legend on
ax.legend(loc=[0.06, 0.08], fontsize=20)


Out[547]:
<matplotlib.legend.Legend at 0x1400dad50>

Run 2: Fixed X0/Y0, Free Ellipticity/PA


In [115]:
outBin2 = imgFile.replace(".fits", "_ellipse_2.bin")
outTab2 = imgFile.replace(".fits", "_ellipse_2.tab")
outCdf2 = imgFile.replace(".fits", "_ellipse_2.cdf")

In [125]:
#Define parameters for the ellipse run
# 1. Initial guess of the central X, Y 
iraf.ellipse.geompar.x0         = avgCenX
iraf.ellipse.geompar.y0         = avgCenY
# 2. Initial guess of the ellipticity and PA of the first ISOPHOTE
iraf.ellipse.geompar.ellip0     = (1.0 - galQ)
iraf.ellipse.geompar.pa0        = galPA
# 3. Initial radius for ellipse fitting
iraf.ellipse.geompar.sma0       = iniSMA2
# 4. The minimum and maximum radius for the ellipse fitting
iraf.ellipse.geompar.minsma     = minSMA2
iraf.ellipse.geompar.maxsma     = maxSMA2
# 5. Parameters about the stepsize during the fitting.
iraf.ellipse.geompar.linear     = "no"
iraf.ellipse.geompar.step       = ellStep2
# 6. Do you want to allow the ellipse to decide the galaxy center during the
iraf.ellipse.geompar.recenter   = "no"
# 7. The next three parameters control the behavior of the fit
iraf.ellipse.controlpar.conver  = 2
iraf.ellipse.controlpar.hcenter = "yes"
iraf.ellipse.controlpar.hellip  = "no"
iraf.ellipse.controlpar.hpa     = "no"
# 8. Parameters about the iterations
# minit/maxit: minimun and maximum number of the iterations
iraf.ellipse.controlpar.minit   = 10
iraf.ellipse.controlpar.maxit   = 200
# 9. Threshold for the object locator algorithm
iraf.ellipse.controlpar.olthresh = 1.00000
# 10. Make sure the Interactive Mode is turned off
iraf.ellipse.interactive         = "no"
# 11. Magnitude Zeropoint 
iraf.ellipse.magpar.mag0         = zpPhoto
# 12. Sampler
iraf.ellipse.samplepar.integrmode  = intMode2  # "bi-linear" "mean" "median"
iraf.ellipse.samplepar.usclip      = uppClip2
iraf.ellipse.samplepar.lsclip      = lowClip2
iraf.ellipse.samplepar.nclip       = nClip2
iraf.ellipse.samplepar.fflag       = fracBad2
# 13. Optional Harmonics
iraf.ellipse.samplepar.harmonics   = '1,2'

In [126]:
# Check and remove outputs from the previous Ellipse run
if os.path.exists(outBin2):
    os.remove(outBin2)
    
# Start the fitting of the 1st round
iraf.ellipse(input=imgFile, output=outBin2)


Running object locator... Done.
#
# Semi-    Isophote      Ellipticity     Position   Grad.  Data Flag Iter. Stop
# major      mean                         Angle      rel.                  code
# axis     intensity                                error
#(pixel)                                 (degree)
#
   6.00     3.47(  0.15) 0.118(0.007)  86.32( 1.78) 0.064   36   0    20    0
   6.60     2.94(  0.14) 0.119(0.009)  84.23( 2.25) 0.074   39   0    10    0
   7.26     2.49(  0.15) 0.113(0.011)  86.02( 2.87) 0.092   43   0    10    0
   7.99     2.13(  0.16) 0.111(0.010) -89.17( 2.81) 0.092   47   0    10    0
   8.78     1.82(  0.14) 0.115(0.010) -81.71( 2.77) 0.089   52   0    10    0
   9.66     1.49(  0.15) 0.107(0.015) -75.20( 4.36) 0.139   57   0    10    0
  10.63     1.26(  0.20) 0.104(0.022) -71.26( 6.36) 0.175   63   0    10    0
  11.69     1.06(  0.19) 0.103(0.022) -73.23( 6.34) 0.178   64   5    10    0
  12.86     0.88(  0.16) 0.096(0.026) -65.11( 7.44) 0.219   70   7    10    0
  14.15     0.73(  0.18) 0.102(0.034) -63.69( 8.89) 0.239   75   9    10    0
  15.56     0.64(  0.14) 0.132(0.034) -73.87( 7.08) 0.245   81  10    10    0
  17.12     0.52(  0.11) 0.112(0.025) -54.53( 6.48) 0.139   91  10    10    0
  18.83     0.47(  0.09) 0.164(0.022) -63.44( 4.13) 0.140  101   7    10    0
  20.71     0.42(  0.11) 0.206(0.033) -65.11( 5.12) 0.227   50   0    10    0
  22.78     0.32(  0.08) 0.229(0.040) -71.98( 5.74) 0.274   46   2    10    0
  25.06     0.23(  0.07) 0.090(0.064) -75.37(21.57) 0.302   56   1    10    0
  27.57     0.21(  0.06) 0.167(0.046) -80.95( 8.57) 0.258   52   0    10    0
  30.33     0.16(  0.05) 0.160(0.054) -81.37(11.22) 0.364   50   3    10    0
  33.36     0.12(  0.05) 0.156(0.043) -83.27( 8.74) 0.262   52   4    10    0
  36.70     0.09(  0.05) 0.096(0.219)  49.45(68.62) 1.693   63   3    10    0
  40.37     0.08(  0.07) 0.096(INDEF)  49.45(INDEF) 1.153   67   6     1    4
  44.40     0.07(  0.08) 0.096(INDEF)  49.45(INDEF) 0.462   73   7     1    4
  48.84     0.05(  0.07) 0.096(INDEF)  49.45(INDEF) 0.494   79   9     1    4
  53.73     0.04(  0.06) 0.096(INDEF)  49.45(INDEF) 0.730   86  11     1    4
  59.10     0.03(  0.05) 0.096(INDEF)  49.45(INDEF) 0.794   96  11     1    4
  65.01     0.02(  0.04) 0.096(INDEF)  49.45(INDEF) INDEF  110   7     1    4
  71.51     0.02(  0.04) 0.096(INDEF)  49.45(INDEF) 0.669  114   9     1    4
  78.66     0.02(  0.06) 0.096(INDEF)  49.45(INDEF) 0.486  113  13     1    4
  86.53     0.01(  0.03) 0.096(INDEF)  49.45(INDEF) 0.411  113  13     1    4
  95.18     0.00(  0.04) 0.096(INDEF)  49.45(INDEF) 0.783  116  10     1    4
 104.70    -0.00(  0.03) 0.096(INDEF)  49.45(INDEF) INDEF  117   9     1    4
 115.17     0.00(  0.03) 0.096(INDEF)  49.45(INDEF) INDEF  115  11     1    4
 126.68     0.01(  0.03) 0.096(INDEF)  49.45(INDEF) INDEF  112  14     1    4
 139.35     0.01(  0.03) 0.096(INDEF)  49.45(INDEF) 1.172  114  12     1    4
 153.29     0.01(  0.03) 0.096(INDEF)  49.45(INDEF) 0.585  118   8     1    4
 168.61     0.00(  0.03) 0.096(INDEF)  49.45(INDEF) 1.069  114  12     1    4
 185.48    -0.00(  0.03) 0.096(INDEF)  49.45(INDEF) INDEF  115  11     1    4
   5.45     4.02(  0.15) 0.107(0.007)  88.06( 2.07) 0.060   33   0    10    0
   4.96     4.61(  0.13) 0.097(0.007)  85.91( 2.12) 0.055   30   0    10    0
   4.51     5.33(  0.12) 0.104(0.006)  82.36( 1.89) 0.048   27   0    10    0
   4.10     6.09(  0.10) 0.101(0.005)  78.87( 1.61) 0.038   25   0    10    0
   3.73     6.99(  0.09) 0.109(0.004)  79.12( 1.20) 0.033   22   0    10    0
   3.39     7.89(  0.08) 0.114(0.005)  79.71( 1.42) 0.030   20   0    10    0
   3.08     8.77(  0.07) 0.092(0.004)  77.88( 1.50) 0.028   19   0    10    0
   2.80     9.83(  0.10) 0.089(0.006)  77.51( 2.11) 0.027   17   0    10    0
   2.54    10.86(  0.10) 0.093(0.007)  79.62( 2.32) 0.031   16   0    10    0
   2.31    11.71(  0.07) 0.068(0.006)  76.13( 2.57) 0.036   14   0    10    0
   2.10    12.70(  0.15) 0.077(0.008)  78.18( 2.97) 0.049   13   0    10    0
   1.91    13.56(  0.18) 0.070(0.011)  81.42( 4.77) 0.070   13   0    10    0
   1.74    14.43(  0.24) 0.064(0.015)  82.08( 6.92) 0.091   13   0    10    0
   1.58    15.22(  0.29) 0.057(0.019)  79.38(10.43) 0.123   13   0    10    0
   1.44    15.95(  0.34) 0.052(0.021)  83.99(12.34) 0.169   13   0    10    0
   1.31    16.56(  0.32) 0.048(0.021)  78.66(13.25) 0.190   13   0    10    0
   1.19    17.13(  0.35) 0.051(0.023)  81.96(13.66) 0.234   13   0    10    0
   1.08    17.68(  0.38) 0.074(0.037)  84.16(15.33) 0.315   13   0    10    0
   0.98    18.08(  0.40) 0.098(0.055)  84.41(17.67) 0.435   13   0    10    0
   0.89    18.37(  0.38) 0.102(0.067)  85.63(20.51) 0.535   13   0    10    0
   0.81    18.62(  0.36) 0.099(0.070)  86.28(22.27) 0.568   13   0    10    0
   0.74    18.85(  0.33) 0.095(0.072)  86.89(23.46) 0.573   13   0    10    0
   0.67    19.06(  0.31) 0.092(0.074)  87.47(24.84) 0.582   13   0    10    0
   0.61    19.26(  0.29) 0.089(0.076)  88.03(26.41) 0.596   13   0    10    0
   0.55    19.44(  0.28) 0.085(0.078)  88.56(28.25) 0.615   13   0    10    0
   0.50    19.60(  0.26) 0.082(0.081)  89.07(30.38) 0.640   13   0    10    0
   0.00    21.32
   2.88  CPU seconds.
   0.03  minutes elapsed.

In [127]:
# Dump the useful information into a ASCII table
iraf.tables()
iraf.ttools()

if os.path.exists(outTab2):
    os.remove(outTab2)
if os.path.exists(outCdf2):
    os.remove(outCdf2)

iraf.tdump.columns=''
iraf.tdump(outBin2, datafil=outTab2, cdfile=outCdf2)


IMAGE    t red_21572_Icut_ori.fits

In [138]:
ellipseOut2 = readEllipseOut(outTab2, zp=zpPhoto, pix=pixelScale, 
                             exptime=expTime, harmonic=True)

nIso2 = len(ellipseOut2)
print "%d elliptical isophotes have been extracted" % nIs2

maxRad2 = np.nanmax(ellipseOut2['sma'])
print "The maximum radius is %7.2f pixels" % maxRad2


64 elliptical isophotes have been extracted
The maximum radius is  185.48 pixels
/Users/songhuang/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:58: RuntimeWarning: invalid value encountered in log10
/Users/songhuang/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:60: RuntimeWarning: invalid value encountered in log10

Overplot the Ellipse on the Image


In [139]:
fig = plt.figure(figsize=(10, 10))
fig.subplots_adjust(hspace=0.1, wspace=0.1,
                    top=0.95, right=0.95)
ax = gca()
fontsize = 14
ax.minorticks_on()

for tick in ax.xaxis.get_major_ticks():
    tick.label1.set_fontsize(fontsize)
for tick in ax.yaxis.get_major_ticks():
    tick.label1.set_fontsize(fontsize)
    
ax.set_title('Elliptical Isophotes (2)', fontsize=25, fontweight='bold')
ax.title.set_position((0.5,1.01))

# NaN-Mask the image
imgMsk = copy.deepcopy(img)
imin, imax = zscale(imgMsk, contrast=0.6, samples=500)
imgMsk[msk > 0] = np.nan 

# Color map
cmap = cubehelix.cmap(start=0.5, rot=-1.0, minSat=1.2, maxSat=1.2, 
                      minLight=0., maxLight=1., gamma=0.5)
cmap.set_bad('k',1.)

# Crop the image
if (galX > maxRad2) and (galY > maxRad2): 
    zoomReg = imgMsk[np.int(galX-maxRad2):np.int(galX+maxRad2), 
                     np.int(galY-maxRad2):np.int(galY+maxRad2)]
    # Define the new center of the cropped images
    xPad = (imgSizeX / 2.0 - maxRad2)
    yPad = (imgSizeY / 2.0 - maxRad2)
else: 
    zoomReg = imgMsk
    xPad = 0 
    yPad = 0

# Show the image
ax.imshow(np.arcsinh(zoomReg), interpolation="none", 
          vmin=imin, vmax=imax, cmap=cmap)

# Get the Shapes
ellipIso2 = convIso2Ell(ellipseOut2, xpad=xPad, ypad=yPad)

# Overlay the ellipses on the image
for e in ellipIso2:
    ax.add_artist(e)
    e.set_clip_box(ax.bbox)
    e.set_alpha(0.9)
    e.set_edgecolor('r')
    e.set_facecolor('none')
    e.set_linewidth(2.0)