Python QubicAcquisition.get_coverage方法代碼示例

# 需要導入模塊: from qubic import QubicAcquisition [as 別名]
# 或者: from qubic.QubicAcquisition import get_coverage [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_coverage [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_coverage [as 別名]
C = acq.get_convolution_peak_operator()
H = H_ga * C

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

# noise
sigma = acq.instrument.detector.nep / np.sqrt(2 * sampling.period)
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()
noise[...] = 0

# map-making
coverage = acq.get_coverage()
mask = coverage / coverage.max() > 0.01

acq_red = acq[..., mask]
H_ga_red = acq_red.get_operator()
# map without covariance matrix
solution1 = pcg(H_ga_red.T * H_ga_red, H_ga_red.T(y + noise),
                M=DiagonalOperator(1/coverage[mask]), disp=True)
x1 = acq_red.scene.unpack(solution1['x'])

# map with covariance matrix
solution2 = pcg(H_ga_red.T * invntt * H_ga_red,
                (H_ga_red.T * invntt)(y + noise),
                M=DiagonalOperator(1/coverage[mask]), disp=True)
x2 = acq_red.scene.unpack(solution2['x'])

# 需要導入模塊: from qubic import QubicAcquisition [as 別名]
# 或者: from qubic.QubicAcquisition import get_coverage [as 別名]
sky = read_map(PATH + 'syn256_pol.fits')
sampling = create_random_pointings([racenter, deccenter], 1000, 10)


all_solutions_fusion = []
all_coverages = []

nbptg = np.linspace(1000,5000,5)
correct_time = 365*86400./(nbptg/1000)
detector_nep = 4.7e-17/np.sqrt(correct_time / len(sampling)*sampling.period)

for i in xrange(len(all_instruments)):
	acq_qubic = QubicAcquisition(150, sampling, nside=nside,
                             detector_nep=detector_nep[i])
	all_coverages.append(acq_qubic.get_coverage())
	convolved_sky = acq_qubic.instrument.get_convolution_peak_operator()(sky)
	acq_planck = PlanckAcquisition(150, acq_qubic.scene, true_sky=convolved_sky)
	acq_fusion = QubicPlanckAcquisition(acq_qubic, acq_planck)

	H = acq_fusion.get_operator()
	invntt = acq_fusion.get_invntt_operator()
	obs = acq_fusion.get_observation()

	A = H.T * invntt * H
	b = H.T * invntt * obs

	solution_fusion = pcg(A, b, disp=True)
	all_solutions_fusion.append(solution_fusion)

# 需要導入模塊: from qubic import QubicAcquisition [as 別名]
# 或者: from qubic.QubicAcquisition import get_coverage [as 別名]
sampling = create_random_pointings([racenter, deccenter], 5000, 5)
scene = QubicScene(nside, kind='I')

#sky = read_map(PATH + 'syn256_pol.fits')[:,0]
sky = np.zeros(12*nside**2)
ip0=hp.ang2pix(nside, np.radians(90.-center[1]), np.radians(center[0]))
v0 = np.array(hp.pix2vec(nside, ip0))
sky[ip0] = 1000
#hp.gnomview(sky, rot=center, reso=5, xsize=800)

instrument = QubicInstrument(filter_nu=150e9,
                            detector_nep=1e-30)

acq_qubic = QubicAcquisition(instrument, sampling, scene, effective_duration=1, photon_noise=False)

coverage = acq_qubic.get_coverage()
observed = coverage > 0.01 * np.max(coverage)
H = acq_qubic.get_operator()
invntt = acq_qubic.get_invntt_operator()
y, sky_convolved = acq_qubic.get_observation(sky, convolution=True)

A = H.T * invntt * H
b = H.T * invntt * y

tol = 5e-6
solution_qubic = pcg(A, b, disp=True, maxiter=maxiter, tol=tol)



################ With only one detector
sampling = create_random_pointings([racenter, deccenter], 500000, 15)

# 需要導入模塊: from qubic import QubicAcquisition [as 別名]
# 或者: from qubic.QubicAcquisition import get_coverage [as 別名]
import matplotlib.pyplot as mp
import numpy as np

racenter = 0.0      # deg
deccenter = -57.0   # deg
maxiter = 1000
tol = 5e-6

sky = read_map(PATH + 'syn256_pol.fits')
nside = 256
sampling = create_random_pointings([racenter, deccenter], 1000, 10)
scene = QubicScene(nside)

acq_qubic = QubicAcquisition(150, sampling, scene, effective_duration=1)
convolved_sky = acq_qubic.instrument.get_convolution_peak_operator()(sky)
cov = acq_qubic.get_coverage()
mask = cov < cov.max() * 0.2
acq_planck = PlanckAcquisition(150, acq_qubic.scene, true_sky=convolved_sky, mask=mask)
acq_fusion = QubicPlanckAcquisition(acq_qubic, acq_planck)

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=maxiter, tol=tol)

acq_qubic = QubicAcquisition(150, sampling, scene, effective_duration=1)
H = acq_qubic.get_operator()

# 需要導入模塊: from qubic import QubicAcquisition [as 別名]
# 或者: from qubic.QubicAcquisition import get_coverage [as 別名]
savefig('azel.png')


ts=5.
nsweeps_el = 300
delta_az = 40.
sampling = create_sweeping_pointings(
        [racenter, deccenter], duration, ts, angspeed, delta_az, nsweeps_el,
        angspeed_psi, maxpsi)


detector_nep = 4.7e-17*np.sqrt(len(sampling) * sampling.period / (365 * 86400))
acq_qubic = QubicAcquisition(150, sampling[np.abs(sampling.elevation-50) < 20], nside=nside,
                             detector_nep=detector_nep)
                            
coverage_map = acq_qubic.get_coverage()
coverage_map = coverage_map / np.max(coverage_map)
angmax = hp.pix2ang(nside, coverage_map.argmax())
maxloc = np.array([np.degrees(angmax[1]), 90.-np.degrees(angmax[0])])

figure(0)
clf()
cov = hp.gnomview(coverage_map, rot=maxloc, reso=5, xsize=800, return_projected_map=True,sub=(2,1,1),min=0.0)
contour(cov,[0,0.01, 0.1])
subplot(2,1,2)
x, y = profile(cov)
x *= 5 / 60
plot(x, y/np.max(y))

maskok = coverage_map > 0.1
fsky=np.sum(coverage_map[maskok]) *1./len(maskok)