SHARE
TWEET

Untitled

a guest Sep 20th, 2019 93 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. def get_delensed_amanda(tt_unlen, te_unlen, ee_unlen, lpp_theo, clpp_theo, L, N_L, klmin, klmax, load_dls=True, unit='uk', lmax_ret=8000):
  2.     import camb.correlations
  3.     import numpy as np
  4.  
  5.     kappa_lmin = int(klmin)
  6.     kappa_lmax = int(klmax)
  7.  
  8.     reload(np)  # necessary to call this function multiple times (?)
  9.  
  10.     # want theory to start at l = 0
  11.     lpp_th = np.zeros(len(lpp_theo)+int(lpp_theo[0]))
  12.     lpp_th[0] = 0.
  13.     lpp_th[1] = 1.
  14.     lpp_th[int(lpp_theo[0]):] = np.array(lpp_theo)
  15.     # fill in unlensed theory and lensing potential
  16.     clpp_th = np.zeros(len(lpp_th))
  17.     tt_th = np.zeros(len(lpp_th))
  18.     te_th = np.zeros(len(lpp_th))
  19.     ee_th = np.zeros(len(lpp_th))
  20.     idx   = int(lpp_theo[0])
  21.     clpp_th[idx:] = np.array(clpp_theo)
  22.     tt_th[idx:] = np.array(tt_unlen)
  23.     te_th[idx:] = np.array(te_unlen)
  24.     ee_th[idx:] = np.array(ee_unlen)
  25.     # fill in the noise curve so it also starts at L=0 and is same size as theory
  26.     NL = np.zeros(lpp_th.shape)
  27.     idx = int(L[0])
  28.     NL[:idx] = np.inf
  29.     if len(N_L) >= len(NL) + idx:
  30.         NL[idx:] = np.array(N_L[:len(NL)+idx])
  31.     else:
  32.         NL[idx:len(N_L)+idx] = np.array(N_L)
  33.         NL[len(N_L)+idx:] = np.inf
  34.     L = lpp_th.copy()  
  35.  
  36.     def _weiner1d_filter(signal, noise):
  37.         return (signal/(signal+noise))*signal
  38.  
  39.     # get residual lensing potential
  40.     clkk           = clpp_th*(2*np.pi/4)
  41.     clkk_filt      = _weiner1d_filter(clkk, NL)
  42.     clpp_filt      = clkk_filt*4/(2*np.pi)
  43.     loc            = np.where(L > kappa_lmax)
  44.     if len(loc[0]) > 0:
  45.         if loc[0][0] > len(clpp_filt):
  46.             clpp_filt[loc[0][0]:] = 0.
  47.     loc            = np.where(L < kappa_lmin)
  48.     clpp_filt[loc] = 0.
  49.  
  50.     clpp_res = clpp_th-clpp_filt
  51.  
  52.     # create array to hold delensed spectra
  53.     dls = np.zeros((len(lpp_th),4))
  54.     dls[:,0] = tt_th
  55.     dls[:,1] = ee_th
  56.     dls[:,2] = np.zeros(len(tt_th))
  57.     dls[:,3] = te_th
  58.  
  59.     reload(camb)
  60.     reload(camb.correlations)
  61.     camb_ret = camb.correlations.lensed_cls(dls, clpp_res, delta_cls=False, sampling_factor=2.)
  62.  
  63.     ret         = {}
  64.     ret['dltt'] = camb_ret[:,0]
  65.     ret['dlee'] = camb_ret[:,1]
  66.     ret['dlbb'] = camb_ret[:,2]
  67.     ret['dlte'] = camb_ret[:,3]
  68.    
  69.     return ret['dltt'][1:], ret['dlte'][1:], ret['dlee'][1:]
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top