notebook.community

Edit and run 





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)



    


















    
























In [174]:



    


# Get the Average Ellipticity and PA 
radInner =  4.0  # Need to have a way to define this limit
radOuter = 20.0

avgEll = np.nanmedian(ellipseOut2['ell'][np.logical_and((ellipseOut2['sma'] >= radInner), 
                                                        (ellipseOut2['sma'] <= radOuter))])
avgPA  = np.nanmedian(ellipseOut2['pa_norm'][np.logical_and((ellipseOut2['sma'] >= radInner), 
                                                            (ellipseOut2['sma'] <= radOuter))])

print "The median ellipticity of the inner region is %5.3f"    % avgEll 
print "The median position angle of the inner region is %7.2f" % avgPA



    


















    






The median ellipticity of the inner region is 0.107
The median position angle of the inner region is   98.29

















In [177]:



    


# Plot the elliptical and PA profiles

fig = plt.figure(figsize=(8, 8))
fig.subplots_adjust(hspace=0.1, wspace=0.1, 
                    top=0.95, right=0.95)

rec1 = [0.1, 0.1, 0.88, 0.4]
rec2 = [0.1, 0.5, 0.88, 0.4]

ax1 = plt.axes(rec1)
ax2 = plt.axes(rec2)

# Ellipticity 
ax1.tick_params(axis='both', which='major', labelsize=20)
ax1.set_xlabel('RSMA (pixel)', fontsize=23)
ax1.set_ylabel('Ellipticity',  fontsize=23)
ax1.set_ylim(-0.05, np.nanmax(ellipseOut2['ell'] + ellipseOut2['ell_err']))

ax1.axhline(avgEll, color='k', linestyle='--', linewidth=2)
ax1.errorbar(ellipseOut2['rsma'], ellipseOut2['ell'], 
            yerr=ellipseOut2['ell_err'], fmt='o', markersize=3, ecolor='r')
ax1.plot(ellipseOut2['rsma'], ellipseOut2['ell'], '-', 
        color='r', linewidth=2.0)
ax1.minorticks_on()


# Position Angle

ax2.tick_params(axis='both', which='major', labelsize=20)
ax2.set_ylabel('PA (degree)',  fontsize=23)
#ax2.set_ylim(-0.05, np.nanmax(ellipseOut2['ell'] + ellipseOut2['ell_err']))

ax2.axhline(avgPA, color='k', linestyle='--', linewidth=2)
ax2.errorbar(ellipseOut2['rsma'], ellipseOut2['pa_norm'], 
            yerr=ellipseOut2['ell_err'], fmt='o', markersize=3, ecolor='r')
ax2.plot(ellipseOut2['rsma'], ellipseOut2['pa_norm'], '-', 
        color='r', linewidth=2.0)
ax2.xaxis.set_major_formatter(NullFormatter())
ax2.minorticks_on()



    


















    
























In [198]:



    


# Plot the A3/B3 and A4/B4

fig = plt.figure(figsize=(8, 8))
fig.subplots_adjust(hspace=0.1, wspace=0.1, 
                    top=0.95, right=0.95)

rec1 = [0.1, 0.1, 0.88, 0.4]
rec2 = [0.1, 0.5, 0.88, 0.4]

ax1 = plt.axes(rec1)
ax2 = plt.axes(rec2)

# Ellipticity 
ax1.tick_params(axis='both', which='major', labelsize=20)
ax1.set_xlabel('RSMA (pixel)', fontsize=23)
ax1.set_ylabel('A3/B3',  fontsize=23)
#ax1.set_ylim(-0.05, np.nanmax(ellipseOut2['b3'] + ellipseOut2['ell_err']))

ax1.errorbar(ellipseOut2['rsma'], ellipseOut2['a3'], 
            yerr=ellipseOut2['a3_err'], fmt='o', markersize=3, ecolor='r')
ax1.plot(ellipseOut2['rsma'], ellipseOut2['a3'], '-', 
        color='r', linewidth=2.0)
ax1.errorbar(ellipseOut2['rsma'], ellipseOut2['b3'], 
            yerr=ellipseOut2['b3_err'], fmt='o', markersize=3, ecolor='b')
ax1.plot(ellipseOut2['rsma'], ellipseOut2['b3'], '-', 
        color='b', linewidth=2.0)
ax1.minorticks_on()


# Position Angle

ax2.tick_params(axis='both', which='major', labelsize=20)
ax2.set_ylabel('A4 / B4',  fontsize=23)

ax2.errorbar(ellipseOut2['rsma'], ellipseOut2['a4'], 
            yerr=ellipseOut2['a4_err'], fmt='o', markersize=3, ecolor='r')
ax2.plot(ellipseOut2['rsma'], ellipseOut2['a4'], '-', 
        color='r', linewidth=2.0)
ax2.errorbar(ellipseOut2['rsma'], ellipseOut2['b4'], 
            yerr=ellipseOut2['b4_err'], fmt='o', markersize=3, ecolor='b')
ax2.plot(ellipseOut2['rsma'], ellipseOut2['b4'], '-', 
        color='b', linewidth=2.0)

ax2.xaxis.set_major_formatter(NullFormatter())
ax2.minorticks_on()

Run 3: Fixed X0/Y0, Ellipticity, and PA





In [178]:



    


outBin3 = imgFile.replace(".fits", "_ellipse_3.bin")
outTab3 = imgFile.replace(".fits", "_ellipse_3.tab")
outCdf3 = imgFile.replace(".fits", "_ellipse_3.cdf")



    











In [299]:



    


# The range for PA allowed by IRAF is [-90, 90]
import angles 
avgPaIraf = angles.normalize(avgPA, -90.0, 90.0)

#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     = avgEll
iraf.ellipse.geompar.pa0        = avgPaIraf
# 3. Initial radius for ellipse fitting
iraf.ellipse.geompar.sma0       = iniSMA3
# 4. The minimum and maximum radius for the ellipse fitting
iraf.ellipse.geompar.minsma     = minSMA3
#iraf.ellipse.geompar.maxsma     = maxSMA3
iraf.ellipse.geompar.maxsma     = 800
# 5. Parameters about the stepsize during the fitting.
iraf.ellipse.geompar.linear     = "no"
iraf.ellipse.geompar.step       = ellStep3
# 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  = "yes"
iraf.ellipse.controlpar.hpa     = "yes"
# 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  = intMode3  # "bi-linear" "mean" "median"
iraf.ellipse.samplepar.usclip      = uppClip3 
iraf.ellipse.samplepar.lsclip      = lowClip3
iraf.ellipse.samplepar.nclip       = nClip3
iraf.ellipse.samplepar.fflag       = fracBad3



    











In [300]:



    


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



    


















    






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.42(  0.14) 0.107(0.006) -81.71( 1.90) 0.064   32   4     1    0
   6.48     2.96(  0.13) 0.107(0.008) -81.71( 2.38) 0.079   35   4     1    0
   7.00     2.62(  0.15) 0.107(0.011) -81.71( 3.19) 0.096   34   8     1    0
   7.56     2.32(  0.15) 0.107(0.010) -81.71( 2.82) 0.091   44   1     1    0
   8.16     2.05(  0.16) 0.107(0.010) -81.71( 2.97) 0.091   49   0     1    0
   8.82     1.79(  0.14) 0.107(0.011) -81.71( 3.09) 0.093   52   0     1    0
   9.52     1.54(  0.15) 0.107(0.013) -81.71( 3.60) 0.107   57   0     1    0
  10.28     1.31(  0.15) 0.107(0.013) -81.71( 3.66) 0.140   56   5     1    0
  11.11     1.14(  0.16) 0.107(0.017) -81.71( 4.70) 0.170   58   8     1    0
  11.99     1.00(  0.17) 0.107(0.021) -81.71( 5.59) 0.183   60  11     1    0
  12.95     0.87(  0.16) 0.107(0.020) -81.71( 5.19) 0.170   64  13     1    0
  13.99     0.74(  0.15) 0.107(0.029) -81.71( 7.31) 0.224   71  12     1    0
  15.11     0.66(  0.16) 0.107(0.032) -81.71( 8.30) 0.217   79  11     1    0
  16.32     0.57(  0.12) 0.107(0.032) -81.71( 8.44) 0.166   86  11     1    0
  17.62     0.49(  0.10) 0.107(0.023) -81.71( 6.05) 0.133   92  12     1    0
  19.03     0.42(  0.07) 0.107(0.024) -81.71( 6.31) 0.141   95  18     1    0
  20.56     0.34(  0.08) 0.107(0.029) -81.71( 8.13) 0.323  106  16     1    0
  22.20     0.32(  0.09) 0.107(0.024) -81.71( 7.01) 0.134  125   6     1    0
  23.98     0.25(  0.07) 0.107(0.041) -81.71(11.18) 0.407  127  15     1    0
  25.89     0.22(  0.07) 0.107(0.051) -81.71(15.34) 0.296   62   8     1    0
  27.97     0.19(  0.07) 0.107(0.044) -81.71(12.79) 0.270   67   3     1    0
  30.20     0.14(  0.07) 0.107(0.067) -81.71(19.68) 0.440   63   7     1    0
  32.62     0.13(  0.09) 0.107(0.057) -81.71(16.34) 0.264   67   3     1    0
  35.23     0.10(  0.04) 0.107(0.046) -81.71(13.86) 0.296   62   8     1    0
  38.05     0.09(  0.05) 0.107(0.113) -81.71(31.83) 0.708   64   6     1    0
  41.09     0.08(  0.07) 0.107(INDEF) -81.71(INDEF) 0.656   65   9     1    4
  44.38     0.06(  0.05) 0.107(INDEF) -81.71(INDEF) 0.966   71   9     1    4
  47.93     0.06(  0.06) 0.107(INDEF) -81.71(INDEF) 0.276   75  11     1    4
  51.76     0.03(  0.06) 0.107(INDEF) -81.71(INDEF) INDEF   82  11     1    4
  55.90     0.04(  0.07) 0.107(INDEF) -81.71(INDEF) 0.474   89  12     1    4
  60.38     0.02(  0.04) 0.107(INDEF) -81.71(INDEF) 0.496   95  14     1    4
  65.21     0.02(  0.04) 0.107(INDEF) -81.71(INDEF) INDEF   99  17     1    4
  70.42     0.01(  0.04) 0.107(INDEF) -81.71(INDEF) INDEF  105  17     1    4
  76.06     0.02(  0.05) 0.107(INDEF) -81.71(INDEF) 0.537  108  17     1    4
  82.14     0.03(  0.10) 0.107(INDEF) -81.71(INDEF) 0.401  109  17     1    4
  88.71     0.01(  0.05) 0.107(INDEF) -81.71(INDEF) 0.449  108  18     1    4
  95.81     0.00(  0.03) 0.107(INDEF) -81.71(INDEF) 0.536  113  13     1    4
 103.47    -0.00(  0.03) 0.107(INDEF) -81.71(INDEF) INDEF  109  17     1    4
 111.75     0.00(  0.03) 0.107(INDEF) -81.71(INDEF) 1.529  113  13     1    4
 120.69     0.00(  0.03) 0.107(INDEF) -81.71(INDEF) INDEF  113  13     1    4
 130.35     0.01(  0.03) 0.107(INDEF) -81.71(INDEF) 0.943  113  13     1    4
 140.78     0.00(  0.03) 0.107(INDEF) -81.71(INDEF) 2.159  112  14     1    4
 152.04    -0.00(  0.02) 0.107(INDEF) -81.71(INDEF) 2.104  109  17     1    4
 164.20     0.00(  0.02) 0.107(INDEF) -81.71(INDEF) 1.420  113  13     1    4
 177.34    -0.01(  0.02) 0.107(INDEF) -81.71(INDEF) INDEF  109  17     1    4
 191.52     0.00(  0.02) 0.107(INDEF) -81.71(INDEF) INDEF  111  15     1    4
 206.84     0.00(  0.04) 0.107(INDEF) -81.71(INDEF) 2.731  115  11     1    4
 223.39     0.00(  0.02) 0.107(INDEF) -81.71(INDEF) 1.165  111  15     1    4
 241.26    -0.00(  0.02) 0.107(INDEF) -81.71(INDEF) 1.977  104  22     1    4
 260.57    -0.00(  0.02) 0.107(INDEF) -81.71(INDEF) 0.998  108  18     1    4
 281.41    -0.01(  0.02) 0.107(INDEF) -81.71(INDEF) INDEF  106  20     1    4
 303.92    -0.00(  0.02) 0.107(INDEF) -81.71(INDEF) INDEF  109  17     1    4
 328.24     0.00(  0.04) 0.107(INDEF) -81.71(INDEF) 0.994  114  12     1    4
 354.50    -0.00(  0.03) 0.107(INDEF) -81.71(INDEF) INDEF  114  12     1    4
 382.86     0.00(  0.02) 0.107(INDEF) -81.71(INDEF) INDEF  108  18     1    4
 413.48     0.00(  0.02) 0.107(INDEF) -81.71(INDEF) 0.669   82  44     1    4
 446.56     0.00(  0.02) 0.107(0.929) -81.71(262.2) 1.810   47  79     1    1
   5.56     3.88(  0.16) 0.107(0.008) -81.71( 2.26) 0.068   31   2     1    0
   5.14     4.38(  0.15) 0.107(0.008) -81.71( 2.21) 0.068   29   2     1    0
   4.76     4.91(  0.15) 0.107(0.006) -81.71( 1.95) 0.066   27   2     1    0
   4.41     5.49(  0.18) 0.107(0.006) -81.71( 1.72) 0.066   24   2     1    0
   4.08     6.14(  0.24) 0.107(0.005) -81.71( 1.41) 0.071   25   0     1    0
   3.78     6.83(  0.26) 0.107(0.004) -81.71( 1.20) 0.083   23   0     1    0
   3.50     7.54(  0.28) 0.107(0.004) -81.71( 1.13) 0.095   21   0     1    0
   3.24     8.32(  0.29) 0.107(0.004) -81.71( 1.29) 0.090   20   0     1    0
   3.00     9.16(  0.32) 0.107(0.005) -81.71( 1.40) 0.097   18   0     1    0
   2.78    10.00(  0.31) 0.107(0.006) -81.71( 1.83) 0.096   17   0     1    0
   2.57    10.82(  0.28) 0.107(0.006) -81.71( 1.82) 0.096   16   0     1    0
   2.38    11.58(  0.27) 0.107(0.006) -81.71( 1.85) 0.102   15   0     1    0
   2.21    12.38(  0.31) 0.107(0.007) -81.71( 2.07) 0.112   14   0     1    0
   2.04    13.17(  0.29) 0.107(0.010) -81.71( 2.91) 0.124   13   0     1    0
   1.89    13.92(  0.30) 0.107(0.012) -81.71( 3.39) 0.124   13   0     1    0
   1.75    14.61(  0.31) 0.107(0.013) -81.71( 3.85) 0.138   13   0     1    0
   1.62    15.24(  0.32) 0.107(0.013) -81.71( 3.79) 0.160   13   0     1    0
   1.50    15.82(  0.34) 0.107(0.013) -81.71( 3.81) 0.177   13   0     1    0
   1.39    16.36(  0.34) 0.107(0.013) -81.71( 3.86) 0.196   13   0     1    0
   1.29    16.84(  0.35) 0.107(0.017) -81.71( 5.07) 0.225   13   0     1    0
   1.19    17.28(  0.36) 0.107(0.026) -81.71( 7.54) 0.252   13   0     1    0
   1.10    17.65(  0.36) 0.107(0.036) -81.71(10.37) 0.302   13   0     1    0
   1.02    17.94(  0.36) 0.107(0.050) -81.71(14.74) 0.377   13   0     1    0
   0.95    18.19(  0.37) 0.107(0.063) -81.71(18.35) 0.447   13   0     1    0
   0.88    18.41(  0.34) 0.107(0.067) -81.71(19.67) 0.496   13   0     1    0
   0.81    18.62(  0.32) 0.107(0.069) -81.71(20.14) 0.513   13   0     1    0
   0.75    18.81(  0.30) 0.107(0.069) -81.71(20.25) 0.516   13   0     1    0
   0.70    18.99(  0.29) 0.107(0.070) -81.71(20.35) 0.521   13   0     1    0
   0.64    19.15(  0.27) 0.107(0.070) -81.71(20.43) 0.529   13   0     1    0
   0.60    19.31(  0.26) 0.107(0.070) -81.71(20.52) 0.539   13   0     1    0
   0.55    19.45(  0.25) 0.107(0.071) -81.71(20.61) 0.552   13   0     1    0
   0.51    19.59(  0.24) 0.107(0.071) -81.71(20.72) 0.569   13   0     1    0
   0.00    21.32
  13.40  CPU seconds.
   0.15  minutes elapsed.

















In [308]:



    


# Dump the useful information into a ASCII table
iraf.tables()
iraf.ttools()

if os.path.exists(outTab3):
    os.remove(outTab3)
if os.path.exists(outCdf3):
    os.remove(outCdf3)

iraf.tdump.columns=''
iraf.tdump(outBin3, datafil=outTab3, cdfile=outCdf3)



    


















    






IMAGE    t red_21572_Icut_ori.fits

















In [312]:



    


ellipseOut3 = readEllipseOut(outTab3, zp=zpPhoto, pix=pixelScale, exptime=expTime)

nIso3 = len(ellipseOut3)
print "%d elliptical isophotes have been extracted" % nIso3

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



    


















    






90 elliptical isophotes have been extracted
The maximum radius is  446.56 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
/Users/songhuang/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:62: RuntimeWarning: invalid value encountered in log10

















In [315]:



    


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 (3)', 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 > maxRad3) and (galY > maxRad3): 
    zoomReg = imgMsk[np.int(galX-maxRad3):np.int(galX+maxRad3), 
                     np.int(galY-maxRad3):np.int(galY+maxRad3)]
    # Define the new center of the cropped images
    xPad = (imgSizeX / 2.0 - maxRad3)
    yPad = (imgSizeY / 2.0 - maxRad3)
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
ellipIso3 = convIso2Ell(ellipseOut3, xpad=xPad, ypad=yPad)

# Overlay the ellipses on the image
for e in ellipIso3:
    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 [505]:



    


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('${\mu}_{i}$ (mag/arcsec$^2$)', fontsize=23)

# Getting the growth curve
isoGrowthCurve, maxIsoRad = ellipseGetGrowthCurve(ellipseOut3, relThreshold=(5E-3))

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(ellipseOut3['rsma'], isoGrowthCurve, '-', color='r', linewidth=2.0, 
        label='$growth$')
ax.axvline(maxIsoRad**0.25, linestyle='-.', color='g', alpha=0.6, 
           linewidth=2, label='max $R$')

isoRsma50 = ellipseGetR50(ellipseOut3['rsma'], isoGrowthCurve, simple=False)
ax.axvline(isoRsma50, linestyle=':', color='c', alpha=0.6, linewidth=2, 
           label='$R_{50}$')
ax.legend(loc=[0.8, 0.05])



    


















    






/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










    
Out[505]:








<matplotlib.legend.Legend at 0x13262fe90>

4 1-D Surface Brightness Profile of PSF





In [213]:



    


# Files for PSF profiles
psfBin = psfFile.replace(".fits", "_ellip.bin")
psfTab = psfFile.replace(".fits", "_ellip.tab")
psfCdf = psfFile.replace(".fits", "_ellip.cdf")



    











In [214]:



    


# Read in the PSF image, get the center of PSF
psf = (fits.open(psfFile))[0].data
psfCenX, psfCenY = (psf.shape[0] / 2.0), (psf.shape[1] / 2.0)
psfMaxRad = psfCenX * np.sqrt(2.0)



    











In [215]:



    


#Define parameters for the ellipse run
# 1. Initial guess of the central X, Y 
iraf.ellipse.geompar.x0         = psfCenX
iraf.ellipse.geompar.y0         = psfCenY
# 2. Initial guess of the ellipticity and PA of the first ISOPHOTE
iraf.ellipse.geompar.ellip0     = 0.05
iraf.ellipse.geompar.pa0        = 0.00
# 3. Initial radius for ellipse fitting
iraf.ellipse.geompar.sma0       = 6.0
# 4. The minimum and maximum radius for the ellipse fitting
iraf.ellipse.geompar.minsma     = 0
iraf.ellipse.geompar.maxsma     = psfMaxRad
# 5. Parameters about the stepsize during the fitting.
iraf.ellipse.geompar.linear     = "no"
iraf.ellipse.geompar.step       = 0.08
# 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 = "yes"
iraf.ellipse.controlpar.hellip  = "yes"
iraf.ellipse.controlpar.hpa     = "yes"
# 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  = "median"  # "bi-linear" "mean"
iraf.ellipse.samplepar.usclip      = 2.5 
iraf.ellipse.samplepar.lsclip      = 2.5
iraf.ellipse.samplepar.nclip       = 3
iraf.ellipse.samplepar.fflag       = 0.50



    











In [216]:



    


# Check and remove outputs from the previous Ellipse run
if os.path.exists(psfBin):
    os.remove(psfBin)
    
# Start the fitting of the 1st round
iraf.ellipse(input=psfFile, output=psfBin)



    


















    






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     0.00(  0.00) 0.050(0.000)   0.00( 0.04) INDEF   37   0     1    0
   6.48     0.00(  0.00) 0.050(0.000)   0.00( 0.03) INDEF   40   0     1    0
   7.00     0.00(  0.00) 0.050(INDEF)   0.00(INDEF) INDEF   43   0     1    4
   7.56     0.00(  0.00) 0.050(INDEF)   0.00(INDEF) INDEF   46   0     1    4
   8.16     0.00(  0.00) 0.050(INDEF)   0.00(INDEF) INDEF   50   0     1    4
   8.82     0.00(  0.00) 0.050(INDEF)   0.00(INDEF) INDEF   54   0     1    4
   9.52     0.00(  0.00) 0.050(INDEF)   0.00(INDEF) INDEF   57   1     1    4
  10.28     0.00(  0.00) 0.050(INDEF)   0.00(INDEF) INDEF   63   0     1    4
  11.11     0.00(  0.00) 0.050(INDEF)   0.00(INDEF) INDEF   62   6     1    4
  11.99     0.00(  0.00) 0.050(INDEF)   0.00(INDEF) INDEF   72   1     1    4
  12.95     0.00(  0.00) 0.050(INDEF)   0.00(INDEF) INDEF   79   0     1    4
  13.99     0.00(  0.00) 0.050(INDEF)   0.00(INDEF) INDEF   86   0     1    4
  15.11     0.00(  0.00) 0.050(INDEF)   0.00(INDEF) INDEF   92   0     1    4
  16.32     0.00(  0.00) 0.050(INDEF)   0.00(INDEF) INDEF   96   4     1    4
  17.62     0.00(  0.00) 0.050(INDEF)   0.00(INDEF) INDEF  106   2     1    4
  19.03     0.00(  0.00) 0.050(INDEF)   0.00(INDEF) INDEF  108   8     1    4
  20.56     0.00(  0.00) 0.050(INDEF)   0.00(INDEF) INDEF   80  46     1    4
  22.20     0.00(  0.00) 0.050(0.001)   0.00( 0.32) INDEF   59  77     1    1
   5.56     0.00(  0.00) 0.050(0.004)   0.00( 2.64) 0.142   34   0     1    0
   5.14     0.00(  0.00) 0.050(0.005)   0.00( 2.98) 0.169   32   0     1    0
   4.76     0.00(  0.00) 0.050(0.007)   0.00( 3.97) 0.185   29   0     1    0
   4.41     0.00(  0.00) 0.050(0.006)   0.00( 3.68) 0.225   27   0     1    0
   4.08     0.00(  0.00) 0.050(0.007)   0.00( 4.16) 0.257   25   0     1    0
   3.78     0.00(  0.00) 0.050(0.009)   0.00( 5.76) 0.284   24   0     1    0
   3.50     0.01(  0.00) 0.050(0.008)   0.00( 4.89) 0.352   22   0     1    0
   3.24     0.01(  0.00) 0.050(0.010)   0.00( 6.20) 0.423   20   0     1    0
   3.00     0.01(  0.00) 0.050(0.010)   0.00( 6.24) 0.444   19   0     1    0
   2.78     0.01(  0.01) 0.050(0.013)   0.00( 7.53) 0.497   17   0     1    0
   2.57     0.01(  0.01) 0.050(0.016)   0.00( 9.34) 0.612   16   0     1    0
   2.38     0.02(  0.01) 0.050(0.013)   0.00( 8.11) 0.836   15   0     1    0
   2.21     0.02(  0.01) 0.050(0.012)   0.00( 7.08) 0.879   14   0     1    0
   2.04     0.02(  0.01) 0.050(0.015)   0.00( 9.36) 0.916   13   0     1    0
   1.89     0.03(  0.01) 0.050(0.015)   0.00( 9.15) 0.988   13   0     1    0
   1.75     0.03(  0.01) 0.050(0.016)   0.00( 9.64) 1.112   13   0     1    0
   1.62     0.03(  0.01) 0.050(0.016)   0.00( 9.62) 1.259   13   0     1    0
   1.50     0.03(  0.01) 0.050(0.017)   0.00(10.50) 1.394   13   0     1    0
   1.39     0.04(  0.01) 0.050(0.020)   0.00(12.17) 1.506   13   0     1    0
   1.29     0.04(  0.01) 0.050(0.025)   0.00(15.26) 1.658   13   0     1    0
   1.19     0.04(  0.01) 0.050(0.033)   0.00(20.01) 1.815   13   0     1    0
   1.10     0.04(  0.01) 0.050(0.063)   0.00(38.33) 2.292   13   0     1    0
   1.02     0.04(  0.01) 0.050(0.087)   0.00(52.66) 2.571   13   0     1    0
   0.95     0.04(  0.01) 0.050(0.103)   0.00(62.30) 3.123   13   0     1    0
   0.88     0.04(  0.01) 0.050(0.105)   0.00(63.88) 3.272   13   0     1    0
   0.81     0.05(  0.01) 0.050(0.103)   0.00(62.79) 3.274   13   0     1    0
   0.75     0.05(  0.01) 0.050(0.102)   0.00(61.88) 3.275   13   0     1    0
   0.70     0.05(  0.01) 0.050(0.101)   0.00(61.14) 3.277   13   0     1    0
   0.64     0.05(  0.01) 0.050(0.100)   0.00(60.54) 3.279   13   0     1    0
   0.60     0.05(  0.01) 0.050(0.099)   0.00(60.06) 3.280   13   0     1    0
   0.55     0.05(  0.01) 0.050(0.098)   0.00(59.66) 3.282   13   0     1    0
   0.51     0.05(  0.01) 0.050(0.098)   0.00(59.35) 3.284   13   0     1    0
   0.00     0.06
   0.07  CPU seconds.
   0.00  minutes elapsed.

















In [217]:



    


# Dump the useful information into a ASCII table
iraf.tables()
iraf.ttools()

if os.path.exists(psfTab):
    os.remove(psfTab)
if os.path.exists(psfCdf):
    os.remove(psfCdf)

iraf.tdump.columns=''
iraf.tdump(psfBin, datafil=psfTab, cdfile=psfCdf)



    


















    






IMAGE    t red_21572_Ipsf.fits

















In [218]:



    


ellipsePsf = readEllipseOut(psfTab, zp=zpPhoto, pix=pixelScale, exptime=expTime)

nIsoPsf = len(ellipsePsf)
print "%d elliptical isophotes have been extracted from PSF image" % nIsoPsf

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



    


















    






51 elliptical isophotes have been extracted from PSF image
The maximum radius is   22.20 pixels










    






/Users/songhuang/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:60: RuntimeWarning: invalid value encountered in log10

















In [236]:



    


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('${\mu}_{i}$ (mag/arcsec$^2$)', fontsize=23)

ax.plot(ellipseOut3['rsma'], ellipseOut3['sbp_uerr'], '-', 
        color='k', linewidth=2.0, alpha=0.6)
ax.plot(ellipseOut3['rsma'], ellipseOut3['sbp_lerr'], '-', 
        color='k', linewidth=2.0, alpha=0.6)
ax.plot(ellipseOut3['rsma'], ellipseOut3['sbp'], '-', 
        color='r', linewidth=2.0, label='${\mu}_{i}$')

sbpBuffer = 0.5  
minSbp, maxSbp = np.nanmin(ellipseOut3['sbp_uerr']), np.nanmax(ellipseOut3['sbp_lerr'])
ax.set_ylim((maxSbp+sbpBuffer), (minSbp-sbpBuffer))

# Overplot the profile for PSF 

# Find the offset at the center 
sbpOffsetPsf = ellipseOut3['sbp'][0] - ellipsePsf['sbp'][0]

ax.plot(ellipsePsf['rsma'], (ellipsePsf['sbp'] + sbpOffsetPsf), '--', 
        color='g', linewidth=2.0, label='PSF')

ax.legend(loc=[0.04,0.04], fontsize=18)



    


















    
Out[236]:








<matplotlib.legend.Legend at 0x124578750>










    
























In [483]:



    


len(ellipseOut3['rsma'])



    


















    
Out[483]:








90

















In [ ]:

Content source: dr-guangtou/hs_galphot

Similar notebooks:

notebook.community | gallery | about