Python QubicAcquisition.get_convolution_peak_operator方法代码示例

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_convolution_peak_operator [as 别名]
def get_qubic_map(instrument, sampling, scene, input_maps, withplanck=True, covlim=0.1):
    acq = QubicAcquisition(instrument, sampling, scene, photon_noise=True, effective_duration=1)
    C = acq.get_convolution_peak_operator()
    coverage = acq.get_coverage()
    observed = coverage > covlim * np.max(coverage)
    acq_restricted = acq[:, :, observed]
    H = acq_restricted.get_operator()
    x0_convolved = C(input_maps)
    if not withplanck:
        pack = PackOperator(observed, broadcast='rightward')
        y_noiseless = H(pack(x0_convolved))
        noise = acq.get_noise()
        y = y_noiseless + noise
        invntt = acq.get_invntt_operator()
        A = H.T * invntt * H
        b = (H.T * invntt)(y)
        preconditioner = DiagonalOperator(1 / coverage[observed], broadcast='rightward')
        solution_qubic = pcg(A, b, M=preconditioner, disp=True, tol=1e-3, maxiter=1000)
        maps = pack.T(solution_qubic['x'])
        maps[~observed] = 0
    else:
        acq_planck = PlanckAcquisition(150, acq.scene, true_sky=x0_convolved)#, fix_seed=True)
        acq_fusion = QubicPlanckAcquisition(acq, acq_planck)
        map_planck_obs=acq_planck.get_observation()
        H = acq_fusion.get_operator()
        invntt = acq_fusion.get_invntt_operator()
        y = acq_fusion.get_observation()
        A = H.T * invntt * H
        b = H.T * invntt * y
        solution_fusion = pcg(A, b, disp=True, maxiter=1000, tol=1e-3)
        maps = solution_fusion['x']
        maps[~observed] = 0
    x0_convolved[~observed] = 0
    return(maps, x0_convolved, observed)    

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_convolution_peak_operator [as 别名]
def get_coverage_onedet(idetector=231, angspeed=1., delta_az=15., angspeed_psi=0., maxpsi=15., nsweeps_el=100, duration=24, ts=1., decrange=2., decspeed=2., recenter=False):
    print('##### Getting Coverage for: ') 
    print('## idetector = '+str(idetector))
    print('## duration = '+str(duration))
    print('## ts = '+str(ts))
    print('## recenter = '+str(recenter))
    print('## angspeed = '+str(angspeed))
    print('## delta_az = '+str(delta_az))
    print('## nsweeps_el = '+str(nsweeps_el))
    print('## decrange = '+str(decrange))
    print('## decspeed = '+str(decspeed))
    print('## angspeed_psi = '+str(angspeed_psi))
    print('## maxpsi = '+str(maxpsi))
    print('##########################')
    nside = 256
    racenter = 0.0
    deccenter = -57.0
    pointing = pointings_modbyJC.create_sweeping_pointings(
        [racenter, deccenter], duration, ts, angspeed, delta_az, nsweeps_el,
        angspeed_psi, maxpsi, decrange=decrange, decspeed=decspeed, recenter=recenter)
    pointing.angle_hwp = np.random.random_integers(0, 7, pointing.size) * 22.5
    ntimes = len(pointing)

    # get instrument model with only one detector
    instrument = QubicInstrument('monochromatic')
    mask_packed = np.ones(len(instrument.detector.packed), bool)
    mask_packed[idetector] = False
    mask_unpacked = instrument.unpack(mask_packed)
    instrument = QubicInstrument('monochromatic', removed=mask_unpacked)
    
    obs = QubicAcquisition(instrument, pointing)
    convolution = obs.get_convolution_peak_operator()
    projection = obs.get_projection_peak_operator(kmax=0)
    coverage = projection.pT1()
    return coverage

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_convolution_peak_operator [as 别名]
def map2TOD(input_map, pointing,kmax=2):
    ns=hp.npix2nside(input_map.size)
    qubic = QubicInstrument('monochromatic,nopol',nside=ns)
    #### configure observation
    obs = QubicAcquisition(qubic, pointing)
    C = obs.get_convolution_peak_operator()
    P = obs.get_projection_peak_operator(kmax=kmax)
    H = P * C
    # Produce the Time-Ordered data
    tod = H(input_map)
    input_map_conv = C(input_map)
    return(tod,input_map_conv,P)

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_convolution_peak_operator [as 别名]
def get_tod(instrument, sampling, scene, input_maps, withplanck=True, covlim=0.1, photon_noise=True): 
    acq = QubicAcquisition(instrument, sampling, scene, photon_noise=photon_noise)
    C = acq.get_convolution_peak_operator()
    coverage = acq.get_coverage()
    observed = coverage > covlim * np.max(coverage)
    acq_restricted = acq[:, :, observed]
    H = acq_restricted.get_operator()
    x0_convolved = C(input_maps)
    pack = PackOperator(observed, broadcast='rightward')
    y_noiseless = H(pack(x0_convolved))
    noise = acq.get_noise()
    y = y_noiseless + noise
    return (y_noiseless, noise, y)

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_convolution_peak_operator [as 别名]
angspeed = 1  # deg/sec
delta_az = 15.
angspeed_psi = 0.1
maxpsi = 45.
nsweeps_el = 300
duration = 24   # hours
ts = 20         # seconds
sampling = create_sweeping_pointings(
    [racenter, deccenter], duration, ts, angspeed, delta_az, nsweeps_el,
    angspeed_psi, maxpsi)

# acquisition model
acq = QubicAcquisition(150, sampling, kind='I', synthbeam_fraction=0.99,
                       detector_sigma=sigma, detector_fknee=fknee,
                       detector_fslope=fslope, detector_ncorr=ncorr)
C = acq.get_convolution_peak_operator()
P = acq.get_projection_operator()
H = P * C

# produce the Time-Ordered data
y = H(x0)

# noise
psd = _gaussian_psd_1f(len(acq.sampling), sigma=sigma, fknee=fknee,
                       fslope=fslope, sampling_frequency=1/ts)
invntt = acq.get_invntt_operator()
noise = acq.get_noise()

# map-making
coverage = P.pT1()
mask = coverage > 10

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_convolution_peak_operator [as 别名]
center_gal = 0, 90
center = gal2equ(center_gal[0], center_gal[1])

# sampling model
np.random.seed(0)
sampling = create_random_pointings(center, 1000, 10)

# scene model
scene = QubicScene(hp.npix2nside(x0.size), kind='I')

# instrument model
instrument = QubicInstrument(filter_nu=150e9)

# acquisition model
acq = QubicAcquisition(instrument, sampling, scene)
x0_convolved = acq.get_convolution_peak_operator()(x0)
H = acq.get_operator()
coverage = H.T(np.ones(H.shapeout))
mask = coverage > 0

# restrict the scene to the observed pixels
acq_restricted = acq[..., mask]
H_restricted = acq_restricted.get_operator()
x0_restricted = x0[mask]
y = H_restricted(x0_restricted)
invntt = acq_restricted.get_invntt_operator()

# solve for x
A = H_restricted.T * invntt * H_restricted
b = H_restricted.T(invntt(y))
solution = pcg(

# 需要导入模块: from qubic import QubicAcquisition [as 别名]
# 或者: from qubic.QubicAcquisition import get_convolution_peak_operator [as 别名]
    [racenter, deccenter], duration, ts, angspeed, delta_az, nsweeps_el,
    angspeed_psi, maxpsi, decrange=decrange, decspeed=decspeed,recenter=True)
pointing.angle_hwp = np.random.random_integers(0, 7, pointing.size) * 22.5
ntimes = len(pointing)

# get instrument model with only one detector
idetector = 231
instrument = QubicInstrument('monochromatic')
instrument.plot()
mask_packed = np.ones(len(instrument.detector.packed), bool)
mask_packed[idetector] = False
mask_unpacked = instrument.unpack(mask_packed)
instrument = QubicInstrument('monochromatic', removed=mask_unpacked)

obs = QubicAcquisition(instrument, pointing)
convolution = obs.get_convolution_peak_operator()
projection = obs.get_projection_peak_operator(kmax=0)
hwp = obs.get_hwp_operator()
polarizer = obs.get_polarizer_operator()
coverage = projection.pT1()
mask = coverage > 0
projection.restrict(mask)
pack = PackOperator(mask, broadcast='rightward')
ra,dec=pointings_modbyJC._hor2equ(pointing.azimuth, pointing.elevation, pointing.latitude, pointing.time/3600)
ns_scan = int(delta_az / angspeed / 0.1)
ns_tot = ns_scan * 2
ichunk = (np.arange(len(pointing)) / ns_tot / nsweeps_el).astype(numpy.int64)
isweep = ((np.arange(len(pointing)) / ns_tot).astype(numpy.int64)) % nsweeps_el
clf()
subplot(2,2,1)
plot(pointing.azimuth,pointing.elevation)