In [1]:
## This notebook shows how to simulate observations from a CMB ground-based experiment.
## The flow is the following: setup your env -> generate scans -> covariance maps -> noise maps

In [2]:
import ConfigParser
import os
import sys

## No need for installation - use directly the module
## Thx Antony!
sys.path.insert(0,os.path.realpath(os.path.join(os.getcwd(),'..')))

from LaFabrique import communication as comm
from LaFabrique import util_CMB
from LaFabrique import scanning_strategy
from LaFabrique import noise

import healpy as hp
import matplotlib.pyplot as pl


Parallel setup OK, rank 0 in 1

In [3]:
#### SET UP YOUR ENVIRONMENT
## Load the environment (paths, etc)
Config = ConfigParser.ConfigParser()
Config.read('../test/setup_env_test.ini')
environment = util_CMB.normalise_env_parser(
        Config._sections['Environment'])

## Initialise paths
environment.outpath_noise = os.path.join(
    environment.out_path, 'noise')
environment.outpath_masks = os.path.join(
    environment.out_path, 'masks')
environment.outpath_foregrounds = os.path.join(
    environment.out_path, 'foregrounds')

## Create folders if necessary
if comm.rank == 0:
    ## Create root folder
    if not os.path.exists(environment.out_path):
        os.makedirs(environment.out_path)

    ## Create folders for noise and masks
    if not os.path.exists(environment.outpath_noise):
        os.makedirs(environment.outpath_noise)
    if not os.path.exists(environment.outpath_masks):
        os.makedirs(environment.outpath_masks)
    if not os.path.exists(environment.outpath_foregrounds):
        os.makedirs(environment.outpath_foregrounds)

In [6]:
#### SIMULATE A SCAN
## Keep calm and relaunch twice if you get errors (from the compilation of C. Weave is quite picky...)
scanning_strategy.generate_scans('../test/setup_scanning_test.ini', environment)

## You should get your generic scan here (hdf5 file)
!ls test/test/masks


setup_scanning.ini test.hdf5

In [7]:
#### SIMULATE COVARIANCE AND NOISE MAPS AT ALL DESIRED FREQUENCIES FROM THE PREVIOUS SCAN

## generate_noise_sims() will simulate noise maps and covariance matrix.
## In the instrument ini file, you can specify to simulate only the covariance matrix.
noise.generate_noise_sims('../test/setup_instrument_test.ini', comm=comm, env=environment)


/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/healpy/pixelfunc.py:306: RuntimeWarning: invalid value encountered in greater
  return np.absolute(m - badval) > atol + rtol * np.absolute(badval)
WARNING: AstropyDeprecationWarning: The new_table function is deprecated and may be removed in a future version.
        Use :meth:`BinTableHDU.from_columns` for new BINARY tables or :meth:`TableHDU.from_columns` for new ASCII tables instead. [astropy.utils.decorators]
WARNING: AstropyDeprecationWarning: The use of header.update() to add new keywords to a header is deprecated.  Instead, use either header.set() or simply `header[keyword] = value` or `header[keyword] = (value, comment)`.  header.set() is only necessary to use if you also want to use the before/after keyword arguments. [astropy.io.fits.header]
NSIDE = 64
ORDERING = RING in fits file
INDXSCHM = EXPLICIT
NSIDE = 64
ORDERING = RING in fits file
INDXSCHM = EXPLICIT
NSIDE = 64
ORDERING = RING in fits file
INDXSCHM = EXPLICIT
NSIDE = 64
ORDERING = RING in fits file
INDXSCHM = EXPLICIT
NSIDE = 64
ORDERING = RING in fits file
INDXSCHM = EXPLICIT
NSIDE = 64
ORDERING = RING in fits file
INDXSCHM = EXPLICIT
NSIDE = 64
ORDERING = RING in fits file
INDXSCHM = EXPLICIT
NSIDE = 64
ORDERING = RING in fits file
INDXSCHM = EXPLICIT

In [8]:
## You should get your 800 noise sims here (100 MC for the 8 frequency channels)
!ls test/test/noise/*.fits | wc 
## And the covariances for all frequency channels + the combination (useful for component separation for example) here
!ls test/test/masks


     800     800   54300
IQU_nside64_test_weights_freq150_HF.fits
IQU_nside64_test_weights_freq150_MF.fits
IQU_nside64_test_weights_freq220_HF.fits
IQU_nside64_test_weights_freq220_UHF.fits
IQU_nside64_test_weights_freq270_UHF.fits
IQU_nside64_test_weights_freq27_LF.fits
IQU_nside64_test_weights_freq39_LF.fits
IQU_nside64_test_weights_freq90_MF.fits
IQU_nside64_test_weights_freq_combined.fits
setup_scanning.ini
test.hdf5

In [9]:
## Let's have a look at one of the Stokes Q noise realisation 
# and the polarised covariance matrix for one frequency channel (that is the inverse noise map in polarisation)
data_27GHz_cov = hp.read_map(
    os.path.join(environment.outpath_masks, 'IQU_nside64_test_weights_freq27_LF.fits'), 1)
data_27GHz_noise = hp.read_map(
    os.path.join(environment.outpath_noise, 'IQU_nside64_test_freq27_LF_white_noise_sim000.fits'), 1)

hp.gnomview(
    data_27GHz_cov, rot=[0,-57.5], xsize=900,reso=6.8, cmap=pl.cm.viridis, 
    title='Covariance ($A^T N^{-1} A$) [27 GHz]', notext=True, unit='$\mu K_{\\rm CMB}^{-2}$', sub=121)

hp.gnomview(
    data_27GHz_noise, rot=[0,-57.5], xsize=900,reso=6.8, cmap=pl.cm.viridis,
    min=-15, max=15,
    title='Stokes Q Noise map [27 GHz]', notext=True, unit='$\mu K_{\\rm CMB}$', sub=122)
hp.graticule(dpar = 1.)
pl.show()


NSIDE = 64
ORDERING = RING in fits file
INDXSCHM = EXPLICIT
NSIDE = 64
ORDERING = RING in fits file
INDXSCHM = EXPLICIT
38.4945419593 141.505458041 -51.5054580407 51.5054580407
The interval between parallels is 10 deg 0.00'.
The interval between meridians is 6 deg 0.00'.
38.4945419593 141.505458041 -51.5054580407 51.5054580407
The interval between parallels is 10 deg 0.00'.
The interval between meridians is 6 deg 0.00'.

In [10]:
## Et voilà! Now you can do the same with different map resolutions, frequency channels, noise levels, etc...