Python QubicAcquisition.get_projection_peak_operator方法代碼示例

# 需要導入模塊: from qubic import QubicAcquisition [as 別名]
# 或者: from qubic.QubicAcquisition import get_projection_peak_operator [as 別名]
def tod2map(tod,pointing,instrument_in,detector_list=False,disp=True,kmax=2,displaytime=False):
    t0=time.time()
    #### Detectors
    mask_packed = np.ones(len(instrument_in.detector.packed), bool)
    if detector_list:
        mask_packed[detector_list] = False
        mask_unpacked = instrument_in.unpack(mask_packed)
        instrument = QubicInstrument('monochromatic', removed=mask_unpacked,nside=instrument_in.sky.nside)
    else:
        instrument = instrument_in
        
    #### Observations
    obs = QubicAcquisition(instrument, pointing)
    projection = obs.get_projection_peak_operator(kmax=kmax)
    coverage = projection.pT1()
    mask = coverage == 0
    projection = pack_projection_inplace(projection, mask)
    hwp = obs.get_hwp_operator()
    polgrid = DenseOperator([[0.5, 0.5, 0],
                             [0.5,-0.5, 0]])
    H = polgrid * hwp * projection
    preconditioner = DiagonalOperator(1/coverage[~mask], broadcast='rightward')
    solution = pcg(H.T * H, H.T(tod), M=preconditioner, disp=disp, tol=1e-3)
    output_map = unpack(solution['x'], mask)
    t1=time.time()
    if displaytime: print(' Map done in {0:.4f} seconds'.format(t1-t0))
    return output_map,coverage

# 需要導入模塊: from qubic import QubicAcquisition [as 別名]
# 或者: from qubic.QubicAcquisition import get_projection_peak_operator [as 別名]
def map2tod(maps,pointing,instrument_in,detector_list=False,kmax=2):
    #### Detectors
    mask_packed = np.ones(len(instrument_in.detector.packed), bool)
    if detector_list:
        mask_packed[detector_list] = False
        mask_unpacked = instrument_in.unpack(mask_packed)
        instrument = QubicInstrument('monochromatic', removed=mask_unpacked,nside=instrument_in.sky.nside)
    else:
        instrument = instrument_in
        
    #### Observations
    obs = QubicAcquisition(instrument, pointing)
    #C = obs.get_convolution_peak_operator()
    #convmaps=np.transpose(np.array([C(maps[:,0]),C(maps[:,1]),C(maps[:,2])]))
    projection = obs.get_projection_peak_operator(kmax=kmax)
    coverage = projection.pT1()
    mask = coverage == 0
    projection = pack_projection_inplace(projection, mask)
    hwp = obs.get_hwp_operator()
    polgrid = DenseOperator([[0.5, 0.5, 0],
                             [0.5,-0.5, 0]])
    #H = polgrid * hwp * projection * C
    H = polgrid * hwp * projection
    x1 = pack(maps, mask)
    y = H(x1)
    #return y,convmaps
    return y

# 需要導入模塊: from qubic import QubicAcquisition [as 別名]
# 或者: from qubic.QubicAcquisition import get_projection_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_projection_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_projection_peak_operator [as 別名]
def TOD2map(tod,pointing,nside,kmax=2,disp=True,P=False,covmin=10):
    qubic = QubicInstrument('monochromatic,nopol',nside=nside)
    #### configure observation
    obs = QubicAcquisition(qubic, pointing)
    if not P:
        P = obs.get_projection_peak_operator(kmax=kmax)
    coverage = P.T(np.ones_like(tod))
    mask=coverage < covmin
    P.matrix.pack(mask)
    P_packed = ProjectionInMemoryOperator(P.matrix)
    unpack = UnpackOperator(mask)
    # data
    solution = pcg(P_packed.T * P_packed, P_packed.T(tod), M=DiagonalOperator(1/coverage[~mask]), disp=disp)
    output_map = unpack(solution['x'])
    output_map[mask] = np.nan
    coverage[mask]=np.nan
    return(output_map,mask,coverage)

# 需要導入模塊: from qubic import QubicAcquisition [as 別名]
# 或者: from qubic.QubicAcquisition import get_projection_peak_operator [as 別名]
def test_load_save():
    info = 'test\nconfig'

    def func_acquisition(obs, info):
        filename_ = os.path.join(outpath, 'config-' + str(uuid1()))
        obs.save(filename_, info)
        obs2, info2 = QubicAcquisition.load(filename_)
        assert_equal(str(obs), str(obs2))
        assert_equal(info, info2)

    def func_observation(obs, tod, info):
        filename_ = os.path.join(outpath, 'obs-' + str(uuid1()))
        obs._save_observation(filename_, tod, info)
        obs2, tod2, info2 = QubicAcquisition._load_observation(filename_)
        assert_equal(str(obs), str(obs2))
        assert_equal(tod, tod2)
        assert_equal(info, info2)

    def func_simulation(obs, input_map, tod, info):
        filename_ = os.path.join(outpath, 'simul-' + str(uuid1()))
        obs.save_simulation(filename_, input_map, tod, info)
        obs2, input_map2, tod2, info2 = QubicAcquisition.load_simulation(
            filename_)
        assert_equal(str(obs), str(obs2))
        assert_equal(input_map, input_map2)
        assert_equal(tod, tod2)
        assert_equal(info, info2)

    for p in ptgs:
        p = np.asarray(p)
        smp = QubicSampling(azimuth=p[..., 0], elevation=p[..., 1],
                            pitch=p[..., 2])
        acq = QubicAcquisition(qubic, smp)
        P = acq.get_projection_peak_operator()
        tod = P(input_map)
        yield func_acquisition, acq, info
        yield func_observation, acq, tod, info
        yield func_simulation, acq, input_map, tod, info

# 需要導入模塊: from qubic import QubicAcquisition [as 別名]
# 或者: from qubic.QubicAcquisition import get_projection_peak_operator [as 別名]
    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)
e2g = SphericalEquatorial2GalacticOperator(degrees=True)