"""A set of routines to assist reduction of Spitzer MIPS (24 micron) data.
To prepare a run (set AORs, frame numbers, etc.) edit :func:`initobs`
To generate empirical flat-fields, see :func:`mipsflat`
To extract photometry (aperture or PSF-fitting), see :func:`mipsphot`
To better array photometry, remove first-frames, etc., see :func:`binphot`
:REQUIREMENTS:
:doc:`pylab`, :doc:`numpy`, :doc:`os`, :doc:`nsdata`,
:doc:`pyfits`,:doc:`phot`,:doc:`tools`, and maybe some
others
:HISTORY:
2009-08-24 15:55 IJC: Created
2011-06-15 09:37 IJMC: Better Sphinx documentation.
"""
#try:
# import psyco
# psyco.full()
#except ImportError:
# print 'Psyco not installed, the program will just run slower'
import pdb
[docs]class baseObject:
"""Empty object container.
"""
# 2010-01-24 15:13 IJC: Added to spitzer.py (from my ir.py)
def __init__(self):
return
def __class__(self):
return 'baseObject'
class obs:
def __init__(self, path=None, files=None, aor=None, expid=None, dce=None, \
target=None, tbins=None):
self.aor = aor
self.expid = expid
self.files = files
self.path = path
self.dce = dce
self.target = target
self.tbins = tbins
[docs]class genericObject:
"""Object generated with whatever attributes are passed to it.
Tested for object types of:
str
list
dict
array <--- _SMALL_ arrays of small values only... limited utility!"""
def __init__(self, **kw):
from numpy import ndarray, array
#print kw
#print kw.keys()
for k in kw.keys():
#print k, kw[k]
if kw[k].__class__==str:
exec('self.'+k+'="'+str(kw[k]))+'"'
elif isinstance(kw[k], ndarray):
temp = str(kw[k]).replace(' ',', ')
string = 'array('+temp+')'
exec('self.'+k+'='+string)
else:
#print 'self.'+k+'='+str(kw[k])
exec('self.'+k+'='+str(kw[k]))
[docs]def mipsfile(*args):
"""Generate a MIPS data file name according to the standard in the
MIPS Data Handbook v3.2:
SPITZER_M+chnlnum_aorkey_expid_dcenum_version_type.suffix
:INPUT: separate values, or a list, or lists of lists, representing the quantities:
chnlnum, aor, expid, dce, ver, type, suffix
ALTERNATIVELY, input can be a (str) filename
:OUTPUT:
either:
(1) list of strings or
(2) string corresponding to all combinations of the input paramters.
ALTERNATIVELY, output is a tuple of type (ch,aor,exp,dce,ver,type,suff)
:EXAMPLES:
::
inp = (1, 26598144, 7, range(3), 1, ['bunc', 'bcd'], 'fits')
spitzer.mipsfile(inp)
# OR
spitzer.mipsfile(1, 26598144, 7, range(3), 1, ['bunc', 'bcd'], 'fits')
# OR
spitzer.mipsfile('SPITZER_M1_26598144_0007_0000_1_bunc.fits')
"""
# 2009-08-24 16:14 IJC: Created.
# 2009-08-27 14:54 IJC: Updated to check for numpy array
# 2009-09-02 16:02 IJC: Added str filename option
# 2009-12-14 11:50 IJC: Chop off root path.
# 2009-12-16 15:35 IJC: Accept either python str or numpy string
from nsdata import intstr
from numpy import ndarray, string_
import os
while len(args)==1:
args = args[0]
if len(args)<7:
print "Need seven arguments passed to MIPSFILE!"
return '-1'
# If a string was input, extract desired parameters from the filename:
if args.__class__==str or args.__class__==string_:
args = os.path.split(args)[1]
ch = int(args[9])
aor = int(args[11:19])
expid = int(args[20:24])
dce = int(args[25:29])
ver = int(args[30])
type = args[32:args.find('.')]
suff = args[args.find('.')+1::]
return (ch, aor, expid, dce, ver, type, suff)
# If not a string, must be a sequence of values: create filename
nargs = 7
args = list(args)
for ii in range(nargs):
thisclass = args[ii].__class__
if isinstance(args[ii],ndarray):
args[ii] = list(args[ii])
elif thisclass<>list and thisclass<>tuple:
args[ii] = [str(args[ii])]
ch, aor, expid, dce, ver, type, suff = args
#print args
ret = []
for c in ch:
for a in aor:
for e in expid:
for d in dce:
for v in ver:
for t in type:
for s in suff:
#print c,a,e,d,v,t,s
if s[0]=='.': s=s[1::]
part1 = "SPITZER_M"+c+"_"+str(a)+"_"
part2 = intstr(e,4)+'_'+intstr(d,4)+'_'
part3 = v+"_"+t+"."+s
ret.append(part1+part2+part3)
if len(ret)==1:
ret = ret[0]
return ret
[docs]def initobs(obsname='upsand1', dce='all', expid='all', verbose=False):
"""Initialize variables for SPITZER data analysis.
:OPTIONAL INPUT:
obsname: (str) Current options are 'upsand1', 'upsand2',
'hd209458_0', 'hd209458', 'hd189733_1', 'hd189733_2'
dce: (list) Which DCEs to load. Defaults to 0-35 inclusive.
:OUTPUT:
a DICT containing the following values, in order:
datalist -- list of BCD data file names to analyse
_bcd -- BCD (processed data) directory
"""
# 2009-08-27 12:51 IJC: Created; much simpler than its NIRSPEC counterpart!
# 2009-08-31 17:21 IJC: Checked that this produces a full and complete list.
# 2009-09-21 14:22 IJC: Added Knutson/Charbonneau HD_189733 phase curves.
# 2009-09-23 15:52 IJC: Added Knutson/Charbonneau HD_209458 phase curves.
# 2009-12-14 10:40 IJC: Updated doc string with current valid obsnames.
from numpy import tile, sort
tbins = None
if obsname=='upsand2':
planet = 'upsilon And b'
_path = '/Users/ianc/proj/spitzer/data/upsand_2/bcd/'
aors = [26574592, 26574848, 26575104, 26575360, 26575616, 26575872,
26576128, 26576384, 26576640, 26576896, 26577152, 26577408,
26577664, 26577920, 26578176, 26578432, 26578688, 26578944,
26579200, 26579456, 26579712, 26579968, 26580224, 26580480,
26580736, 26580992, 26581248, 26581504, 26581760, 26582016,
26582272, 26582528, 26582784, 26583040, 26583296, 26583552,
26583808, 26584064, 26584320, 26584576, 26584832, 26585088,
26585344, 26585600, 26585856, 26586112, 26586368, 26586624,
26586880, 26587136, 26587392, 26587648, 26587904, 26588160,
26588416, 26588672, 26588928, 26589184, 26589440, 26589696,
26589952, 26590208, 26590464, 26590720, 26590976, 26591232,
26591488, 26591744, 26592000, 26592256, 26592512, 26592768,
26593024, 26593280, 26593536, 26593792, 26594048, 26594304,
26594560, 26594816, 26595072, 26595328, 26595584, 26595840,
26596096, 26596352, 26596608, 26596864, 26597120, 26597376,
26597632, 26597888, 26598144, 26598400, 26598656, 26598912,
26599168, 26599424]
long = range(8)
short = range(4)
if expid<>'all':
long = list(set(long).intersection(set(expid)))
short = list(set(short).intersection(set(expid)))
lls = [long,long,short]
six = [list(el) for el in list(tile(long,(6,1)))]
big = [list(el) for el in list(tile(long,(71,1)))]
expids = lls+lls+lls+lls+big + lls+lls+lls+six
versions = [1]*len(expids)
dce_normal = range(36)
if dce<>'all':
dce_normal = list( set(dce_normal).intersection(set(dce)) )
dce_normal = [list(sort(dce_normal))]
dce = dce_normal*len(aors)
tbins = [0,.7,1.2,1.5,1.8,2.4,2.7,2.79,2.912,3.034,3.2,3.373,3.51,3.67,3.809,4.2,4.7195,5.2, 6,6.4,7]
elif obsname=='upsand2cal':
planet = 'upsilon And b'
_path = '/Users/ianc/proj/spitzer/data/upsand_2/bcd/'
aors = [26598144, 26598400, 26598656, 26598912, 26599168, 26599424]
naor = len(aors)
expid_normal = range(8)
if expid<>'all':
expid_normal = list(set(expid_normal).intersection(set(expid)))
expids = [expid_normal]*naor
versions = [1]*naor
dce_normal = range(36)
if dce<>'all':
dce_normal = list( set(dce_normal).intersection(set(dce)) )
dce_normal = [list(sort(dce_normal))]
dce = dce_normal*naor
elif obsname=='upsand1':
planet = 'upsilon And b'
_path = '/Users/ianc/proj/spitzer/data/upsand_1/bcd/'
aors = [14037760, 14038016, 14038272, 14041344, 14041600]
short = range(4)
if expid<>'all':
short = list(set(short).intersection(set(expid)))
expids = [short] * 5
versions = [3]*len(expids)
dce_normal = range(43)
if dce<>'all':
dce_normal = list( set(dce_normal).intersection(set(dce)) )
dce_normal = [list(sort(dce_normal))]
dce = dce_normal*len(aors)
tbins = [0., 1.0, 2.0, 2.5, 3.0, 5.0]
elif obsname=='tauboo':
planet = 'tau Boo b'
_path = '/Users/ianc/proj/spitzer/data/tauboo/bcd/'
aors = [14036992, 14037248, 14037504, 14040832, 14041088]
short = range(4)
if expid<>'all':
short = list(set(short).intersection(set(expid)))
expids = [short] * 5
versions = [5]*len(expids)
dce_normal = range(43)
if dce<>'all':
dce_normal = list( set(dce_normal).intersection(set(dce)) )
dce_normal = [list(sort(dce_normal))]
dce = dce_normal*len(aors)
tbins = [0., 1.0, 2.0, 2.5, 3.0, 5.0]
elif obsname=='51peg':
planet = '51 Peg b'
_path = '/Users/ianc/proj/spitzer/data/51peg/bcd/'
aors = [14038528, 14038784, 14039040, 14041856, 14042112]
short = range(4)
if expid<>'all':
short = list(set(short).intersection(set(expid)))
expids = [short] * 5
versions = [8]*len(expids)
dce_normal = range(43)
if dce<>'all':
dce_normal = list( set(dce_normal).intersection(set(dce)) )
dce_normal = [list(sort(dce_normal))]
dce = dce_normal*len(aors)
tbins = [0., 1.0, 2.0, 2.5, 3.0, 5.0]
elif obsname=='hd75289':
planet = 'HD 75289 b'
_path = '/Users/ianc/proj/spitzer/data/hd75289/bcd/'
aors = [14040064, 14040320, 14040576]
short = range(4)
if expid<>'all':
short = list(set(short).intersection(set(expid)))
expids = [short] * 3
versions = [8]*len(expids)
dce_normal = range(43)
if dce<>'all':
dce_normal = list( set(dce_normal).intersection(set(dce)) )
dce_normal = [list(sort(dce_normal))]
dce = dce_normal*len(aors)
tbins = [0., 1.0, 2.0, 2.5, 3.0, 5.0]
elif obsname=='hd179949':
planet = 'HD 179949 b'
_path = '/Users/ianc/proj/spitzer/data/hd179949/bcd/'
aors = [14039296, 14039552, 14039808]
short = range(4)
if expid<>'all':
short = list(set(short).intersection(set(expid)))
expids = [short] * 3
versions = [8]*len(expids)
dce_normal = range(43)
if dce<>'all':
dce_normal = list( set(dce_normal).intersection(set(dce)) )
dce_normal = [list(sort(dce_normal))]
dce = dce_normal*len(aors)
tbins = [0., 1.0, 2.0, 2.5, 3.0, 5.0]
elif obsname=='hd189733_1':
planet = 'HD 189733 b'
_path = '/Users/ianc/proj/spitzer/data/hd189733_1/bcd/'
aors = [16343040, 16343296]
long = range(18)
if expid<>'all':
long = list(set(long).intersection(set(expid)))
long = [long]
expids = 2*long
versions = [4, 4]
dce_normal = range(43)
if dce<>'all':
dce_normal = list( set(dce_normal).intersection(set(dce)) )
dce_normal = [list(sort(dce_normal))]
dce = dce_normal*2
elif obsname=='hd209458_0':
planet = 'HD 209458 b'
_path = '/Users/ianc/proj/spitzer/data/hd209458_0/bcd/'
aors = [10902272, 10902528, 10902784, 10903040, 10903296, 10903552, \
10903808, 10904064, 10904320, 10904576, 10904832, 10905088]
naors = len(aors)
# all 8x36
long = range(8)
if expid<>'all':
long = list(set(long).intersection(set(expid)))
long = [long]
expids = long*naors
versions = [5,5,5,7,7,7,7,7,7,7,7,7]
dce_normal = range(36)
if dce<>'all':
dce_normal = list( set(dce_normal).intersection(set(dce)) )
dce_normal = [list(sort(dce_normal))]
dce = dce_normal*naors
elif obsname=='hd209458_1':
planet = 'HD 209458 b'
_path = '/Users/ianc/proj/spitzer/data/hd209458_1/bcd/'
aors = [14909952, 14910208, 14910464, 14910720, 14910976, \
14911232, 14911488, 14911744, 14912000]
naors = len(aors)
# all 8x36
long = range(8)
if expid<>'all':
long = list(set(long).intersection(set(expid)))
long = [long]
expids = long*naors
versions = [9,9,9, 8,8,8,8,8,8]
dce_normal = range(36)
if dce<>'all':
dce_normal = list( set(dce_normal).intersection(set(dce)) )
dce_normal = [list(sort(dce_normal))]
dce = dce_normal*naors
elif obsname=='hd189733_2':
planet = 'HD 189733 b'
_path = '/Users/ianc/proj/spitzer/data/hd189733_2/bcd/'
aors = [22863872, 22864128, 22864384, \
22864640, 22864896, 22865152, 22865408, 22865664, \
22865920, 22866176, 22866432, 22866688, 24412928]
long = range(24)
medium = range(19)
short = range(12)
shorter = range(6)
if expid<>'all':
long = list(set(long).intersection(set(expid)))
medium = list(set(medium).intersection(set(expid)))
short = list(set(short).intersection(set(expid)))
shorter = list(set(shorter).intersection(set(expid)))
long = [long]
medium = [medium]
short = [short]
shorter = [shorter]
#expids = 7*long+ short+ 11*long+ short+ long*4 + \
# medium+ long*4+ short
expids = 7*long+short+4*long+shorter
versions = [1]*13 #+[3]*19
dce_short = range(8)
dce_normal = range(36)
dce_long = range(50)
if dce<>'all':
dce_short = list( set(dce_short).intersection(set(dce)) )
dce_normal = list( set(dce_normal).intersection(set(dce)) )
dce_long = list( set(dce_long).intersection(set(dce)) )
dce_short = [list(sort(dce_short))]
dce_normal = [list(sort(dce_normal))]
dce_long = [list(sort(dce_long))]
dce = dce_normal*7+dce_long+dce_normal*4+dce_short
elif obsname=='hd209458':
planet = 'HD 209458 b'
_path = '/Users/ianc/proj/spitzer/data/hd209458/bcd/'
aors =[22866944, 22867200, 22867456, 22867712, 22867968, 22868224, \
22868480, 22868736, 22868992, 22869248, 22869504, \
22869760, 22870016, 22870272, 22870528, 22870784, \
22871040, 22871296, 22873344] #, 24412928]
long = range(24)
medium = range(19)
short = range(12)
if expid<>'all':
long = list(set(long).intersection(set(expid)))
medium = list(set(medium).intersection(set(expid)))
short = list(set(short).intersection(set(expid)))
long = [long]
medium = [medium]
short = [short]
#expids = 7*long+ short+ 11*long+ short+ long*4 + \
# medium+ long*4+ short
expids = 7*long + short + 9*long + short + long
versions = [3]*19 #+[3]*19
dce_short = range(8)
dce_medium = range(29)
dce_normal = range(36)
dce_long = range(50)
if dce<>'all':
dce_short = list( set(dce_short).intersection(set(dce)) )
dce_medium = list( set(dce_medium).intersection(set(dce)) )
dce_normal = list( set(dce_normal).intersection(set(dce)) )
dce_long = list( set(dce_long).intersection(set(dce)) )
dce_short = [list(sort(dce_short))]
dce_medium = [list(sort(dce_medium))]
dce_normal = [list(sort(dce_normal))]
dce_long = [list(sort(dce_long))]
dce = dce_normal*7 + dce_long + dce_normal*4 + dce_medium + \
dce_normal*4 + dce_long + dce_normal
filelist = []
if verbose:
print "aors>>"+str(aors)
print "expids>>"+str(expids)
print "versions>>"+str(versions)
for aor, exp, dce0, ver in zip(aors, expids, dce, versions):
if verbose: print "(aor, exp, dce0, ver)>>" + str((aor, exp, dce0, ver))
filelist.append(mipsfile(1,aor,exp,dce0,ver,'bcd','fits'))
return obs(path=_path, files=filelist, aor=aors, expid=expids, dce=dce, \
target=planet, tbins=tbins)
[docs]def mipspos(dcenum, expid, retfnum=False):
"""Return the approximate x,y location of a star on the MIPS detector.
:INPUT:
expid -- (int) exposure ID number.
dcenum -- (int or list of ints) DCE, 0-35 inclusive.
:OPTIONAL INPUT:
retfnum -- (bool). Whether to return the frame position number
(from 0-13), as listed in Table 8.7 of the Spitzer
Observer's Manual.
:OUTPUT:
(xpos, ypos) -- tuple of floats, if retfnum==False
( (xpos, ypos), fnum) -- tuple of (tuple, floats)
:EXAMPLE:
::
mipspos(0,0)
mipspos(35,1, retfnum=True)
[mipspos(range(36), expid) for expid in [0,1]]
:NOTE:
Based on the prescription of even DCE on the LHS, odd DCE on
the RHS, and per Table 8.7 in the Spitzer Observer's Manual.
"""
# 2009-09-23 13:58 IJC: Made it work with dce>36
# 2009-09-23 16:22 IJC: Corrected confusion of EXPID & DCE
# 2009-11-16 12:50 IJC: Switched EXPID & DCE again...
from numpy import tile, array, concatenate, arange,max,ceil
if (dcenum % 2)==0:
xoff = 39.0
else:
xoff = 89.5
numberOfReps = ceil(max(expid)/7.0)
yoff = [ 80.7, 33.5, 85.3, 38.1, 89.9, 42.7, 94.5]
yoff = list(tile(yoff,numberOfReps)) + [yoff[0]]
expid = array(expid).astype(int)
yret = array(yoff)[expid]
if yret.shape==():
ret = (xoff, yret)
else:
ret = (tile(xoff, len(yret)), yret)
if retfnum:
frameNumberOffset = 7*(dcenum%2)
frameOrder = concatenate((tile(arange(7),numberOfReps), [0]))
frameNumber = frameNumberOffset+frameOrder[expid]
ret = (ret, frameNumber)
return ret
[docs]def pxsize(fn, key1='PXSCAL1', key2='PXSCAL2'):
"""Compute size of pixels (in steradians) for a given FITS image.
:INPUT:
fn -- str. filename
:OPTIONAL INPUT:
key1 -- str. FITS Header Key of first pixel scale value, in arcsec/pix
key2 -- str. FITS Header Key of second pixel scale value, in arcsec/pix
:OUTPUT:
a -- float. area of a pixel, in steradians.
"""
# 2009-09-10 15:16 IJC: Created
# 2014-09-03 14:06 IJMC: Added PyFITS/astropy check.
try:
from astropy.io.fits import getval
except:
from pyfits import getval
scale1 = getval(fn, key1)
scale2 = getval(fn, key2)
return abs(scale1*scale2*2.4129606682083946e-11)
[docs]def mipsaperphot(fn, dap=(11,13,41), mask=None, fix=False, verbose=False, center=None):
"""Perform aperture photometry on a MIPS file or list of MIPS files.
:INPUT:
fn -- (str) filename, or (list) list of filenames
:OPTIONAL INPUT:
dap -- (3-sequence) pixel diameters: target aperture, inner sky, outer sky
mask -- (array) mask to pass to phot.aperphot. Otherwise, generate one from dap.
:NOTES:
** if no center location is passed, target positions will be taken
from Table 8.7 of the S.O.M. based on the MIPS filename.
** this works if you pass in an array or a filename, but
passing a filename enables the function to do a more
complete job (grabbing object name,RA/DEC, HJD fields, etc.)
"""
# 2009-09-14 12:59 IJC: Created
# 2009-09-27 14:49 IJC: accepts MASK input to pass to phot.aperphot
# 2009-11-16 12:52 IJC: Swapped exp & DCE in mipspos function call
import os, phot, pyfits
from numpy import ndarray
def getpositions(fn, center, mask):
"""Get photometry positions based on parameters passed"""
if isinstance(fn,str):
ch, aor, dce, exp, ver, type, suff = mipsfile(os.path.split(fn)[1])
if center==None and mask==None:
(x0,y0), fnum = mipspos(dce, exp, retfnum=True)
elif center<>None:
x0,y0 = center
fnum = None
else:
x0=None
y0=None
fnum=None
elif isinstance(fn,ndarray):
fnum = None
if center==None:
x0=None
y0=None
else:
x0,y0 = center
return x0,y0,fnum
if isinstance(fn, list):
obs = [mipsaperphot(file, dap=dap, mask=mask, fix=fix, verbose=verbose, center=center) for file in fn]
elif isinstance(fn, str) or isinstance(fn,ndarray):
x0,y0,fnum = getpositions(fn, center, mask)
if verbose: print '(x0,y0)>>' + str((x0,y0))
if isinstance(fn,str):
if verbose: print 'fn>> %s' % fn
pixelsize = pxsize(fn)*1e6
photunit = "Jy"
elif isinstance(fn,ndarray):
if center==None:
print "must enter a center location for the aperture!"
return -1
pixelsize = 1
photunit = "MJy pix / sr"
obs = phot.aperphot(fn, pos=(x0,y0), mask=mask, dap=dap, timekey='HMJD_OBS', verbose=verbose)
obs.phot = obs.phot*pixelsize
obs.photunit=photunit
obs.timestr='heliocentric modified julian date'
obs.ditherposition=fnum
obs.ditherpositionstr="add 1 (one) for index in Table 8.7 in the Spitzer Observer's Manual."
if isinstance(fn, str):
obs.ra = pyfits.getval(fn, 'RA_REF')
obs.rastr = '[deg] Commanded RA (J2000) of ref. position '
obs.dec = pyfits.getval(fn, 'DEC_REF')
obs.decstr = '[deg] Commanded DEC (J2000) of ref. position '
return obs
[docs]def mipsphot(dataset, dithernum=0, dce='all', expid='all', aper=[13,23,41], \
sigma=5, verbose=False, plotalot=False, \
psffile='/Users/ianc/proj/spitzer/psf/psfx100_6200K_sm.fits', \
pmask='/Users/ianc/proj/spitzer/psf/mips24_pmask.fits', scale=100.,\
fileind='all', dflatfield=None, sflatfield=None, offset=[0,0]):
"""
Analyze a Spitzer/MIPS 24um photometric dataset.
:INPUTS:
dataset -- (str) 'upsand1', 'upsand2', 'upsand2cal', 'hd209458_0', 'hd209458', 'hd189733_1', 'hd189733_2'
dce & expid -- if 'all', use all dce & expid values valid for the specified dithernum.
aper -- (tuple) [targ, sky_inner, sky_outer]
aperture diameters for photometry, in pixels. For
PSF-fitting, this must be _odd_.
sigma -- sigma-clipping value when computing means
psffile -- model PSF to use for PSF-fitting photometry (from TinyTim).
pmask -- filename for MIPS bad pixel mask
scale -- scaling of model PSF relative to physical MIPS detector.
fileind -- either set to 'all', or a sequence of which files to extract
dflatfield -- filename of a FITS flat-field by which to divide all frames before photometry.
sflatfield -- filename of a FITS flat-field by which to subtract all frames before photometry.
dithernum -- should be any value from 0 to 13, inclusive.
offset -- pixel offset from nominal MIPS 14-dither position
:EXAMPLE:
::
import spitzer
outputs = []
for ii in range(14):
outputs.append(spitzer.mipsphot(dataset='upsand1', dithernum=ii))
# OR use different model PSFs for each position:
psfs = ['/Users/ianc/proj/spitzer/psf/psfx100_6200K_%s_sm.fits' %s \\
for s in ['39.0x87.6','39.0x38.1','89.5x87.6','89.5x38.1']]
dnum = [[0,2,4,6],[1,3,5],[7,9,11,13],[8,10,12]]
outputs = []
for ii in range(4):
for jj in dnum[ii]:
outputs.append(spitzer.mipsphot(dataset='upsand1',dithernum=jj,psffile=psfs[ii]))
"""
# 2009-12-14 10:03 IJC: Created function and added to spitzer.py
# 2010-01-19 11:04 IJC: Now save 'chisq' and FITS headers
# 2010-01-26 10:31 IJC: Headers were saved incorrectly: now updated.
# 2010-02-22 13:49 IJC: Redirected PSF file.
# 2010-06-02 17:19 IJC: Added bad pixel masking per J. Harrington's advice.
# 2010-06-11 11:09 IJC: fixed quadrant-mask; added flat-field option.
# 2010-07-29 09:22 IJC: Now report PSF-fitting photometry in Jy
# 2011-05-11 09:33 IJMC: Replace fmin_powell with fmin; Only do
# aperture photometry on selected files
# 2011-05-22 13:39 IJMC: Added options for either division or
# subtraction by flat fields.
# 2012-07-18 15:35 IJMC: Switched 'subreg' calls to 'subreg2'
# 2012-09-27 19:47 IJMC: Now save a few of the more important options.
import phot, pyfits
import analysis as an
from nsdata import collapse_objlist
from tools import flatten, keylist
from time import time
from numpy import arange, meshgrid, array, isnan, isfinite, zeros, mean, nan, unique, Inf,median, abs, random
from scipy.optimize import fmin_powell, fmin
import os
import pdb
if dce=='all':
dce = arange(0,60,7)+ (dithernum % 7)
if expid=='all':
expid = arange(0,30,2) + (dithernum >= 7)
keys = ['CSM_PRED', 'CSM_SKY', 'RMS_JIT', 'RMS_JITY', 'RMS_JITZ', 'ZODY_EST', 'ISM_EST', \
'CSM_POS', 'CSM_POSF', 'CSM_POSC', 'ANCSMGN', 'CSM_OFFS', 'AD24TMPA', \
'AD24TMPB', 'ACSMMTMP', 'ACEBOXTM', 'APROFF1V', 'APROFF2V', 'APROFF3V', \
'APROFF4V', 'AD24ANLI', 'PXSCAL1', 'PXSCAL2','DRIBKGND']
dap_targ, dap_skyinner, dap_skyouter = aper
sz = (dap_targ,dap_targ)
dap = (dap_targ, dap_skyinner, dap_skyouter)
fileind = array(fileind, copy=True)
setup = initobs(dataset, dce=dce, expid=expid)
planet = an.getobj(setup.target)
filelist = [setup.path + el for el in flatten(setup.files)]
unc_filelist = [setup.path + el.replace('_bcd','_bunc') for el in flatten(setup.files)]
bbmask_filelist = [setup.path + el.replace('_bcd','_bbmsk') for el in flatten(setup.files)]
ch, aor, dc, ex, ver, type, suff = mipsfile(os.path.split(filelist[0])[1])
nominal_center = mipspos(dc, ex)
new_center = [nc + off for nc, off in zip(nominal_center, offset)]
if verbose:
print "Using as center: ", new_center
obs1 = mipsaperphot(filelist[0], center=new_center)
cen = obs1.position[1],obs1.position[0]
if verbose:
print "Using as center: ", cen
if dflatfield is None and sflatfield is None:
flat = 1.+zeros((128,128),float)
divFlat = False
subFlat = False
elif dflatfield is not None:
flat = pyfits.getdata(dflatfield)
divFlat = True
subFlat = False
elif sflatfield is not None:
flat = pyfits.getdata(sflatfield)
divFlat = False
subFlat = True
tbins = setup.tbins
x = range(128)
xx,yy = meshgrid(x,x)
if verbose:
print "filelist[0]>>",filelist[0]
print "%i files in list" % len(filelist)
print "dce[0]>>", dce[0]
print "expid[0]>>", expid[0]
# Load the MIPS bad-pixel mask:
if pmask is not None:
try:
badpixelmask = pyfits.getdata(pmask)
except:
print "pmask file %s could not be found!" % str(pmask)
badpixelmask = zeros((128.,128.),float)
else:
print "no pmask file entered!"
badpixelmask = zeros((128.,128.),float)
pixelmask = badpixelmask==0 # all zero-values pixels are now "True"
if dithernum in set((0,2,4,6)):
quadmask = pixelmask*phot.makemask(xx,yy,shape='quad',params=[1,2,3])
elif dithernum in set((1,3,5)):
quadmask = pixelmask*phot.makemask(xx,yy,shape='quad',params=[0,2,3])
elif dithernum in set((7,9,11,13)):
quadmask = pixelmask*phot.makemask(xx,yy,shape='quad',params=[0,1,3])
elif dithernum in set((8,10,12)):
quadmask = pixelmask*phot.makemask(xx,yy,shape='quad',params=[0,1,2])
else:
quadmask = []
targmask = pixelmask*phot.makemask(xx,yy,shape='circle',params=list(cen)+[dap_targ])
skymask_in = pixelmask*phot.makemask(xx,yy,shape='circle',params=list(cen)+[dap_skyinner])
skymask_out = pixelmask*phot.makemask(xx,yy,shape='circle',params=list(cen)+[dap_skyouter])
apermask = targmask+2*(skymask_out-skymask_in)
obs = []
quadbg_list = []
equadbg_list = []
tic=time()
nfiles = len(filelist)
if fileind=='all':
fileind = arange(nfiles)
nfiles2 = len(fileind)
# for jj in range(nfiles2):
# ii = fileind[jj]
for jj in range(nfiles2): # Only do aperture photometry on selected files
ii = fileind[jj]
file = filelist[ii]
obs.append(mipsaperphot(file, center=new_center)) #, mask=apermask, center=cen))
if divFlat:
frame = pyfits.getdata(file)/flat;
elif subFlat:
frame = pyfits.getdata(file) - flat
else:
frame = pyfits.getdata(file)
an.fixval(frame,-100)
#print dithernum,frame.shape, flat.shape, quadmask.shape, median(frame),median(flat), median(quadmask)
thisbg, thisebg = phot.estbg(frame,mask=quadmask)
quadbg_list.append(thisbg)
equadbg_list.append(thisebg)
if verbose: print "%i sec to do aperture photometry" % (time()-tic)
objlist = collapse_objlist(obs, 'object')
targetind = array(objlist)==objlist[0] #'hd209458b' #upsilon andromeda'
calind = True - targetind #array(objlist)=='hd9712'
pho = array(collapse_objlist(obs, 'phot'))[targetind]
bg = array(collapse_objlist(obs, 'bg'))[targetind]
quadbg = array(quadbg_list)[targetind]
hjd = an.mjd(array(collapse_objlist(obs, 'time')))[targetind]
exptime = array(collapse_objlist(obs, 'exptime'))[targetind]
fnum = array(collapse_objlist(obs, 'ditherposition'))[targetind]
t = hjd-int(hjd[0])
phase = planet.phase(hjd)
cens = array(collapse_objlist(obs, 'position'))
cens = cens[:,[1,0]]
# Loop through the stack, flagging 3-sigma discrepant pixels
# Apply the flat-field, if any. The stack creation is inside the
# IF block to conserve memory.
if divFlat:
small_flat = phot.subreg2(flatfield, cens, sz)
stack = phot.subreg2(array(filelist)[fileind][targetind], cens, sz)/small_flat
elif subFlat:
small_flat = phot.subreg2(sflatfield, cens, sz)
stack = phot.subreg2(array(filelist)[fileind][targetind], cens, sz) - small_flat
else:
#small_flat = 1.+zeros((1,)+sz,float)
stack = phot.subreg2(array(filelist)[fileind][targetind], cens, sz)#/small_flat
unc_stack = phot.subreg2(array(unc_filelist)[fileind][targetind], cens, sz)
bbmask_stack = phot.subreg2(array(bbmask_filelist)[fileind][targetind], cens, sz)
stack[isnan(stack)] = 0.0
unc_stack[True-isfinite(unc_stack)] = 1e10
cleanmean=zeros(stack.shape[1::], float)
goodind = zeros(stack.shape, bool)
for ii in range(sz[0]):
for jj in range(sz[1]):
temp = an.removeoutliers(stack[:,ii,jj].ravel(),sigma, retind=True)
cleanmean[ii,jj] = mean(temp[0])
goodind[:,ii,jj] = temp[1]
goodind = goodind * (bbmask_stack==0) # additionally flag all BCD-flagged problem pixels
weights_stack = goodind.copy() # / unc_stack**2
if plotalot:
from pylab import figure, subplot, plot, errorbar, legend, xlabel, ylabel, ylim, grid
figure()
subplot(211)
plot([0,1],[2,2], '.-r', [0,1],[2,2], '.-b')
errorbar(t2,pho2/an.meanr(pho2, sigma), epho2/an.meanr(pho2, sigma), ehjd2, fmt='.r', linewidth=2)
errorbar(t2,bg2/an.meanr(bg2, sigma), ebg2/an.meanr(bg2, sigma), ehjd2, fmt='.b', linewidth=2)
leg = legend(['ups and (aper. phot.)', 'detrended annular background'])
xlabel('HJD - %i'% hjd[0], fontsize=20)
ylabel('normalized flux',fontsize=20)
ylim([.99,1.01]); grid()
subplot(212)
plot([0,1],[2,2], '.-m', [0,1],[2,2], '.-g')
errorbar(t2,divCorr2/an.meanr(divCorr2,sigma),edivCorr2/an.meanr(divCorr2,sigma),ehjd2,fmt='.m', linewidth=1.5)
errorbar(t2,qdivCorr2/an.meanr(qdivCorr2,sigma),eqdivCorr2/an.meanr(qdivCorr2,sigma),ehjd2,fmt='.g', linewidth=1.5)
leg=legend(['detrended [ (ups and) / (annular background) ]', 'detrended [ (ups and) / (3-quadrant background) ]'])
xlabel('HJD - %i'% hjd[0], fontsize=20)
ylabel('normalized flux',fontsize=20)
ylim([.99,1.01]); grid()
fullmedian = an.meanr(pho,sigma)
bgmedian = an.meanr(bg,sigma)
phoCopy = pho.copy()
bgCopy = bg.copy()
for pos in unique(fnum):
phoCopy[fnum==pos] = phoCopy[fnum==pos]*fullmedian/an.meanr(phoCopy[fnum==pos],sigma)
bgCopy[fnum==pos] = bgCopy[fnum==pos]*bgmedian/an.meanr(bgCopy[fnum==pos],sigma)
cleanopt = dict(nsigma=3, niter=Inf)
# try some PSF-fitting!
psf = pyfits.getdata(psffile)
psfscale = psf.sum()
headers = zeros(nfiles,object)+nan
chisq = zeros(nfiles,float)+nan
xoff = zeros(nfiles,float)+nan
yoff = zeros(nfiles,float)+nan
locs = zeros((nfiles,2),float)+nan
fitback = zeros(nfiles,float)+nan
fitflux = zeros(nfiles,float)+nan
tic=time()
for ii in range(nfiles2):
if verbose: print "starting file %i of %i." % (ii+1, nfiles2)
#ii = fileind[jj]
frame = stack[ii]
headers[ii] = keylist([filelist[ii]], keys)
headers[ii].name = filelist[ii]
if len(xoff[isfinite(xoff)])==0: # no positions found yet
guess = [0,0]
else:
guess = [xoff[isfinite(xoff)][-1], yoff[isfinite(yoff)][-1]]
# If MOST pixels are bad, throw away the frame:
if goodind[ii].sum()<(0.5*sz[0]*sz[1]):
chisq[ii] = nan
fitback[ii] = nan
fitflux[ii] = nan
xoff[ii] = nan
yoff[ii] = nan
else:
if min(abs(guess))>0:
nonzdelt = 2./min(guess) # minimum discrete step size
else:
nonzdelt = 1. # Guesses are zero, so nonzdelt doesn't matter.
#pdb.set_trace()
decentFit = False
Ntries = 0
while not decentFit and Ntries<10:
out = an.fmin(phot.psffiterr, guess, args=(psf, stack[ii],weights_stack[ii],100,dap_targ),xtol=0.5,ftol=0.25,disp=verbose, zdelt=2, nonzdelt=nonzdelt)
xoff[ii] = round(out[0])
yoff[ii] = round(out[1])
out2 = phot.psffit(psf, frame, w=weights_stack[ii], scale=scale, dframe=dap_targ,\
xoffs=[xoff[ii]],yoffs=[yoff[ii]], verbose=verbose)
chisq[ii] = out2[7]
fitback[ii] = out2[8]
if ii < 10:
decentFit = True # We don't know any better yet...
elif out2[7] > (an.medianr(chisq[0:ii]) + 6):
print "Step %i, initially found a bad chisq (%1.4f), trying again (step %i)" % (ii, chisq[ii], Ntries)
# Try a new guess
guess = random.normal([xoff[isfinite(xoff)][-1], yoff[isfinite(yoff)][-1]], \
[an.stdr(xoff[0:ii]), an.stdr(yoff[0:ii])])
Ntries += 1
else:
decentFit = True
#print headers[ii]
#print dir(headers[ii])
#print 'object', headers[ii].PXSCAL1[0]
this_pixsize_sr = abs(headers[ii].PXSCAL1[0]*headers[ii].PXSCAL2[0])/206265.**2
fitflux[ii] = out2[9]*psfscale*this_pixsize_sr*1e6
if verbose:
print "jj, ii, filename>>",jj,ii,filelist[ii]
print ("%i of %i" % (ii+1,nfiles)),
if verbose: print "%i sec to finish PSF-fitting photometry" % (time()-tic)
options = dict(dithernum=dithernum, aper=aper, sigma=sigma, psffile=psffile)
mydata = dict(hjd=hjd, tbins=tbins, pho=pho, quadbg=quadbg, fitflux=fitflux, fitback=fitback, xoff=xoff, yoff=yoff, obs=obs, fnum=fnum, target=setup.target, dataset=dataset,chisq=chisq,headers=headers,headerkeys=keys, dflatfield=dflatfield, sflatfield=sflatfield, options=options)
return mydata
[docs]def binphot(outputs, dumpframes=None, sigma=3, tbins=None, nbins=None, target=None, binfactor=None, returnstats=False):
"""
Bin down MIPSPHOT time-series photometry.
:INPUT:
outputs -- EITHER:
(list) outputs from MIPSPHOT, one element for
each dither position (a list of (dict) items)
OR:
(str) filename to a pickle containing such a list.
:OPTIONAL INPUT:
dumpframes -- (sequence) TBD -- frames to throw away.
sigma -- (float) value for sigma-clipping
tbins -- (sequence) boundaries of time bins over which to
average the data. If tbins=='file', taken from
datafile; otherwise specified here, or overridden via 'nbins'
nbins -- (int) instead of using 'tbins' flag, automatically
subdivide the binned sequence into 'nbins' equally-spaced
bins. Empty bins will be returned as 'nan'
target -- for datasets with a second calibrator star
(e.g. hd9712), specify that target -- otherwise only
main planetary data is used.
binfactor -- see ANALYSIS.ERRXY
returnstats -- print numbers about # of observations, and # of deleted frames
:EXAMPLE:
::
import spitzer
outputs = []
for ii in range(14):
outputs.append(spitzer.mipsphot(dataset='upsand1', dithernum=ii))
out2 = spitzer.binphot(outputs)
:NOTES:
This calculates several different types of values, based partly
on values extracted from FITS headers and partly on
calculations. If "out2" is the output object, "out2.keys" will
contain a list of all the header keywords used. Each key values
will be binned down by the appropriate factor, and the standard
error on that measurement contained in a related keyword. For
example, out2.AD24TMPA and out2.eAD24TMPA.
Other pertinent values calculated and returned are out2.flux.
out2.fitflux should NOT be used, since it does not properly
normalize for first-frame effects or dither position sensitivity
variations. Similarly for out3.quadbg.
:SEE ALSO: :func:`mipsphot`, :func:`analysis.errxy`
"""
# 2009-12-22 08:49 IJC: Added orbital phase to output
# 2009-12-14 10:46 IJC: Functionalized and added to spitzer.py.
# 2009-12-29 22:35 IJC: Compute correct orbital phase for
# transiting systems. Added 'ephase' to outputs.
# 2010-01-29 12:57 IJC: Updated to use extra FITS header keys.
# Write header keys to output object.
# 2010-02-04 13:05 IJC: Works for 'outputs' inputs of various lengths
# 2010-02-05 14:05 IJC: Output 'dataset' keyword as well.
# 2010-04-29 09:47 IJC: Added printstats for paper-writing details.
# 2011-05-11 11:52 IJC: Updated w/new planet object attributes.
# 2011-10-26 20:41 IJMC: Added slight error-checking for header keywords
# 2012-09-08 21:41 IJMC: Fixed bug in binning xoff and yoff (they
# are now ravel'd)
# 2012-09-27 12:48 IJMC: Now filenames and dce-numbers are saved.
import analysis as an
import copy
import tools as tl
from nsdata import collapse_objlist
from numpy import array, ones, isnan, dot, Inf, polyval, mean, isfinite, linspace, unique
outputs = copy.deepcopy(outputs)
######################################################
def scaleflux(o, key):
""" Helper function to scale fluxes at different positions to
a common mode. Assume the FINAL position is the most
settled:"""
n_o = len(o)
scalefactors = ones(n_o,dtype=float)
for ii in range(n_o-1):
fidflux = o[-1][key]
thisflux = o[ii][key]
ind = isfinite(fidflux) * isfinite(thisflux)
scalefactors[ii] = dot(fidflux[ind], thisflux[ind]) / dot(thisflux[ind], thisflux[ind])
return scalefactors
######################################################
if outputs.__class__<>list:
outputs = [outputs]
cleanopt = dict(nsigma=sigma, niter=Inf)
p = an.getobj(outputs[0]['target'])
while not hasattr(p, 'name'):
newname = raw_input("Enter planet name: ")
p = an.getobj(newname)
#setup = initobs(outputs[0]['dataset'])
init_count = [0]*len(outputs)
targ_count = [0]*len(outputs)
final_count = [0]*len(outputs)
keys = outputs[0].keys()
if tbins=='file':
tbins = outputs[0]['tbins']
if target==None:
target = outputs[0]['obs'][0].object
# Construct DCE arrays, and trim non-primary target observations:
for ii,oo in enumerate(outputs):
init_count[ii] = len(oo['hjd'])
thisexpid = []
thisdce = []
thisfilename = []
for thisobs in oo['obs']:
ch, aor, expid, dce, ver, type, suff = mipsfile(thisobs.filename)
thisexpid.append(expid)
thisdce.append(dce)
thisfilename.append(thisobs.filename)
targs = array(collapse_objlist(oo['obs'], 'object'))
targindex = (targs==target)
thisaperbg = collapse_objlist(oo['obs'],'bg')
oo['dce'] = array(thisdce)[targindex]
oo['expid'] = array(thisexpid)[targindex]
oo['fn'] = array(thisfilename)[targindex]
oo['bga'] = array(thisaperbg)[targindex]
targ_count[ii] = targindex.sum()
# If headers exist, get their names and construct a set of dicts for each dither position
if outputs[0].has_key('headers') and outputs[0].has_key('headerkeys'):
hkeys = outputs[0]['headerkeys']
else:
hkeys = []
hdrs = []
for oo in outputs:
hdrs.append(dict())
for hh,oo in zip(hdrs,outputs):
for key in hkeys:
hh[key] = [] #array(collapse_objlist(list(oo['headers']), key))
for o in oo['headers']:
try:
hh[key].append(getattr(o, key))
except:
hh[key].append(0.)
for key in ['dce', 'pho', 'fnum', 'fitflux', 'bga', 'hjd', 'xoff', 'yoff', 'quadbg', 'fitback','chisq', 'fn', 'expid']:
if oo.has_key(key):
hh[key] = oo[key]
for key in ['dce','pho', 'fnum', 'fitflux', 'bga', 'hjd', 'xoff', 'yoff', 'quadbg', 'fitback','chisq', 'fn', 'expid']:
hkeys.append(key)
fns = []
for oo in outputs:
fns.append([ooo.filename for ooo in oo['obs']])
# Trim first-frames from each EXPID:
for hh in hdrs:
firstframes = hh['dce']==0
for key in hkeys:
#print len(hh[key][True-firstframes])
try:
hh[key] = array(hh[key]).ravel()[True-firstframes]
except:
hh[key] = ones(min(targ_count)) -1.
scalefactors = scaleflux(hdrs, 'fitflux')
scalefactorsbg3 = scaleflux(hdrs, 'quadbg') #3-quadrant background
scalefactorsbgf = scaleflux(hdrs, 'fitback') #psf-fit background
scalefactorsbga = scaleflux(hdrs, 'bga') #aperture photometry background
h0 = int(hdrs[0]['hjd'][0])
hjd1 = array([h['hjd'] for h in hdrs])
if p.transit==True:
tephemeris = p.tt+ p.per*int((hjd1.min()-p.tt)/p.per)
rveph = p.tt+ p.per*int((hjd1.min()-p.tt)/p.per)
phase = (hjd1-tephemeris)/p.per
print 'transiting planet ephemeris for phase'
else:
rveph = p.rveph(h0) #an.rveph(p,h0)
phase = (hjd1-rveph)/p.per
print 'RV planet ephemeris for phase; RVeph=%f' % rveph
xoff = array([h['xoff'] for h in hdrs])
yoff = array([h['yoff'] for h in hdrs])
xoff = xoff - an.meanr(xoff, nsigma=sigma,niter=Inf)
yoff = yoff - an.meanr(yoff, nsigma=sigma,niter=Inf)
t1 = (hjd1 - hjd1.min()).ravel()
flux1 = array([h['fitflux']*s for h,s in zip(hdrs, scalefactors)])
bg31 = array([h['quadbg']*s for h,s in zip(hdrs, scalefactorsbg3)])
bgf1 = array([h['fitback']*s for h,s in zip(hdrs, scalefactorsbgf)])
bga1 = array([h['bga']*s for h,s in zip(hdrs, scalefactorsbga)])
# Detrend backgrounds:
bg3coef = an.polyfitr(t1,bg31.ravel(),1,sigma)
bgfcoef = an.polyfitr(t1,bgf1.ravel(),1,sigma)
bgacoef = an.polyfitr(t1,bga1.ravel(),1,sigma)
detbg3 = bg31.ravel() - polyval(bg3coef,t1) + mean(an.removeoutliers(bg31.ravel(),sigma))
detbgf = bgf1.ravel() - polyval(bgfcoef,t1) + mean(an.removeoutliers(bgf1.ravel(),sigma))
detbga = bga1.ravel() - polyval(bgacoef,t1) + mean(an.removeoutliers(bga1.ravel(),sigma))
fitfluxBgRatio3 = (flux1/bg31).ravel()
fitfluxBgRatiof = (flux1/bgf1).ravel()
fitfluxBgRatioa = (flux1/bga1).ravel()
ffbrcoef3 = an.polyfitr(t1,fitfluxBgRatio3,1,sigma)
ffbrcoeff = an.polyfitr(t1,fitfluxBgRatiof,1,sigma)
ffbrcoefa = an.polyfitr(t1,fitfluxBgRatioa,1,sigma)
ffbrCorr3 = fitfluxBgRatio3/polyval(ffbrcoef3,t1) * mean(an.removeoutliers(fitfluxBgRatio3,sigma))
ffbrCorrf = fitfluxBgRatiof/polyval(ffbrcoeff,t1) * mean(an.removeoutliers(fitfluxBgRatiof,sigma))
ffbrCorra = fitfluxBgRatioa/polyval(ffbrcoefa,t1) * mean(an.removeoutliers(fitfluxBgRatioa,sigma))
# For each key, output the binned-down version into the output object
out = genericObject()
for key in hkeys:
#print key, array([h[key] for h in hdrs]).shape
if hdrs[0].has_key(key) and len(hdrs[0][key])>0 and hdrs[0][key][0] is not None:
thisdata = array([h[key] for h in hdrs]).ravel()
hjd2,val2,ehjd2,err2 = tl.errxy(hjd1.ravel(),thisdata,tbins,clean=cleanopt,binfactor=binfactor)
exec ("out.%s = val2" % key) in None
exec ("out.e%s = err2" % key) in None
else:
exec ("out.%s = []" % key) in None
exec ("out.e%s = []" % key) in None
# if nbins<>None:
# tbins = linspace(hjd1.min()-1,hjd1.max()+1,nbins+1)
if returnstats:
bigtempval = tl.errxy(hjd1.ravel(),flux1.ravel(),tbins,clean=cleanopt,binfactor=binfactor,returnstats=returnstats)
hjd2,flux2,ehjd2,eflux2,cleanedstats = bigtempval
else:
hjd2,flux2,ehjd2,eflux2 = tl.errxy(hjd1.ravel(),flux1.ravel(),tbins,clean=cleanopt,binfactor=binfactor,returnstats=returnstats)
hjd2,bg32,ehjd2,ebg32 = tl.errxy(hjd1.ravel(),bg31.ravel(),tbins,clean=cleanopt,binfactor=binfactor)
hjd2,bgf2,ehjd2,ebgf2 = tl.errxy(hjd1.ravel(),bgf1.ravel(),tbins,clean=cleanopt,binfactor=binfactor)
hjd2,bga2,ehjd2,ebga2 = tl.errxy(hjd1.ravel(),bga1.ravel(),tbins,clean=cleanopt,binfactor=binfactor)
hjd2,detbg32,ehjd2,edetbg32 = tl.errxy(hjd1.ravel(),detbg3.ravel(),tbins,clean=cleanopt,binfactor=binfactor)
hjd2,detbgf2,ehjd2,edetbgf2 = tl.errxy(hjd1.ravel(),detbgf.ravel(),tbins,clean=cleanopt,binfactor=binfactor)
hjd2,detbga2,ehjd2,edetbga2 = tl.errxy(hjd1.ravel(),detbga.ravel(),tbins,clean=cleanopt,binfactor=binfactor)
hjd2,ffbrCorr32,ehjd2,effbrCorr32 = tl.errxy(hjd1.ravel(),ffbrCorr3,tbins,clean=cleanopt,binfactor=binfactor)
hjd2,ffbrCorrf2,ehjd2,effbrCorrf2 = tl.errxy(hjd1.ravel(),ffbrCorrf,tbins,clean=cleanopt,binfactor=binfactor)
hjd2,ffbrCorra2,ehjd2,effbrCorra2 = tl.errxy(hjd1.ravel(),ffbrCorra,tbins,clean=cleanopt,binfactor=binfactor)
hjd2,xoff2,ehjd2,exoff2 = tl.errxy(hjd1.ravel(),xoff.ravel(),tbins,clean=cleanopt,binfactor=binfactor)
hjd2,yoff2,ehjd2,eyoff2 = tl.errxy(hjd1.ravel(),yoff.ravel(),tbins,clean=cleanopt,binfactor=binfactor)
if p.transit==True:
tephemeris = p.tt+ p.per*int((hjd1.min()-p.tt)/p.per)
phase2 = (hjd2-tephemeris)/p.per
print 'transiting planet ephemeris for phase'
else:
phase2 = (hjd2-rveph)/p.per
print 'RV planet ephemeris for phase; RVeph=%f' % rveph
t2 = hjd2 - int(hjd1.min())
# Write other pertinent values to the object
out.planet = p
out.dataset = outputs[0]['dataset']
out.phase = phase2
out.ephase = ehjd2/p.per
out.h0 = int(hjd1.min())
out.t = hjd2-out.h0
out.flux = flux2
out.eflux = eflux2
out.bg = bg32
out.ebg = ebg32
out.bg2 = bgf2
out.ebg2 = ebgf2
out.bg3 = bga2
out.ebg3 = ebga2
out.detbg = detbg32
out.edetbg = edetbg32
out.detbg2 = detbgf2
out.edetbg2 = edetbgf2
out.detbg3 = detbga2
out.edetbg3 = edetbga2
out.detratio = ffbrCorr32
out.edetratio = effbrCorr32
out.detratio2 = ffbrCorrf2
out.edetratio2 = effbrCorrf2
out.detratio3 = ffbrCorra2
out.edetratio3 = effbrCorra2
out.keys = hkeys
out.filename = fns
#out.hjd = hjd2
#out.ehjd = ehjd2
#out.xoff = xoff2
#out.exoff = exoff2
#out.yoff = yoff2
#out.eyoff = eyoff2
if returnstats:
ret = out, (init_count, targ_count, cleanedstats)
else:
ret = out
return ret
[docs]def mipsflat(files, norm=True, nmax=99999, retall=False):
"""
Make a flat field from a set of MIPS filenames. Use a stack
median. In creating the flat only use the half-side of the detector
that the target star is not falling on.
If nmax is too large, you may run out of memory. Set it lower!
:INPUTS:
norm=True -- whether to normalize the final flat by the chip-wide median.
nmax -- the number of files to use, per
:EXAMPLE:
::
init = spitzer.initobs('upsand2')
files = [init.path+f for f in init.files[0]]
flat = spitzer.mipsflat(files)
"""
#2010-06-11 09:37 IJC: Created
import spitzer, phot
from numpy import array, zeros, median,mean,sort
# Load the filenames and parameters, and initialize:
fileids = array([spitzer.mipsfile(f)[1:4] for f in files])
aor = fileids[:,0]
dce = fileids[:,1]
expid = fileids[:,2]
lhs = (dce % 2)==0
rhs = (dce % 2)==1
nfiles = len(fileids)
nmax = min(nfiles,nmax)
nflats = 1.*nfiles/nmax
full_flat = zeros((nflats,128,128),float)
nper_flat = []
#pdb.set_trace()
# Create the sub-flat-fields
for ii in range(nflats):
start = ii*nmax
end = min((ii+1)*nmax,nfiles)
stack = phot.subreg2(files[start:end],center=[64,64],dim=[128,128])
left_flat = median(stack[lhs[start:end],:,:],0)
right_flat = median(stack[rhs[start:end],:,:],0)
full_flat[ii,:,0:64] = right_flat[:,0:64]
full_flat[ii,:,64:128] = left_flat[:,64:128]
nper_flat.append(stack.shape[0])
# Properly weight the sub-flats, in case the last one used fewer files:
nper_flat = array(nper_flat).reshape(nflats,1,1)
final_flat = (full_flat*nper_flat).sum(0)/nper_flat.sum()
if norm:
final_flat = final_flat/median(sort(final_flat.ravel()))
if retall:
ret = final_flat, full_flat
else:
ret = final_flat
return ret
[docs]def dither14(cparams):
"Returns: cparams[0:14] as an arra of shape (14,1)"
# 2011-10-27 16:04 IJMC: Created
from numpy import reshape
return reshape(cparams[0:14], (14,1))
[docs]def feps_flux_loss(offset, version='S18.18.0', order='SL1', path='~/proj/spitzer/data/irs/jeroen_20131121/', verbose=False):
"""
Compute flux throughput of the IRS instrument using a lookup table.
:INPUTS:
offset: scalar
offset of target point source from the center of slit [units???]
version: string
order: string
path: string
Directory in which the lookup table is located.
Alternatively, this can be the cal file in question, loaded by
idlsave.read('IRSX_%s_%s_flcf.sav')
:EXAMPLE:
::
import tools
offset_sl1 = 0.0
wave, fraction = spitzer.feps_flux_loss(0, version="S18.18.0",order='SL1')
"""
#2013-11-22 IJMC: Received from J. Bouwman. Added minimal comments.
#2013-12-02 00:08 IJMC: Converted from IDL to Python. Replaced
# call to FEPS_INTERPOL with numpy.interp()
#2014-04-30 17:06 IJMC: Can now pass in pre-loaded cal file via 'path'.
import idlsave # or: from scipy import io
import numpy as np
from tools import feps_interpol
import os
if not isinstance(offset, np.ndarray):
offset = np.array([offset]).ravel()
if isinstance(path, str):
path = os.path.expanduser(path)
cal_file = idlsave.read('%s/IRSX_%s_%s_flcf.sav' % (path,order, version) , verbose=verbose)
else:
cal_file = path
xoffset = cal_file['flcf']['xoffset'][0]
wave = cal_file['flcf']['wave'][0]
flcf = cal_file['flcf']['flcf'][0]
fraction = np.zeros(wave.size)
for iwave in xrange(wave.size):
#fraction[iwave] = feps_interpol(xoffset, flcf[:,iwave], offset, linear=True)
fraction[iwave] = np.interp(offset, xoffset, flcf[:, iwave])
return wave, fraction
[docs]def model_irs_slitloss(offsets, nobs, nchan, normalize=False, version='S18.18.0', order='SL1', path='~/proj/spitzer/data/irs/jeroen_20131121/', verbose=False, cal_file=None):
"""
Compute flux throughput of the IRS instrument using a lookup table.
:INPUTS:
offset: 1D NumPy array
offset of target point source from the center of slit [units???]
nobs: int
number of obesrvations
nchan: int
number of wavelength channels
normalize: bool
Whether to divide result by stack median
cal_file: None, 2-tuple, or IDLSave object
pass calibration file in directly (or tuple [offset, flcf]),
to save time in a highly repetitive analysis.
:EXAMPLE:
::
import spitzer
slitloss = spitzer.model_irs_slitloss(y_offset_pu_a, nobs, nchan, False)
"""
# 2013-12-10 16:06 IJMC: Created.
import numpy as np
from tools import feps_interpol
import os
if not isinstance(offsets, np.ndarray):
offsets = np.array([offsets]).ravel()
if isinstance(cal_file, tuple):
xoffset, flcf = cal_file[0:2]
else:
if cal_file is None:
import idlsave # or: from scipy import io
path = os.path.expanduser(path)
cal_file = idlsave.read('%s/IRSX_%s_%s_flcf.sav' % (path,order, version) , verbose=verbose)
xoffset = cal_file['flcf']['xoffset'][0]
flcf = cal_file['flcf']['flcf'][0]
slitloss = np.zeros((nobs, nchan), dtype=float)
offset_inds = np.array([((xoffset - offsets[iobs]) < 0).nonzero()[0][-1] for iobs in xrange(nobs)])
xoi = xoffset[offset_inds]
doff = (xoi - xoffset[offset_inds+1])
frac0 = (1-(offsets - xoi) / doff).reshape(nobs, 1)
slitloss = flcf[offset_inds-1] * (1.-frac0) + flcf[offset_inds] * (frac0)
if normalize:
ret = slitloss / np.median(slitloss, 0)
else:
ret = slitloss
return ret
[docs]def model_irs_eclipse(params, nobs, nchan, rampfunc, npar_ramp, tparams, time, phase, normalize=False, cal_file=None):
"""
Generate a full spectrophotometric model of an IRS-observed eclipse.
:INPUTS:
offset: 1D NumPy array
offset of target point source from the center of slit [units???]
nobs: int
number of obesrvations
nchan: int
number of wavelength channels
rampfunc: function
Ramp function: e.g., :func:`phasecurves.ramp5n`, et al.
npar_ramp: int
Number of parameters required for rampfunc.
tparams: sequence
Transiting planet parameters, to be passed to
:func:`transit.modeleclipse_simple`
time: 1D NumPy array
Time of observations, to be passed to
:func:`transit.modeleclipse_simple`
phase: 1D NumPy array
Orbital phase of observations, to be passed to 'rampfunc.'
normalize: bool
Whether to normalize result by stack median.
cal_file: None or IDLSave object
pass calibration file in directly, to :func:`model_irs_slitloss`
:EXAMPLE:
::
test = model_irs_eclipse(params, nobs, nchan, rampfunc, npar_ramp, transit_params, this_mbjd, phase, )
"""
# 2013-12-11 14:30 IJMC: Created.
from numpy import hstack, vstack, median, ones
from transit import modeleclipse_simple, occultuniform
offsets = params[0:nobs]
ramp_param = params[nobs:nobs+npar_ramp*nchan].reshape(nchan, npar_ramp)
eclipse_param = hstack((params[nobs+npar_ramp*nchan:].reshape(nchan, 2), ones((nchan,1))))
slitloss = model_irs_slitloss(offsets, nobs, nchan, normalize=False, cal_file=cal_file)
ramps = vstack([rampfunc(par, phase) for par in ramp_param]).T
lightcurves = vstack([modeleclipse_simple(epar, tparams, occultuniform, time) for epar in eclipse_param]).T
ret = slitloss * lightcurves * ramps
if normalize:
ret /= median(ret, 0)
return ret