Source Code for Module dark_models

  1  """Define models to represent the observed ACA CCD dark current distribution. 
  2  Provide a function to degrade the CCD for a delta time.""" 
  3   
  4  import numpy as np 
  5  from numpy import exp, log, arange 
  6   
  7  # Define a common fixed binning of dark current distribution 
  8  from darkbins import x0, x1, dx, bins as xbins 
  9   
 10  # Some constants and globals.  Done this way to support sherpa fitting. 
 11  # Needs to be re-worked to be nicer. 
 12   
 13  # Fixed gaussian for smoothing the broken power law 
 14  dx = 0.1 
 15  sigma = 0.25                            # Gaussian sigma in log space 
 16  xg = arange(-2.5*sigma, 2.5*sigma, dx, dtype=float) 
 17  yg = exp(-0.5 * (xg/sigma)**2) 
 18  yg /= np.sum(yg) 
 19   
 20  NPIX = 1024**2 
 21   
 22  # Fixed 
 23  xall = (xbins[:-1] + xbins[1:]) / 2.0 
 24  imin = 0 
 25  imax = len(xall) 
 26   
 27 -def broken_pow(pars, x): 
 28      """Broken power-law.  Pars are same as bpl1d: 
 29      1: gamma1 
 30      2: gamma2 
 31      3: x_b (break point) 
 32      4: x_r (normalization reference point) 
 33      5: ampl1""" 
 34      (gamma1, gamma2, x_b, x_r, ampl1) = pars 
 35      ampl2 = ampl1 * (x_b / x_r)**(gamma2 - gamma1) 
 36      ok = x > x_b 
 37      y = ampl1 * (x / x_r)**(-gamma1) 
 38      y[ok] = ampl2 * (x[ok] / x_r)**(-gamma2) 
 39      return y 
 40   
 41 -def smooth_broken_pow(pars, x): 
 42      """Smoothed broken power-law.  Pars are same as bpl1d (NOT + gaussian sigma): 
 43      1: gamma1 
 44      2: gamma2 
 45      3: x_b (break point) 
 46      4: x_r (normalization reference point) 
 47      5: ampl1 
 48      #   NOT 6: sigma (bins)""" 
 49      (gamma1, gamma2, x_b, x_r, ampl1) = pars 
 50      ampl2 = ampl1 * (x_b / x_r)**(gamma2 - gamma1) 
 51      ok = xall > x_b 
 52      y = ampl1 * (xall / x_r)**(-gamma1) 
 53      y[ok] = ampl2 * (xall[ok] / x_r)**(-gamma2) 
 54      return np.convolve(y, yg, mode='same')[imin:imax] 
 55   
 56 -def smooth_twice_broken_pow(pars, x): 
 57      """Smoothed broken power-law.  Pars are same as bpl1d (NOT + gaussian sigma): 
 58      1: gamma1 
 59      2: gamma2 
 60      3: x_b (break point) 
 61      4: x_r (normalization reference point) 
 62      5: ampl1 
 63      #   NOT 6: sigma (bins)""" 
 64      gamma0 = -4.0 
 65      x_b0 = 5 
 66   
 67      (gamma1, gamma2, x_b, x_r, ampl1) = pars 
 68      y = ampl1 * (x / x_r)**(-gamma1) 
 69   
 70      ok = x > x_b 
 71      ampl2 = ampl1 * (x_b / x_r)**(gamma2 - gamma1) 
 72      y[ok] = ampl2 * (x[ok] / x_r)**(-gamma2) 
 73   
 74      ok = x < x_b0 
 75      ampl0 = ampl1 * (x_b0 / x_r)**(gamma0 - gamma1) 
 76      y[ok] = ampl0 * (x[ok] / x_r)**(-gamma0) 
 77   
 78      return np.convolve(y, yg, mode='same') 
 79   
 80 -def zero_dark_ccd(): 
 81      return np.zeros(NPIX) 
 82   
 83 -def pristine_ccd(): 
 84      """Generate and return a 'pristine' (zero-year) dark current map.  This was 
 85      empirically derived.""" 
 86      nlow = 350000 
 87      dark = np.append(np.random.lognormal(log(4.5), 0.25, NPIX-nlow), 
 88                       np.random.lognormal(log(2.5), 0.45, nlow)) 
 89      return dark 
 90   
 91 -def nompars(dyear): 
 92      """Return nominal degradation parameters for dyear interval""" 
 93      gamma1 = 0.14 
 94      gamma2 = 3.15 
 95      x_b = 125. 
 96      x_r = 155. 
 97      ampl1 = 1032. * dyear 
 98      return (gamma1, gamma2, x_b, x_r, ampl1) 
 99   
100 -def degrade_ccd(dark, dyear): 
101      """Degrade CCD (dark current map) for dyear years.  The input 'dark' map is 
102      updated in place.""" 
103      pars = nompars(dyear) 
104      darkmodel = smooth_twice_broken_pow(pars, xall) 
105   
106      ndegrade = 800000 
107      ipix = np.random.randint(0, NPIX, ndegrade) 
108      dark[ipix] += np.random.poisson(dyear * 0.7, ndegrade) 
109   
110      darkran = np.random.poisson(darkmodel) 
111      for i, n in enumerate(darkran): 
112          # Generate n log-uniform variates within bin 
113          if n > 0: 
114              logdark = np.random.uniform(log(xbins[i]), log(xbins[i+1]), n) 
115              ipix = np.random.randint(0, NPIX, n) 
116              dark[ipix] += exp(logdark) 
117   
118 -def temp_scalefac(T_ccd): 
119      """Return the multiplicative scale factor to convert a CCD dark map from 
120      the nominal -19C temperature to the temperature T.  Based on an overall reduction 
121      of 0.62 after changing from -15 to -19.   
122      """ 
123      return exp(log(0.62) / 4.0 * (-19.0 - T_ccd)) 
124