Mismatch Bandpasses

mbp.py

import numpy as np
from planck import LFI
import matplotlib.pyplot as plt

lfi = LFI.LFI()

def horn2freq(horn):
    if horn <= 23:
        return 70
    elif horn <= 26:
        return 44
    else: 
        return 30

# a factors are:
# a = (nu_Side - nu_Main)/(2nu0).
a = np.loadtxt('afacdx7.txt') # released a factors by Adam
a_factors = {}
for row in a:
    a_factors[int(row[0])] = row[1]

# DX9 update
a_factors[27] = .0047
a_factors[28] = -0.0103

# as assumed in the analysis
nu0 = { 30 : 28.45733,
        44 : 44.09855,
        70 : 70.3240 }

def create_tophat(low, high, frequency_axis):
    bandstep = frequency_axis[1] - frequency_axis[0]
    bandwidth_gain = 1/(high-low) * bandstep
    return ( (frequency_axis > low) & (frequency_axis < high) ) * bandwidth_gain

centerfreq = {}
bandwidth = {}

#shift is positive for side arm and negative for main arm
SIGN = {'M':-1, 'S':1} 

neq = np.loadtxt('neq_compression_corrected.txt')

bwtypes =  ['raa bw', '5%-95% bw', 'qucs noise bw', 'neq bw']

print('|* Chan *|* nu_cen *|* nu_min *|* nu_max *|')
for horn in range(18, 28+1):
    freq = horn2freq(horn)

    shift = a_factors[horn] * 2 * nu0[freq] # nu_Side - nu_Main

    for side in ['M','S']:
        
        chtag = 'LFI%d%s' % (horn, side)
        bandwidth[chtag] = neq[horn-18]
        centerfreq[chtag] = nu0[freq] + SIGN[side] * shift/2
        print('| %s | %.3f | %.3f | %.3f |' % (chtag, centerfreq[chtag], centerfreq[chtag]-bandwidth[chtag]/2,  centerfreq[chtag]+bandwidth[chtag]/  2))

        #plt.figure()
        #plt.plot(qucs_bp['WAVENUMBER'], qucs_bp['TRANSMISSION'], label='QUCS')
        #for t in ['5%-95% bw', 'qucs noise bw', 'neq bw']:
        #    plt.plot(qucs_bp['WAVENUMBER'], create_tophat(nu0[ch.f.freq]-bandpass[t][chtag]/2., nu0[ch.f.freq]+bandpass[t][chtag]/2., qucs_bp['WAVENUMBER']), label=t)
        #plt.plot(qucs_bp['WAVENUMBER'], create_tophat(cross_low, cross_high, qucs_bp['WAVENUMBER']), label='cut')
        #plt.legend(); plt.grid(); plt.xlabel('Frequency'); plt.title(ch.tag)
        #plt.vlines(center_freq, np.min(qucs_bp['TRANSMISSION']), np.max(qucs_bp['TRANSMISSION']))

frequency_axis = np.arange(25, 80, .01)
for horn in range(18, 28+1):
    for side in ['M','S']:
        chtag = 'LFI%d%s' % (horn, side)
        plt.plot(frequency_axis, create_tophat(centerfreq[chtag]-bandwidth[chtag]/  2,  centerfreq[chtag]+bandwidth[chtag]/ 2, frequency_axis), label=chtag)
plt.xlabel('Freq [GHz]')

#print('Bandpasses GHz')
#print('|* Ch*|* ' + '*|*'.join(bandpass.keys()) + '*|')
#for k in pl.ch:
#    print(('|%s| ' % k.tag) + (' | '.join([('%.2f' % bandpass[t][k.tag]) for t in bandpass.keys()])) + '|')
#
#hornax = np.arange(18, 28+1, .5)
#
#isM = {True:'M', False:'S'}
#for t in bandpass.keys():
#    plt.plot(hornax, [bandpass[t]['LFI%d%s' % (int(ha), isM[int(ha)==ha])] for ha in hornax], label=t)
#ax = plt.gca()
##ax.xaxis.set_major_formatter(
#plt.xlabel('Radiometers')
#plt.xlim([17.5, 28.5])
#plt.legend(loc=0)
#plt.title("LFI Bandwidths")
#plt.grid()
#

def qucs_band(chtag):
    band = np.loadtxt("radiometer_bandpasses/bandpass_%s.dat" % chtag)
    return np.average(band[:,0], weights=band[:,1])

print('|* Chan *|* nu_cen(a_fac) *|* nu_cen(idb) *|* nu_cen(qucs) *|* nu_cen(mix) *|')
for horn in range(18, 28+1):
    freq = horn2freq(horn)

    shift = a_factors[horn] * 2 * nu0[freq] # nu_Side - nu_Main

    qucs_horn_center = (lfi["LFI%dM"%horn].get_instrument_db_field("nu_cen") + lfi["LFI%dS"%horn].get_instrument_db_field("nu_cen"))/2.e9
    #qucs_horn_center = (qucs_band("LFI%dM"%horn)+qucs_band("LFI%dS"%horn))/2.
    for side in ['M','S']:
        
        chtag = 'LFI%d%s' % (horn, side)
        ch = lfi[chtag]
        centerfreq[chtag] = nu0[freq] + SIGN[side] * shift/2
        print('| %s | %.2f | %.2f | %.2f | %.2f |' % (chtag, centerfreq[chtag], ch.get_instrument_db_field("nu_cen")/1.e9, qucs_band(chtag), qucs_horn_center  + SIGN[side] * shift/2
))
print('|* horn *|* nu_cen(qucs)-ref *|')
for horn in range(18, 28+1):
    freq = horn2freq(horn)

    shift = a_factors[horn] * 2 * nu0[freq] # nu_Side - nu_Main

    qucs_horn_center = (lfi["LFI%dM"%horn].get_instrument_db_field("nu_cen") + lfi["LFI%dS"%horn].get_instrument_db_field("nu_cen"))/2.e9


    print('| %d | %.2f |' % (horn, qucs_horn_center - nu0[freq]))


print('|* horn *|* DX9 a-factors *|* QUCS a-factors *|')
a_fac_qucs = {}
for horn in range(18, 28+1):
    freq = horn2freq(horn)
    qucs_horn_center = (qucs_band("LFI%dM"%horn)+qucs_band("LFI%dS"%horn))/2.
    a_fac_qucs[horn] = (qucs_band("LFI%dS"%horn) - qucs_band("LFI%dM"%horn))/(2. * qucs_horn_center)
    #a_fac_qucs[horn] = (centerfreq["LFI%dS"%horn] - centerfreq["LFI%dM"%horn])/(2. * nu0[freq])
    print('| %d | %.2e | %.2e |' % (horn, a_factors[horn], a_fac_qucs[horn]))