Doc notes

[cmbant]

I managed to run it on windows via Ubuntu WSL!

Things that were a bit confusing:

Starting point: whether to use plancklens or lensit to make noise forecasts..

qe_key is missing from https://plancklens.readthedocs.io/en/latest/qresp.html Valid values for "ksource", "qe_key" are not obviously documented.

plot_noiselevels is not advertised on home page. Would be easier to just provide a general more-obviously-named N0 function. I made this (which is not totally general -not all keys/legs and no general N_L - but seems to work)

def get_N0(beam_fwhm = 1.4, nlev_t = 5.,nlev_p = None,  lmax_CMB= 3000, 
       lmin_CMB =100, lmax_out=None,   cls_len = None, cls_weight = None, 
       joint_TP=True, ksource = 'p', fname = None): 

      if nlev_p is None:
          nlev_p = nlev_t*np.sqrt(2)

      lmax_ivf = lmax_CMB
      lmin_ivf = lmin_CMB
      lmax_qlm = lmax_out or lmax_ivf

      cls_path = os.path.join(os.path.dirname(os.path.abspath(plancklens.__file__)), 'data', 'cls')

      cls_len = cls_len or utils.camb_clfile(os.path.join(cls_path, 'FFP10_wdipole_lensedCls.dat'))
      cls_weight = cls_weight or utils.camb_clfile(os.path.join(cls_path, 'FFP10_wdipole_lensedCls.dat'))

      assert ksource in ['p', 'f']
      qe_keys = [ksource + 'tt', ksource+'_p'] 
      if not joint_TP:
          qe_keys.append(ksource)

      transf = hp.gauss_beam(beam_fwhm / 60. / 180. * np.pi, lmax=lmax_ivf)

      Noise_L_T =  (nlev_t / 60. /180. * np.pi) ** 2 / transf ** 2
      Noise_L_P =  (nlev_p / 60. /180. * np.pi) ** 2 / transf ** 2

      cls_dat = {
          'tt': (cls_len['tt'][:lmax_ivf + 1] + Noise_L_T),
          'ee': (cls_len['ee'][:lmax_ivf + 1] + Noise_L_P),
          'bb': (cls_len['bb'][:lmax_ivf + 1] + Noise_L_P),
          'te':  np.copy(cls_len['te'][:lmax_ivf + 1]) }

      # 1/(C+N) filter spectra
      fal_sepTP = {spec: utils.cli(cls_dat[spec]) for spec in ['tt','ee','bb']}

      cls_ivfs_sepTP = {'tt':fal_sepTP['tt'].copy(),
                        'ee':fal_sepTP['ee'].copy(),
                        'bb':fal_sepTP['bb'].copy(),
                        'te':cls_len['te'][:lmax_ivf + 1] * fal_sepTP['tt'] * fal_sepTP['ee']}

      fal_jtTP = utils.cl_inverse(cls_dat)
      cls_ivfs_jtTP = utils.cl_inverse(cls_dat)

      for cls in [fal_sepTP, fal_jtTP, cls_ivfs_sepTP, cls_ivfs_jtTP]:
          for cl in cls.values():
              cl[:max(1, lmin_ivf)] *= 0.

      N0s = {}
      N0_curls = {}
      for qe_key in qe_keys:
          NG, NC, NGC, NCG = nhl.get_nhl(qe_key, qe_key, cls_weight, cls_ivfs_sepTP, lmax_ivf, lmax_ivf, lmax_out=lmax_qlm)
          RG, RC, RGC, RCG = qresp.get_response(qe_key, lmax_ivf, ksource, cls_weight, cls_len, fal_sepTP, lmax_qlm=lmax_qlm)

          N0s[qe_key] = utils.cli(RG ** 2) * NG
          N0_curls[qe_key]= utils.cli(RC ** 2) * NC

      if joint_TP:
          NG, NC, NGC, NCG = nhl.get_nhl(ksource, ksource, cls_weight, cls_ivfs_jtTP, lmax_ivf, lmax_ivf, lmax_out=lmax_qlm)
          RG, RC, RGC, RCG = qresp.get_response(ksource, lmax_ivf, ksource, cls_weight, cls_len, fal_jtTP, lmax_qlm=lmax_qlm)
          N0s[ksource] = utils.cli(RG ** 2) * NG
          N0_curls[ksource] = utils.cli(RC ** 2) * NC

      return N0s, N0_curls

[carronj]

Thanks Antony for this. I have updated the readme, the doc and added a n0s.py module with your script and the iterative N0 calculation.