""" CMB-S4 LAT Survey Simple Noise Model - v181016 This product has been brought to you by the Noise Tiger Team. Based on the Simons Observatory V3 Noise Calculator (public version). See the accompanying PDF for modeling details. """ from __future__ import print_function import numpy as np import matplotlib matplotlib.use('pdf') matplotlib.rc('font', family='serif', serif='cm10') matplotlib.rc('text', usetex=True) fontProperties = {'family':'sans-serif', 'weight' : 'normal', 'size' : 16} import matplotlib.pyplot as plt def get_atmosphere_C(freqs, version=1, el=None): """ Returns atmospheric noise power at ell=1000, for an ACTPol optics tube. In units of [uK^2 sec]. This only works for a few special frequencies. Basic model assumes el=50. A simple rescaling (proportional to csc(el)) is applied for other values of el. version=0: This atmospheric model was used in SO V3 forecasts but contains an error. version=1: This atmospheric model is better than version=0, in that at least one error has been corrected. The relative atmospheric powers have been measured in AT model, and then calibrated to ACT. Low frequency results are inflated somewhat because ACT sees more power at 90 GHz than predicted by this modeling. """ # This el_correction is quite naive; atmosphere may be less # structured (i.e. smoothed out) at lower elevation because of the # increased thickness. if el is None: el = 50. el_correction = np.sin(50.*np.pi/180) / np.sin(el*np.pi/180) if version == 0: data_bands = np.array([ 27., 39., 93., 145., 225., 280.]) data_C = np.array([200., 77., 1800., 12000., 68000., 124000.]) data = {} for b,C in zip(data_bands, data_C): data[b] = C # Jam in a fake value for 20 GHz. data[20.] = 200. elif version == 1: data_bands = np.array([ 20., 27., 39., 93., 145., 225., 280.]) data_C = np.array([ 3.45565594e+02 * 1.4, 6.00734484e+01 * 1.4, 1.45291936e+01 * 1.4, 6.81426710e+02 * 1.4, 1.20000000e+04, 2.66015359e+05, 1.56040514e+06,]) data = {} for b,C in zip(data_bands, data_C): data[b] = C return np.array([data[f] * el_correction**2 for f in freqs]) """See lower for subclasses of SOLatType -- instrument specific parameters are configured in __init__. """ class SOLatType: def __init__(self, *args, **kwargs): raise RuntimeError('You should subclass this.') def get_bands(self): return self.bands.copy() def get_beams(self): return self.beams.copy() def precompute(self, N_tubes, N_tels=1): if self.el is not None: el_lims = el_data['valid'] assert(el_lims[0] <= self.el) and (self.el <= el_lims[1]) band_idx = np.array([np.argmin(abs(el_data['bands'] - b)) for b in self.bands]) assert(np.all(abs(np.array(el_data['bands'])[band_idx] - self.bands) < 5)) white_noise_el_rescale = np.array( [np.polyval(el_data['polys'][i][::-1], self.el) for i in band_idx]) print('Elevation mapping speed hits: ', white_noise_el_rescale**-2) else: white_noise_el_rescale = np.array([1.] * len(self.bands)) # Accumulate total white noise level and atmospheric # covariance matrix for this configuration. band_weights = np.zeros(self.n_bands) for (tube_name, tube_count) in N_tubes: tube_noise = self.tube_configs[tube_name] * white_noise_el_rescale s = (tube_noise != 0) print(s, tube_noise) band_weights[s] += tube_count * N_tels * tube_noise[s]**-2 self.band_sens = np.zeros(self.n_bands) + 1e9 s = (band_weights > 0) self.band_sens[s] = band_weights[s]**-0.5 # Special for atmospheric noise model. self.Tatmos_C = get_atmosphere_C(self.bands, self.atm_version, el=self.el) * self.FOV_mod self.Tatmos_ell = 1000. + np.zeros(self.n_bands) self.Tatmos_alpha = -3.5 + np.zeros(self.n_bands) # Compute covariant weight matrix (atmosphere parameters). cov_weight = np.zeros((self.n_bands,self.n_bands)) pcov_weight = np.zeros((self.n_bands,self.n_bands)) atm_rho = 0.9 for (tube_name, tube_count) in N_tubes: # Get the list of coupled bands; e.g. [1,2] for MF. nonz = self.tube_configs[tube_name].nonzero()[0] for i in nonz: for j in nonz: w = {True: 1., False: atm_rho}[i==j] assert(cov_weight[i,j] == 0.) # Can't do overlapping # tubes without weights. cov_weight[i,j] += tube_count * N_tels / ( w * ( self.Tatmos_C[i] * self.Tatmos_C[j])**.5 ) pcov_weight[i,j] = w # Reciprocate non-zero elements. s = (cov_weight!=0) self.Tatmos_cov = np.diag([1e9]*self.n_bands) self.Tatmos_cov[s] = 1./cov_weight[s] # Polarization is simpler... self.Patmos_ell = 700. + np.zeros(self.n_bands) self.Patmos_alpha = -1.4 + np.zeros(self.n_bands) self.Patmos_cov = pcov_weight def get_survey_time(self): t = self.survey_years * 365.25 * 86400. ## convert years to seconds return t * self.survey_efficiency def get_survey_spread(self, f_sky, units='arcmin2'): # Factor that converts uK^2 sec -> uK^2 arcmin^2. A = f_sky * 4*np.pi if units == 'arcmin2': A *= (60*180/np.pi)**2 elif units != 'sr': raise ValueError("Unknown units '%s'." % units) return A / self.get_survey_time() def get_white_noise(self, f_sky, units='arcmin2'): return self.band_sens**2 * self.get_survey_spread(f_sky, units=units) def get_noise_curves(self, f_sky, ell_max, delta_ell, deconv_beam=True, full_covar=False): ell = np.arange(2, ell_max, delta_ell) W = self.band_sens**2 # Get full covariance; indices are [band,band,ell] ellf = (ell/self.Tatmos_ell[:,None])**(self.Tatmos_alpha[:,None]) T_noise = self.Tatmos_cov[:,:,None] * (ellf[:,None,:] * ellf[None,:,:])**.5 # P noise is tied directly to the white noise level. P_low_noise = (2*W[:,None]) * (ell / self.Patmos_ell[:,None])**self.Patmos_alpha[:,None] P_noise = (self.Patmos_cov[:,:,None] * (P_low_noise[:,None,:] * P_low_noise[None,:,:])**.5) # Add in white noise on the diagonal. for i in range(len(W)): T_noise[i,i] += W[i] P_noise[i,i] += W[i] * 2 # Deconvolve beams. if deconv_beam: beam_sig_rad = self.get_beams() * np.pi/180/60 / (8.*np.log(2))**0.5 beams = np.exp(-0.5 * ell*(ell+1) * beam_sig_rad[:,None]**2) T_noise /= (beams[:,None,:] * beams[None,:,:]) P_noise /= (beams[:,None,:] * beams[None,:,:]) # Diagonal only? if not full_covar: ii = range(self.n_bands) T_noise = T_noise[ii,ii] P_noise = P_noise[ii,ii] sr_per_arcmin2 = (np.pi/180/60)**2 return (ell, T_noise * self.get_survey_spread(f_sky, units='sr'), P_noise * self.get_survey_spread(f_sky, units='sr')) class SOLatV3(SOLatType): def __init__(self, sensitivity_mode=None, N_tubes=None, survey_years=5., survey_efficiency = 0.2*0.85, el=None): # Define the instrument. self.n_bands = 6 self.bands = np.array([ 27., 39., 93., 145., 225., 280.]) self.beams = np.array([ 7.4, 5.1, 2.2, 1.4, 1.0, 0.9]) # Set defaults for survey area, time, efficiency self.survey_years = survey_years self.survey_efficiency = survey_efficiency if sensitivity_mode is None: sensitivity_mode = 1 # baseline self.sensitivity_mode = sensitivity_mode # Sensitivities of each kind of optics tube, in uK rtsec, by # band. 0 represents 0 weight, not 0 noise... nar = np.array self.tube_configs = { 0: { # "threshold" 'LF': nar([ 61., 30., 0, 0, 0, 0 ]), 'MF': nar([ 0, 0, 6.5, 8.1, 0, 0 ])*4**.5, 'UHF': nar([ 0, 0, 0, 0, 17., 42. ])*2**.5, }, 1: { # "baseline" 'LF': nar([ 48., 24., 0, 0, 0, 0 ]), 'MF': nar([ 0, 0, 5.4, 6.7, 0, 0 ])*4**.5, 'UHF': nar([ 0, 0, 0, 0, 15., 10. ])*2**.5, }, 2: { # "goal" 'LF': nar([ 35., 18., 0, 0, 0, 0 ]), 'MF': nar([ 0, 0, 3.9, 4.2, 0, 0 ])*4**.5, 'UHF': nar([ 0, 0, 0, 0, 10., 25. ])*2**.5, }, }[sensitivity_mode] # Save the elevation request. self.el = el # Factor by which to attenuate atmospheric power, given FOV # relative to ACT? self.FOV_mod = 0.5 # Version of the atmospheric model to use? self.atm_version = 0 # The reference tube config. ref_tubes = [('LF', 1), ('MF', 4), ('UHF', 2)] if N_tubes is None: N_tubes = ref_tubes else: N_tubes = [(b,x) for (b,n),x in zip(ref_tubes, N_tubes)] self.precompute(N_tubes) el_data = {'valid': (9., 61.), 'bands': [20,27,39,93,145,225,280], 'polys': [ [1.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00], [1.0865e+00, -4.0827e-03, 6.5708e-05, -3.7321e-07], [1.0433e+00, -2.0940e-03, 3.4543e-05, -1.9931e-07], [1.1805e+00, -8.5858e-03, 1.3929e-04, -7.9535e-07], [1.4914e+00, -2.3383e-02, 3.7955e-04, -2.1686e-06], [2.4305e+00, -6.7918e-02, 1.1002e-03, -6.2787e-06], [2.7908e+00, -8.5208e-02, 1.3833e-03, -7.9075e-06], ] } class S4LatV1(SOLatType): def __init__(self, sensitivity_mode=None, N_tubes=None, N_tels=None, survey_years=7., survey_efficiency=0.23, el=None): # Define the instrument. self.n_bands = 7 self.bands = np.array([ 20., 27., 39., 93., 145., 225., 280.]) self.beams = np.array([ 10., 7.4, 5.1, 2.2, 1.4, 1.0, 0.9]) # Set defaults for survey area, time, efficiency self.survey_years = survey_years self.survey_efficiency = survey_efficiency # Sensitivities of each kind of optics tube, in uK rtsec, by # band. 0 represents 0 weight, not 0 noise... nar = np.array self.tube_configs = { 'ULF': nar([ 60., 0, 0, 0, 0, 0, 0 ]), 'LF': nar([ 0, 28.6, 16.0, 0, 0, 0, 0 ]), 'MF': nar([ 0, 0, 0, 6.6, 6.8, 0, 0 ]), 'UHF': nar([ 0, 0, 0, 0, 0, 12.5, 30.0 ]), } # Save the elevation request. self.el = el # Factor by which to attenuate atmospheric power, given FOV # relative to ACT? self.FOV_mod = 0.5 # Version of the atmospheric model to use? self.atm_version = 1 # The reference tube config. ref_tubes = [('ULF', 1), ('LF', 2), ('MF', 12), ('UHF', 4)] if N_tels is None: N_tels = 2 if N_tubes is None: N_tubes = ref_tubes else: N_tubes = [(b,x) for (b,n),x in zip(ref_tubes, N_tubes)] self.precompute(N_tubes, N_tels) #################################################################### #################################################################### ## demonstration of the code #################################################################### target = 'S4' #target = 'SO' mode=2 suffix='png' fsky=0.4 ellmax=1e4 el=50. colors=['b','r','g','m','k','y'] corr_colors = ['orange', 'fuchsia', 'springgreen'] if target == 'S4': dset_label = 'S4\\_2LAT' lat = S4LatV1(mode, N_tels=2, el=el) corr_pairs = [(1,2),(3,4),(5,6)] ylims = (5e-8,1e-1) colors.insert(0, 'gray') elif target == 'SO': dset_label = 'SO\\_V3' lat = SOLatV3(mode, el=el) corr_pairs = [(0,1),(2,3),(4,5)] ylims = (5e-7,1e0) print(dset_label) bands = lat.get_bands() print("band centers: ", lat.get_bands(), "[GHz]") print("beam sizes: " , lat.get_beams(), "[arcmin]") N_bands = len(bands) ell, N_ell_LA_T_full,N_ell_LA_P_full = lat.get_noise_curves( fsky, ellmax, 1, full_covar=True, deconv_beam=True) WN_levels = lat.get_white_noise(fsky)**.5 N_ell_LA_T = N_ell_LA_T_full[range(N_bands),range(N_bands)] N_ell_LA_Tx = [N_ell_LA_T_full[i,j] for i,j in corr_pairs] N_ell_LA_P = N_ell_LA_P_full[range(N_bands),range(N_bands)] N_ell_LA_Px = [N_ell_LA_P_full[i,j] for i,j in corr_pairs] print("white noise levels: " , WN_levels, "[uK-arcmin]") ## plot the temperature noise curves plt.clf() for i in range(N_bands): plt.loglog(ell,N_ell_LA_T[i], label='%i GHz (%s)' % (bands[i], dset_label), color=colors[i], ls='-', lw=2.) # include correlated atmospheric noise across frequencies for _c,(i,j) in enumerate(corr_pairs): plt.loglog(ell, N_ell_LA_T_full[i,j], label=r'$%i \times %i$ GHz atm.' % (bands[i],bands[j]), color=corr_colors[_c], lw=1.5) plt.title(r"$N(\ell$) Temperature", fontsize=18) plt.ylabel(r"$N(\ell$) [$\mu$K${}^2$]", fontsize=16) plt.xlabel(r"$\ell$", fontsize=16) plt.ylim(*ylims) plt.xlim(100,10000) plt.legend(loc='lower left', ncol=2, fontsize=8) plt.grid() plt.savefig('%s_mode%i_fsky%.2f_LAT_T.%s' % (target, mode, fsky, suffix)) plt.close() ## plot the polarization noise curves plt.clf() for i in range(N_bands): plt.loglog(ell,N_ell_LA_P[i], label='%i GHz (%s)' % (bands[i], dset_label), color=colors[i], ls='-', lw=2.) # include correlated atmospheric noise across frequencies for _c,(i,j) in enumerate(corr_pairs): plt.loglog(ell, N_ell_LA_P_full[i,j], label=r'$%i \times %i$ GHz atm.' % (bands[i],bands[j]), color=corr_colors[_c], lw=1.5) plt.title(r"$N(\ell$) Polarization", fontsize=18) plt.ylabel(r"$N(\ell$) [$\mu$K${}^2$]", fontsize=16) plt.xlabel(r"$\ell$", fontsize=16) plt.ylim(*ylims) plt.xlim(100,10000) plt.legend(loc='upper left', ncol=2, fontsize=9) plt.grid() plt.savefig('%s_mode%i_fsky%.2f_LAT_P.%s' % (target, mode, fsky, suffix)) plt.close() #################################################################### ####################################################################