Example of running pysm on partial sky

There are two main modules to import when using PySM, pysm and models:



In [ ]:

    
import pysm

import numpy as np
import healpy as hp
import matplotlib.pyplot as plt
%matplotlib inline



In [ ]:

    
from pysm.nominal import models

Sky

Sky configuration dictionary. The keys specify which components will be present; the items are lists containing configuration dictionaries for each population of a component. Here we use some of the predefined models provided with PySM which can be accessed from the models module.

pixel_indices is a numpy array of integer pixel indices in RING ordering.



In [ ]:

    
nside = 64
npix = hp.nside2npix(nside)
pixel_indices = np.arange(npix//2-4000, npix//2+4000, dtype=np.int)
sky_config = {
    'synchrotron' : models("s1", nside, pixel_indices=pixel_indices),
    'dust' : models("d8", nside, pixel_indices=pixel_indices),
    'freefree' : models("f1", nside, pixel_indices=pixel_indices),
    'cmb' : models("c1", nside, pixel_indices=pixel_indices),
    'ame' : models("a1", nside, pixel_indices=pixel_indices),
}



In [ ]:

    
# initialise Sky 
sky = pysm.Sky(sky_config)



In [ ]:

    
len(pixel_indices)



In [ ]:

    
len(sky.ame(nu=40)[0])



In [ ]:

    
def build_partial_map(pixel_indices, pixel_values):
    partial_map = np.nan * np.empty(hp.nside2npix(nside))
    partial_map[pixel_indices] = pixel_values
    return partial_map



In [ ]:

    
ame_map = build_partial_map(pixel_indices, sky.ame(nu=40)[0])
hp.mollview(ame_map, title="AME partial map", min=0, max=3000)



In [ ]:

    
# Sky then has attributes corresponding to each of the sky
# components which return that component for a given frequency:
cmb = sky.cmb(nu = 23.)

fig = plt.figure(figsize = (13, 8))
hp.mollview(build_partial_map(pixel_indices, cmb[0]), min = -250, max = 250, title = 'cmb I', sub = (131))
hp.mollview(build_partial_map(pixel_indices, cmb[1]), min = -5, max = 5, title = 'cmb Q', sub = (132))
hp.mollview(build_partial_map(pixel_indices, cmb[2]), min = -5, max = 5, title = 'cmb U', sub = (133))
plt.show()



In [ ]:

    
dust = sky.dust(300.)
fig = plt.figure(figsize = (13, 8))
hp.mollview(build_partial_map(pixel_indices, dust[0]), min = 0, max = 400, title = 'Dust I', sub = (131))
hp.mollview(build_partial_map(pixel_indices, dust[1]), min = -5, max = 10, title = 'Dust Q', sub = (132))
hp.mollview(build_partial_map(pixel_indices, dust[2]), min = -5, max = 10, title = 'Dust U', sub = (133))
plt.show()



In [ ]:

    
# or the total signal as a function of frequency:
total = sky.signal()(100.)
fig = plt.figure(figsize = (13, 8))
hp.mollview(build_partial_map(pixel_indices, total[0]), min = -250, max = 250, title = 'Total I', sub = (131))
hp.mollview(build_partial_map(pixel_indices, total[1]), min = -5, max = 5, title = 'Total Q', sub = (132))
hp.mollview(build_partial_map(pixel_indices, total[2]), min = -5, max = 5, title = 'Total U', sub = (133))
plt.show()

We can also evaluate the signal for a vector of frequencies:

Instrument

If we want to add instrument effects then use the Instrument class. Since this isn't easily a function of frequency it writes maps to file for a given set of observation frequencies, noise, bandpass etc.

Instantiating the Instrument object requires a configuration dictionary, as in the case of Sky.

Delta Bandpasses

First we will consider the case of delta bandpasses, and so we use the frequencies key to specify the observation frequencies.

We need to pass the pixel_indices array also to the instrument configuration dictionary.



In [ ]:

    
instrument_delta_bpass = {
    'frequencies' : np.logspace(1., 3., 20),
    'beams' : np.ones(20)*70.,
    'sens_I' : np.ones(20),
    'sens_P' : np.ones(20),
    'nside' : nside,
    'noise_seed' : 1234,
    'use_bandpass' : False,
    'add_noise' : True,
    'output_units' : 'uK_RJ',
    'use_smoothing' : True,
    'output_directory' : './',
    'output_prefix' : 'test',
    'pixel_indices' : pixel_indices
}

instrument = pysm.Instrument(instrument_delta_bpass)
instrument.observe(sky)



In [ ]:

    
%ls *.fits



In [ ]:

    
hp.mollview(hp.read_map("test_nu0010p00GHz_total_nside0064.fits"), min=0, max=4e4, title="10 GHz")



In [ ]:

    
hp.mollview(hp.read_map("test_nu1000p00GHz_total_nside0064.fits"), min=0, max=4e4, title="1000 GHz")

Arbitrary Bandpasses

Next we consider tophat bandpasses, which we give as a list of tuples: [(frequencies_1, weights_1), (frequencies_2, weights_2) ... ] to the channel key. We must also provide names of these channels as a list of strings to the channel_names key. In this case we do not need to set frequencies since all outputs are defined by the bandpasses.



In [ ]:

    
instrument_bpass = {
    'use_smoothing' : True,
    'beams' : np.array([30., 40., 50.]), # beam fwhm in arcmin
    'nside' : nside,
    'add_noise' : True,
    'sens_P' : np.array([4.5, 4.2, 3.1]), # channel sensitivities in uK_CMB amin
    'sens_I' : np.array([3.5, 3.2, 2.1]),
    'use_bandpass' : True,
    'channels' : [(np.linspace(20, 25, 10), np.ones(10)), (np.linspace(30, 35, 10), np.ones(10)), (np.linspace(40, 45, 10), np.ones(10))],
    'channel_names' : ['channel_1', 'channel_2', 'channel_3'],
    'output_units' : 'uK_RJ',
    'output_directory' : './',
    'output_prefix' : 'test',
    'noise_seed' : 1234,
    'pixel_indices' : pixel_indices
}

instrument = pysm.Instrument(instrument_bpass)
instrument.observe(sky)