# tophat_bandpass.py # ------------------ # Calculate equivalent tophat bandpasses for the bandpass files used in # Victor's CMB-S4 optimization and forecasting. # # Bandpass files can be found in Table 1 of this posting: # http://users.physics.harvard.edu/~buza/20160531_fisher/ # # Colin Bischoff (colin.bischoff@uc.edu) 2016-11-01 import numpy as np import matplotlib matplotlib.use('agg') import matplotlib.pyplot as plt import os # Physical constants. h = 6.626e-34 # m^2 kg s^{-1} k = 1.381e-23 # m^2 kg s^{-2} K^{-1} hk = h / k * 1.e9 Tcmb = 2.72548 # K def dust_scaling(nu, bandpass, beta=1.59, Tdust=19.6, nu0=353.0): """ Calculate the relative brightness of a dust signal observed with the specified bandpass, compared to a reference frequency. Params ------ nu : (N,) ndarray Array of frequency values, in GHz, used for bandpass definition. bandpass : (N,) ndarray Bandpass for which to calculate dust scale factor. The bandpass convention is defined relative to the response of a single-moded antenna to a beam-filling source. Normalization of the bandpass doesn't matter. beta : float, optional Dust emissivity parameter. Default value is 1.59. Tdust : float, optional Dust temperature in Kelvin. Default value is 19.6 K. nu0 : float, optional Reference frequency for dust brightness calculation. Default is 353 GHz. """ # Calculate greybody scaling. gb = (bandpass * nu**(1.0 + beta) / (np.exp(hk * nu / Tdust) - 1.0)).sum() gb0 = nu0**(1.0 + beta) / (np.exp(hk * nu0 / Tdust) - 1.0) # Calculate Rayleigh-Jeans to thermodynamic temperature conversion. th = (bandpass * nu**2 / (np.exp(hk * nu / Tcmb) - 1.0)).sum() th0 = nu0**2 / (np.exp(hk * nu0 / Tcmb) - 1.0) # Ratio of greybody scaling divided by ratio of thermodynamic conversion. # Normalization factors from integrals over bandpass cancel out. return (gb / gb0) / (th / th0) def sync_scaling(nu, bandpass, beta=-3.0, nu0=23.0): """ Calculate the relative brightness of a synchrotron signal observed with the specified bandpass, compared to a reference frequency. Params ------ nu : (N,) ndarray Array of frequency values, in GHz, used for bandpass definition. bandpass : (N,) ndarray Bandpass for which to calculate dust scale factor. The bandpass convention is defined relative to the response of a single-moded antenna to a beam-filling source. Normalization of the bandpass doesn't matter. beta : float, optional Synchrotron power-law parameter. Default value is -3.0. nu0 : float, optional Reference frequency for dust brightness calculation. Default is 23 GHz. """ # Calculate power-law scaling. pl = (bandpass * nu**beta).sum() pl0 = nu0**beta # Calculate Rayleigh-Jeans to thermodynamic temperature conversion. th = (bandpass * nu**2 / (np.exp(hk * nu / Tcmb) - 1.0)).sum() th0 = nu0**2 / (np.exp(hk * nu0 / Tcmb) - 1.0) # Ratio of power-law scaling divided by ratio of thermodynamic conversion. # Normalization factors from integrals over bandpass cancel out. return (pl / pl0) / (th / th0) # Atmospheric spectra are courtesy of Denis Barkats, generated using am: # https://www.cfa.harvard.edu/~spaine/am/index.html # Files have four columns: # Column 0 = frequency, in GHz # Column 1 = opacity, in nepers # Column 2 = Planck brightness temperature, in Kelvin # Column 3 = transmittance spole = np.genfromtxt('spole_winter_0p5mm.out') spole = spole[1:,:] # Drop nu=0 point. atacama = np.genfromtxt('chajnantor_1mm.out') atacama = atacama[1:,:] # Drop nu=0 point. nu = spole[:,0] # Same for both files. dnu = nu[1] - nu[0] # Band centers and widths from Victor's table. # Increase 215->220 GHz for band 7 center. # Increase fractional band with for band 8 from 18% to 22%. nband = 8 center = np.array([30.0, 40.0, 85.0, 95.0, 145.0, 155.0, 220.0, 270.0]) width = center * np.array([0.3, 0.3, 0.24, 0.24, 0.22, 0.22, 0.22, 0.22]) bandpass = np.zeros((len(nu), nband)) scl_dust = np.zeros((nband,)) scl_sync = np.zeros((nband,)) tsky_pole = np.zeros((nband,)) tsky_atacama = np.zeros((nband,)) for ii in range(nband): # Tophat bandpass. nu1 = center[ii] - width[ii] / 2.0 nu2 = center[ii] + width[ii] / 2.0 bandpass[:,ii] = np.zeros(nu.shape) bandpass[(nu >= nu1) & (nu < nu2),ii] = 1.0 # Calculate dust and sync scale factors. scl_dust[ii] = dust_scaling(nu, bandpass[:,ii]) scl_sync[ii] = sync_scaling(nu, bandpass[:,ii]) # Calculate band-averaged brightness temperature for South Pole and # Atacama. tsky_pole[ii] = (bandpass[:,ii] * spole[:,2]).sum() / bandpass[:,ii].sum() tsky_atacama[ii] = (bandpass[:,ii] * atacama[:,2]).sum() / bandpass[:,ii].sum() # Also get the actual BICEP2/Keck spectra. BK_files = ('K95_frequency_spectrum_20150309.txt', 'B2_frequency_spectrum_20141216.txt', 'K220_frequency_spectrum_20160120.txt') nband_bk = 3 bandpass_bk = np.zeros((len(nu), nband_bk)) center_bk = np.zeros((nband_bk,)) width_bk = np.zeros((nband_bk,)) scl_dust_bk = np.zeros((nband_bk,)) scl_sync_bk = np.zeros((nband_bk,)) tsky_pole_bk = np.zeros((nband_bk,)) tsky_atacama_bk = np.zeros((nband_bk,)) for ii in range(nband_bk): bkdata = np.genfromtxt(BK_files[ii], delimiter=',') if ii == 2: # For 220 GHz bandpass, need to manually calculate third column. bkdata[:,2] = bkdata[:,1] / bkdata[:,0]**2 bkdata[:,2] = bkdata[:,2] / np.max(bkdata[:,2]) # Interpolate to fine-binned frequency array. bandpass_bk[:,ii] = np.interp(nu, bkdata[:,0], bkdata[:,2]) bandpass_bk[nu < bkdata[0,0], ii] = 0. bandpass_bk[nu > bkdata[-1,0], ii] = 0. # Calculate band center and width. center_bk[ii] = (nu * bandpass_bk[:,ii]).sum() / bandpass_bk[:,ii].sum() width_bk[ii] = (bandpass_bk[:,ii].sum())**2 / (bandpass_bk[:,ii]**2).sum() * dnu # Calculate foreground scale factors. scl_dust_bk[ii] = dust_scaling(nu, bandpass_bk[:,ii]) scl_sync_bk[ii] = sync_scaling(nu, bandpass_bk[:,ii]) # Calculate band-averaged sky brightness. tsky_pole_bk[ii] = (bandpass_bk[:,ii] * spole[:,2]).sum() / bandpass_bk[:,ii].sum() tsky_atacama_bk[ii] = (bandpass_bk[:,ii] * atacama[:,2]).sum() / bandpass_bk[:,ii].sum() # Make figure showing bandpasses on top of atmospheric spectra. fig = plt.figure(figsize=(8, 4)) plt.plot(spole[:,0], spole[:,2], 'k', label='South Pole, 0.5mm') plt.plot(atacama[:,0], atacama[:,2], color='grey', label='Atacama, 1mm') for ii in range(nband): if np.mod(ii, 2) == 0: clr = 'red' else: clr = 'blue' plt.plot(nu, bandpass[:,ii] * tsky_pole[ii], color=clr, linewidth=2) for ii in range(nband_bk): # To get same normalization as the tophats, scale BK bandpass so that the integral # is equal to tsky * bandwidth. plt.plot(nu, (bandpass_bk[:,ii] * tsky_pole_bk[ii] * width_bk[ii] / bandpass_bk[:,ii].sum() / dnu), color='green') plt.xlim([0, 325]) plt.ylim([0, 50]) plt.legend() plt.xlabel('frequency [GHz]') plt.ylabel('brightness temperature [K]') plt.tight_layout() fig.savefig('tophat_bandpass.png') plt.close(fig) # Print results to text file. fid = open('tophat_bandpass.txt', 'w') for ii in range(nband): fid.write('{:6.2f} {:5.2f} {:6.4f} {:6.4f} {:4.1f} {:4.1f}\n'. format(center[ii], width[ii], scl_dust[ii], scl_sync[ii], tsky_pole[ii], tsky_atacama[ii])) for ii in range(nband_bk): fid.write('{:6.2f} {:5.2f} {:6.4f} {:6.4f} {:4.1f} {:4.1f}\n'. format(center_bk[ii], width_bk[ii], scl_dust_bk[ii], scl_sync_bk[ii], tsky_pole_bk[ii], tsky_atacama_bk[ii])) fid.close()