Source Code for Module wfpc2tools.wfpc2cte

  1  """ WFPC2CTE - Module for computing the CTE degradation for WFPC2 images 
  2   
  3  This module updates the header of the input WFPC2 image with standardized 
  4  computations of the effect of CTE based on the algorithm published by Dolphin 
  5  (2004, http://purcell.as.arizona.edu/wfpc2_calib/2004_12_20.html).   
  6   
  7  ASSUMPTIONS for the COMPUTATION:  
  8      1. The CTE gets computed for a source at the chip center. 
  9      2. The background (in electrons) gets defined by the clipped mode of the  
 10          central 200x200 pixels from the image. 
 11      3. The source is assumed to have 100 electrons, 1000 electrons and  
 12          10000 electrons in the aperture. 
 13      4. The reported CTE is the sum of the XCTE and YCTE computed from  
 14          Dolphin's algorithm. 
 15   
 16  INPUT: 
 17  The sole input for this task is the filename of the WFPC2 image.  
 18       
 19      If the input image is in GEIS format, it will convert it to 
 20      a multi-extension FITS formatted file, then update the FITS file while 
 21      leaving the GEIS image un-modified.   
 22       
 23      If the input image is already multi-extension FITS, it will update the 
 24      header directly. 
 25       
 26      If the input image is waivered FITS, it will quit with a message  
 27      telling the user to first convert the file to GEIS. The user can then 
 28      provide the GEIS image as input. 
 29   
 30  OUTPUT:     
 31  The keywords which get updated are: 
 32      CTE_1E2  - CTE for a source with an intensity of 100 electrons  
 33      CTE_1E3  - CTE for a source with an intensity of 1000 electrons 
 34      CTE_1E4  - CTE for a source with an intensity of 10000 electrons 
 35   
 36  SYNTAX: 
 37  This task can be run on an input WFPC2 image using either of the following  
 38  calls: 
 39      wfpc2cte.compute_CTE(filename,quiet=True) 
 40    -or- 
 41      wfpc2cte.run(filename,quiet=True) 
 42   
 43  where the filename is the name of the input WFPC2 image. 
 44   
 45   
 46  EXAMPLE: 
 47   
 48  The syntax for running this task on a WFPC2 file named 'u40x0102m.c0h': 
 49      import wfpc2cte 
 50      wfpc2cte.run('u40x0102m.c0h') 
 51       
 52  The command to print out this help file: 
 53      wfpc2cte.help()  
 54  """ 
 55  from __future__ import division # confidence medium 
 56   
 57  import numpy as np 
 58  import pyfits 
 59  import pytools 
 60  from pytools import readgeis,fileutil 
 61  import imagestats 
 62  __version__ = '1.2.4 (2-July-2008)' 
 63   
 64  # This contains the default values in electrons for the CTE sources 
 65  DEFAULT_COUNTS = np.array([100,1000,10000],np.float32) 
 66   
 67  chip_bg_factor = 10 
 68  chip_lbg_factor = 1 
 69  chip_lct_factor = 7 
 70   
 71 -def compute_chip_values(extn,gain,nclip=3): 
 72      chip_values = {} 
 73       
 74      # Determine the center of the chip as (Y,X) 
 75      chip_shape = extn.data.shape 
 76      chip_center = (chip_shape[0]/2.,chip_shape[1]/2.) 
 77      chip_values['center'] = (chip_center[0] / 800.,chip_center[1]/800.) 
 78      center_slice = [slice(chip_center[0]-100.,chip_center[0]+100.), 
 79                      slice(chip_center[1]-100.,chip_center[1]+100.)] 
 80      center_region = extn.data[center_slice] 
 81      # Compute the background as the clipped mode of the region. 
 82      chip_stats = imagestats.ImageStats(center_region,fields='mode', 
 83                                      lower=-99.0,upper=4096., 
 84                                      nclip=nclip,binwidth=0.01) 
 85      if chip_stats.mode >= 0: 
 86          chip_background = np.sqrt(np.power(chip_stats.mode,2) + 1) 
 87      else: 
 88          chip_background = 1.0 
 89      chip_background *= gain 
 90      chip_values['bg_raw'] = chip_background 
 91      chip_values['lbg'] = np.log(chip_background) - chip_lbg_factor 
 92      chip_values['bg'] = chip_background - chip_bg_factor     
 93   
 94      # Compute e^(counts) for all the cases 
 95      chip_values['lct'] = np.log(DEFAULT_COUNTS) - chip_lct_factor 
 96       
 97      return chip_values 
 98       
 99 -def compute_XCTE(xpos,bg): 
100   
101      return 0.0021*np.exp(-0.234*bg)*xpos 
102       
103 -def compute_YCTE(chip_values,yr,xcte): 
104   
105      lbg = chip_values['lbg'] 
106      bg = chip_values['bg']                                 
107      ypos = chip_values['center'][0] 
108   
109      #adjust original source counts by XCTE 
110      lct = chip_values['lct'] 
111      lct += 0.921*xcte 
112   
113      c1 = 0.0114*(0.670*np.exp(-0.246*lbg)+0.330*np.exp(-0.0359*bg))*(1.+0.335*yr-0.0074*yr*yr)*ypos 
114      c2 = 3.55*np.exp(-0.474*lct) 
115      ycte = np.log(np.exp(c1)*(1+c2)-c2)/0.436  
116   
117      return ycte 
118       
119       
120 -def update_CTE_keywords(hdr, cte,quiet=False,update=True): 
121      # Start by checking to see if the keywords to be updated already exist 
122      # If not, print a warning and insure quiet=False so the results get 
123      # reported to STDOUT at the very least. 
124      if not hdr.has_key('CTE_1E2'): 
125          print "WARNING: CTE keywords not found in %s,%d header."%(hdr['extname'],hdr['extver']) 
126          print "         Adding the keywords exclusively to " 
127          print "         the FITS file's extension header." 
128          quiet = False 
129   
130      if update: 
131          hdr.update('CTE_1E2',cte[0]) 
132          hdr.update('CTE_1E3',cte[1],after='CTE_1E2') 
133          hdr.update('CTE_1E4',cte[2],after='CTE_1E3') 
134   
135      if not quiet: 
136          print 'CTE_1E2   = ',cte[0] 
137          print 'CTE_1E3  = ',cte[1] 
138          print 'CTE_1E4 = ',cte[2] 
139       
140 -def compute_CTE(filename,quiet=True,nclip=3,update=True): 
141   
142      newname = None 
143      if filename.find('.fits') < 0: 
144          # We are working with a GEIS image 
145          newname = fileutil.buildFITSName(filename) 
146      else: 
147          # We are working with a FITS image, so update it directly 
148          newname = filename 
149   
150      update_mode = 'update' 
151      # Open the image in update mode. 
152      # If it is a GEIS image on input, convert to FITS using 'newname'    
153      # then update the FITS file only... 
154      fimg = fileutil.openImage(filename,mode=update_mode,fitsname=newname) 
155      fimg.info() 
156      if isinstance(fimg[1],pyfits.TableHDU): 
157          fimg.close() 
158          print 'Input image is in the unsupported "waivered" FITS format. ' 
159          print '   Please convert to GEIS or multi-extension FITS format.' 
160          print 'Exiting...' 
161          return 
162            
163      obs_date = fimg[0].header['expstart'] 
164      obs_mjd = (obs_date - 50193)/365.25 
165   
166      obs_gain = fimg[0].header['atodgain'] 
167   
168      for extn in fimg[1:]: 
169          if extn.header['extname'] == 'SCI': 
170              chip_values = compute_chip_values(extn,obs_gain,nclip=nclip) 
171                   
172              # Compute XCTE and YCTE for all sources 
173              xcte = compute_XCTE(chip_values['center'][1],chip_values['bg']) 
174              ycte = compute_YCTE(chip_values, obs_mjd,xcte) 
175              if not quiet: 
176                  print 'Background computed to be: ',chip_values['bg_raw'] 
177               
178              # Based on Workshop 2002 paper, after equation 8... 
179              total_cte = xcte + ycte 
180   
181              update_CTE_keywords(extn.header, total_cte,quiet=quiet,update=update) 
182   
183      if not quiet and update: 
184          print 'Updating keywords in: ',newname 
185           
186      fimg.flush() 
187   
188      # Close the file 
189      fimg.close() 
190              
191 -def run(filename,quiet=True,nclip=3): 
192      compute_CTE(filename,quiet=quiet,nclip=nclip) 
193   
194 -def help(): 
195      print __doc__ 
196