noise.py 1.52 KB
Newer Older
Maude Le Jeune's avatar
Maude Le Jeune committed
1 2 3
""" noise.py

Generate a noise map from hitcount map.
Maude Le Jeune's avatar
Maude Le Jeune committed
4 5
Compute coupling matrix. 
Compute noise power for debiasing
Maude Le Jeune's avatar
Maude Le Jeune committed
6 7 8 9
"""

import healpy as hp
import pylab as pl
Maude Le Jeune's avatar
Maude Le Jeune committed
10 11
import numpy as npy
import spherelib as sp
Maude Le Jeune's avatar
Maude Le Jeune committed
12 13

### Define some global parameters
Maude Le Jeune's avatar
Maude Le Jeune committed
14 15
lst_par = ['sigma', 'input_hitcount', 'sim_id', 'lmax', 'nside', 'noise_unit', 'noise_power']
lst_tag = lst_par[:-4]
Maude Le Jeune's avatar
Maude Le Jeune committed
16

Maude Le Jeune's avatar
Maude Le Jeune committed
17 18 19
nside  = seg_input.values()[0][0] ## pushed from main
sim_id = seg_input.values()[0][1] ## pushed from main
lmax   = 2*nside
Maude Le Jeune's avatar
Maude Le Jeune committed
20

Maude Le Jeune's avatar
Maude Le Jeune committed
21 22 23 24
### Noise parameters
input_hitcount = "hitmap-wmap5-W.fits"
noise_unit     = "mK_CMB"
sigma          = 6.538
Maude Le Jeune's avatar
Maude Le Jeune committed
25 26 27 28 29 30

### Generate a noise map
hitcount = hp.read_map(input_hitcount)
hitcount = hp.ud_grade(hitcount, nside)
hitnside = pl.sqrt(len(hitcount)/12)
hitcount = hitcount * (hitnside/nside)**2
Maude Le Jeune's avatar
Maude Le Jeune committed
31
noise_map= npy.random.standard_normal(12*nside*nside)*sigma / pl.sqrt(hitcount)
Maude Le Jeune's avatar
Maude Le Jeune committed
32 33 34 35 36 37 38

### Save it to disk
noise_fn = get_data_fn("map_noise.fits")
hp.write_map(noise_fn, noise_map)

### Compute mll'
mll    = sp.mask2mll (hitcount, lmax)
Maude Le Jeune's avatar
Maude Le Jeune committed
39
mll    = mll[0:lmax+1,0:lmax+1]
Maude Le Jeune's avatar
Maude Le Jeune committed
40 41 42
mll_fn = get_data_fn("mll.dat")
save_products (mll_fn, globals(), ["mll"])

Maude Le Jeune's avatar
Maude Le Jeune committed
43 44 45
### Compute noise power
noise_power = 4*pl.pi/len(hitcount)**2 * sigma**2 * pl.sum( 1/hitcount) 
  
Maude Le Jeune's avatar
Maude Le Jeune committed
46 47
### Make a plot
noise_fig = noise_fn.replace('.fits', '.png')
Maude Le Jeune's avatar
Maude Le Jeune committed
48 49 50 51
hp.mollview(noise_map, title="noise in %s"%noise_unit)
pl.savefig (noise_fig)
noise_fig = get_data_fn(input_hitcount).replace('.fits', '.png')
hp.mollview(hitcount, title="hitcount")
Maude Le Jeune's avatar
Maude Le Jeune committed
52 53
pl.savefig (noise_fig)

Maude Le Jeune's avatar
Maude Le Jeune committed
54 55
### Set output : forward as many childs as sim ids
seg_output = [sim_id]