cbasspy package

Submodules

cbasspy.angles module

Common angle conversions

latLon2ThetaPhi(lat, lon, deg=False)

Convert latitude and longitude to colatitude (theta) and longitude (phi) Default angles all in radians

Parameters:

lat : array_like

Latitude

lon : array_like

Longitude

deg : optional, {True,False}

Whether angles in deg or rad
raDec2GlatGlon(RA, DEC, deg=False, EPOCH=36525.0)

Convert RA and DEC to galactic lat and lon. Default angles in radians Default epoch J2000

Parameters:

RA : float

Right ascension

DEC : float

Declination

deg : optional, {True,False}

Whether angles in deg or rad

EPOCH : ephem object, optional {ephem.J2000}

Epoch of conversions
Returns:

glat : float

Galactic latitude

glon : float

Galactic longitude

cbasspy.bin_llcl module

Function written by PolSPICE

bin_llcl(llcl_in, bin, flatten=False, uniform=False)
turns llcl_in (= continuous l*(l+1)*cl/2Pi) into a binned version with
a constant or variable binwidth ‘bin’
Parameters:

llcl_in : array_like

input l*(l+1)*Cl/2Pi, 1D vector, defined for each l from l=0

bin : array_like,

can be either a scalar = dl or a vector defining the bins boundaries : [low0, low1, low2, …,low(n-2), low(n-1)+1]

flatten : optional, {False,True}

If set, the input C(l) is multiplied internally by l*(l+1)/2Pi before being binned. By default, the input C(l) is binned as is.

uniform : optional, {False, True}

if set, each l is given the same weight in the bin. By default, a weighting propto (2*l+1) (inverse cosmic variance) is applied to each l. In any cases, the output x is the same
Returns:

x : array_like

center of bins

y : array_like

binned l*(l+1)*Cl/2Pi

dx : float (?)

width of each bin

dy : float

returns on output the rms of C(l) for a full sky observation

= C(l) * sqrt( 2/ 2l+1 / dl)

cbasspy.clusteringCode module

Python module to cluster pixels in the sky into regions.

compute_naive_betas(map1_name, map2_name, freq1, freq2, offset1=0.0, offset2=0.0, make_P_from_inputs=False)

Compute the simple pixel-by-pixel spectral index. Input maps should be in the same coordinate system, have the same Nside and be in the same units.

Parameters:

map1_name : str

Name of first map file

map2_name : str

Name of second map file

freq1 : float

Frequency of first map [GHz]

freq2 : float

Frequency of second map [GHz]

offset1 : float, 0.

Zero level offset to add to map1

offset2 : float, 0.

Zero level offset to add to map2
Returns:

beta : array_like

Healpix map of spectral indices
compute_pixel_vectors(Nside)

Computes the vectors to the pixel centres for every pixel in a map.

Returns:

vec : array_like

Array of vectors to the pixel centres
find_meanshift_clusters(beta12, beta23, vec, goodIdx, bandwidth=0.13, min_bin_freq=20)

Assign each pixel to a region using the meanshift algorithm.

Parameters:

X : array_like

Array of values to cluster on

goodIdx : array_like

The indices to use

bandwidth : float, 0.13

Bandwidth to pass to clustering algorith

min_bin_freq : int, 20

Minimum number of pixels in each region
Returns:

labels_map : array_like

Healpix map, with each pixel assigned a region.
find_voronoi_clusters(S_map_name, N_map=None, S_map_offset=0, threshold_score=2.5, min_pixels=20)

Assign each pixel to a region using the Voronoi binning algorithm.

Parameters:

S_map_name : str

Name of input signal map

N_map : str, None

Name of input noise map, if None then noise is set to unity.

S_map_offset : float, 0

Zero level offset to add to signal map

threshold_score : float, 2.5

Target S/N for each region

min_pixels : int, 20

Target minimum number of pixels for each region
Returns:

voronoi_map : array_like

Healpix map, with each pixel assigned a region.

distance_map : array_like

Healpix map, with the distance from each pixel to its region centroid.

cbasspy.const module

TCMB = <cbasspy.const.constant object>

CMB mean temperature

2.7255 K

ap_eff = <cbasspy.const.constant object>

Aperture efficiency

0.55

area_eff = <cbasspy.const.constant object>

Effective dish area

ap_eff.value*dish_area.value m**2

c = <cbasspy.const.constant object>

Speed of light in vacuum

2.99792458e8 m s**-1

class constant(value, unit_string, name)

Bases: object

A numerical constant

Methods

displayValue()

Print the value of the constant and it’s units

dish_area = <cbasspy.const.constant object>

Dish area

(6.1/2)**2*np.pi m**2

h = <cbasspy.const.constant object>

Planck constant.

6.62607004e-34 m**2 kg s**-1

kb = <cbasspy.const.constant object>

Boltzmann constant

1.38064852e-23 m**2 kg s**-2 K**-1

sample_freq = <cbasspy.const.constant object>

Telescope sampling frequency

  1. Hz

cbasspy.conv module

Module containing functions to convert between brightness, temperature, SI units and Jy

brightness2Trj(B, nu, SI=True)

Convert from brightness to brightness temperature in K_RJ.

Brightness temperature is also known as Reyleigh Jeans temperature and sometimes Antenna temperature By default B is assumed to be in SI units, set SI=False for Jy/sr nu in Hz.

Parameters:

B : float

Brightness, in SI units by default.

nu : float

Frequency, in Hz

SI : {True, False}

If True, assumes brightness in SI units. If False, assumes brightness in Jy/sr.
Returns:

T_RJ : Float

Rayleight-Jeans brightness temperature in Kelvin.

See also

trj2Brightness
Inverse conversion
fluxDensityJy2Si(F_Jy)

Converts a flux density from Jy to SI units W m^{-2} Hz^{-1}

Parameters:

F_Jy : Float

Flux density in Jy.
Returns:

F_SI : Float

Flux density in W m^{-2} Hz^{-1}

See also

fluxDensitySi2Jy
Inverse conversion
fluxDensitySi2Jy(F_SI)

Converts a flux density from SI units W m^{-2} Hz^{-1} to Jy

Parameters:

F_SI : Float

Flux density in W m^{-2} Hz^{-1}
Returns:

F_Jy : Float

Flux density in Jy.

See also

fluxDensityJy2Si
Inverse conversion
tcmb2Trj(T_CMB, nu)

Convert from CMB thermodynamic temperature, K_CMB to brightness temperature in K_RJ.

Sometimes K_CMB is called thermodynamic temperature but this could be confusing nu in Hz.

Parameters:

T_CMB : Float

Thermodynamic or CMB temperature in Kelvin.

nu : Float

Frequency in Hz.
Returns:

T_RJ : Float

Rayleigh-Jeans brightness temperature in Kelvin.

See also

trj2Tcmb
Inverse conversion
trj2Brightness(T_RJ, nu, SI=True)

Convert from brightness temperature in K_RJ to brightness.

By default returns brightness in SI units, set SI=False for Jy/sr Brightness temperature also known as Reyleigh Jeans temperature and sometimes Antenna temperature By default B is returned in SI units, set SI=False for Jy/sr nu in Hz.

Parameters:

T_RJ : float

Rayleight-Jeans brightness temperature in Kelvin.

nu : float

Frequency, in Hz

SI : {True, False}

If True, returns brightness in SI units. If False, returns brightness in Jy/sr.
Returns:

B : Float

Brightness, in SI units by default.

See also

brightness2Trj
Inverse conversion
trj2Tcmb(T_RJ, nu)

Convert from brightness temperature in K_RJ to K_CMB.

Sometimes K_CMB is called thermodynamic temperature but this could be confusing. nu in Hz.

Parameters:

T_RJ : Float

Rayleigh-Jeans brightness temperature in Kelvin.

nu : Float

Frequency in Hz.
Returns:

T_CMB : Float

CMB or thermodynamic temperature in Kelvin.

See also

tcmb2Trj
Inverse conversion.

cbasspy.createClusters module

cbasspy.diffuseComponentModels module

Python module to simulate the frequency spectra of various emission mechanisms.

add_white_noise(T, rms)

Adds white noise to temperatures

Parameters:

T : array_like

signal to add noise to

rms : array_like

sigma of gausian from which to draw noise
Returns:

Tsim : array_like

T+noise realization
cmb_model(amp, T0, f)

CMB model

T_RJ = amp * x^2*exp(x)/(exp(x)-1^2 where x = h*nu/(k*T0)

Parameters:

amp : array_like

amplitude of CMB fluctuation [K_CMB].

T_0 : array_like

Temperature of CMB [K]

f : array_like

frequency at which to return amplitude [Hz]
Returns:

T_RJ : array_like

[K_RJ] at frequency f.
curved_power_law_model(amp, beta, C, f0, f, factorHalf=False)

Curved Power Law model

if factorHalf==False:
T_RJ = amp * (nu/nu0)^(beta+C*log(f/f0))
elif factorHalf==True:
T_RJ = amp * (nu/nu0)^(beta+0.5*C*log(f/f0))
Parameters:

amp : array_like

amplitude at reference frequency

beta : array_like

spectral index

C : array_like

Curvature

f0 : array_like

reference frequency [Hz]

f : array_like

frequency at which to return amplitude [Hz]

factorHalf : {False, True}

Whether to multiply C by 0.5.
Returns:

T_RJ : array_like

Amplitude at frequency f in same units as amp.
free_free_model(Te, EM, f)

Free-free model

Two parameter model of Draine - put in a reference!

Parameters:

Te : array_like

Effective electron temperature in beam [K]

EM : array_like

Effective emission measure in beam [cm^6 pc^-1]

f : array_like

Frequency at which to return amplitude [Hz]
Returns:

T_RJ : array_like

[K] at frequency f.
interpolater = <scipy.interpolate.interpolate.interp1d object>

scipy interp1d function for the SPDUST2 model

power_law_model(amp, beta, f0, f)

Power Law model

T_RJ = amp * (nu/nu0)^beta

Parameters:

amp : array_like

amplitude at reference frequency

beta : array_like

spectral index

f0 : array_like

reference frequency [Hz]

f : array_like

frequency at which to return amplitude [Hz]
Returns:

T_RJ : array_like

Amplitude at frequency f in same units as amp.
spdust2_model(fPeak, amp, f0, f)

Spinning dust model

SPDUST2 in a cold neutral medium shifted in frequency and amplitude space. -Add a reference-

Parameters:

fPeak : float

Frequency of peak emission [Hz]

amp : float

Amplitude at reference frequency f0 [K]

f : array_like

Freq at which to return amplitude [Hz]
Returns:

T_RJ : array_like

[K] at frequency f.
synch_moment_expansion(amp, beta, moments, f0, f)

Created – Richard David Paradorn Grumitt (09/05/2018)

Adapted by Jaz (22/05/2018)

This model implements the moment expansion SED for a power law superposition described in J. Chluba et. al. 2017. Model employed here assumes power law amplitudes and spectral indices are not correlated.

Adapted from the PySM synchrotron class method, which can be used to generate whole sky simulations with instrumental and noise effects.

Parameters:

amp : float

Amplitude

beta : float

spectral index

moments : array_like

moments of the synchrotron SED up to desired order, as defined in J. Chluba et. al. 2017. – numpy.ndarray

f0 : float

Reference frequency in Hz

f : array_like

frequency at which to evualate spectrum
Returns:

T_RJ : array_like

Power law moment expansion SED for I, Q and U, evaluated at frequency nu
thermal_dust_model(amp, beta, T, f0, f)

Thermal dust model

T_RJ = amp*(f/f0)^(beta+1)*(exp(gamma*f0)-1)/(exp(gamma*f)-1.) where gamma = h/(kb*T)

Parameters:

amp : array_like

amplitude of dust [X] at reference frequency

beta : array_like

emissivity index

T : array_like

Dust temperature [K]

f0 :array_like

Reference frequency [Hz]

f : array_like

frequency at which to return amplitude [Hz]
Returns:

T_RJ : array_like

[X] at frequency f.

cbasspy.fileListBash module

A slightly ugly set of wrappers that call commonly useful bash commands for dealing with file lists.

createDailyFileLists(start_date_str, number_of_days, fits_dir, suffix='')

Create daily file lists from all the fits files in one directory.

Outputs files as DD-MMM-YYYY_files[suffix].txt

Parameters:

start_date_str : str

start date in format DD-MMM-YYYY, e.g. ‘10-Sep-2010’

number_of_days : int

Number of subsequent days to create file lists for, e.g. 20

fits_dir : str

Directory containing TOD .fits files

suffix : str, optional {‘’}

extra suffix to add to output file names
createMonthlyFileLists(start_date_str, number_of_months, fits_dir, suffix='')

Create monthly file lists from all the fits files in one directory.

Outputs files as MMM-YYYY_files[suffix].txt

Parameters:

start_date_str : str

start date in format MMM-YYYY, e.g. Sep-2010

number_of_months : int

Number of subsequent months to create file lists for, e.g. 20

fits_dir : str

Directory containing TOD .fits files

suffix : str, optional {‘’}

extra suffix to add to output file names
excludeDateRange(file_list_name, output_name, exclusion_min_str, exclusion_max_str)

Exclude a date range from a file list.

Parameters:

file_list_name : str

Name of file list

output_name : str

Where to save the outputs. Don’t include the file extension .txt, the code adds it.

exclusion_min_str : str

Earliest date to exclude, format e.g. 01-Aug-2014:00:00:01

exclusion_max_str : str

Latest date to exclude, formate e.g. 01-Oct-2014:00:00:01
excludeLinesWithMatches(fileListName, matchesListName, outListName, verbose=False)

Function to remove files from fileListName that contain any string in the file matchesListName.

Parameters:

file_list_name : str

Name of file list

matchesListName : str

Name of list of matches

outListName: str

What to save the output to

verbose : bool, False

Whether to print messages to screen.
makeRaDecList(fileList, RAmin, RAmax, DECmin, DECmax, SAMPLESmin, verbose=0)

Give functiona a list of files and an RA/DEC region (in radians). Fits files with more than SAMPLESmin number of samples within the region are outputted in a goodFileList.

Parameters:

fileList : list

list of file names

RAmin : float

Minimum RA to accept [radians]

RAmax : float

Maximum RA to accept [radians]

DECmin : float

Minimum DEC to accept [radians]

DECmax : float

Maximum DEC to accept [radians]

SAMPLESmin : int

Minimum number of samples within the region for the file to be passed to the good list.

verbose : optional {0}

If >=1 then print summary to screen.
Returns:

goodFileList : list

List of files that have met the criteria
openFileList(fileListName)

Open the .txt file fileListName. Return a list of strings, each element a line from the file

This is only useful if you can’t remember how to open a txt file with python!

Parameters:

fileListName : str

Name of file to open.
Returns:

fileList : list

List of file names

See also

saveFileList

prependString(string2prepend, inFileName, outFileName, verbose=1)

prependString will prepend the string string2prepend to the start of each line of a file in inFileName. The output is saved to outFileName. If verbose>0 then a summary will be printed to the screen.

Parameters:

string2prepend : str

String to prepend to each line of input file

inFileName : str

Name (including path) of file to append string to

outFileName : str

Name (including path) of file to save output to

verbose : optional {1,0}

If >1 then a summary will be printed to screen.
saveCommonLines(inFileName1, inFileName2, outFileName, verbose=1.0)

saveCommonLines looks for matching lines in inFileName1 and inFileName2. Matching lines will be written to outFileName. If verbose>0 then a summary will be printed to the screen.

Parameters:

inFileName1 : str

Name (including path) of first file

inFileName2 : str

Name (including path) of second file

outFileName : str

Name (including path) of file to save output to

verbose : optional {1,0}

If >1 then a summary will be printed to screen.
saveFileList(outListName, fileList)

Save a list of files to a txt file called outListName

This is only useful it you can’t remember how to save a txt file with python!

Parameters:

outListName : str

Where to save the file list

fileList : list

list of lines to save

See also

openFileList

splitFileListCentralDate(file_list_name, output_name)

Split a file list into two (potentially) unequal halves. The central date, halfway between the first and last dates in the file, is the boundary between files that are outputted into one file and another.

Parameters:

file_list_name : str

Name of file list

output_name : str

Where to save the outputs. Don’t include the file extension .txt, the code adds it.
splitFileListSortDate(file_list_name, output_name, number_of_splits=2)

Split a file list into equal sections. Each containing 1/number_of_splits of the chronologically oredered files.

Parameters:

file_list_name : str

Name of file list

output_name : str

Where to save the outputs. Don’t include the file extension .txt, the code adds it.

number_of_splits : int, 2

Number of splits to make to the data
splitOddAndEvenLines(inFileName, evenFileName, oddFileName, verbose=1)

Split a file list into two files, one listing even lines and the other listing odd lines.

Parameters:

inFileName : str

File to split

evenFileName : str

where to save the even lines

oddFileName : str

where to save the odd lines

verbose : optional {1}

If >1 print a summary to screen

cbasspy.fit1line module

TruncatedOneOverX_like(value, a, b)

Depricated. In future this will be replaced by cbasspy.pymcCustomPriors.BoundedOneOverX_like().

fit1line(x_obs, y_obs, x_sigma, y_sigma, init_slope_1=0.1, niter=500000.0, nburn=100000.0, nthin=20, ntune=500.0, seed=1234, rangeLower_x=1e-06, rangeUpper_x=2000.0, rangeLower_y=1e-06, rangeUpper_y=2000.0, locLower_x=-2000.0, locUpper_x=1000.0, locLower_y=-2000.0, locUpper_y=1000.0, save2db=False, dbName='mcmcChains_fit1line.db3')

Fit a straight line with outliers to data with errors in x and y. The problem is parametrized according to Gull 1989, this method produces unbiased estimates of the parameters of a straight line. It assumes errors on x_obs and y_obs are Gaussian. The line is reparametrized as hat{Y}=pm hat{X}, hat{X}=(x-x0)/Rx and hat{Y}=(y-y0)/Ry the loc parameters x0 and y0 are given uniform priors, the range parameters Rx and Ry are given 1/R priors. The true locations of the x values are marginalised over. (See Gull 1989 for details)

The outlies are modelled using a mixtude model; y = a*x+b OR 2D Gaussian with no covariance. The 2D Gaussian is paramterized by mean (x0_3, y0_3) and covariance ((Rx_3^2,0),(0,Ry_3^2))

Parameters:

x_obs : array_like

x values

y_obs : array_like

y values

x_sigma : array_like

errorbar on x

y_sigma : array_like

errorbar on y

init_slope_1 : float {0.1}

Depricated, kept for back compatablility. Parameter was originally Starting guess for slope of line, now does nothing.

niter : int, optional {5e5}

Number of steps in MCMC chain

nburn : int, optional {1e5}

Number of steps to burn

nthin : int, optional {20}

Factor by which to thin the chain

ntune : int, optional {500}

Number of steps between tuning of step size (only during burn in)

seed : int, optional {1234}

Random seed

rangeLower_x : float, optional {1e-6}

Lower limit of x range parameter [Must be >0]

rangeUpper_x : float, optional {2e3}

Upper limit of x range parameter [Must be >0]

rangeLower_y : float, optional {1e-6}

Lower limit of y range parameter [Must be >0]

rangeUpper_y : float, optional {2e3}

Upper limit of y range parameter [Must be >0]

locLower_x : float, optional {-2e3}

Lower limit of x location parameter [float]

locUpper_x : float, optional {1e3}

Upper limit of x location parameter [float]

locLower_y : float, optional {-2e3}

Lower limit of y location parameter [float]

locUpper_y : float, optional {1e3}

Upper limit of y location parameter [float]

save2db : boolean {False}

Whether to save output as a sqlite database

dbName : str {‘mcmcChains_fit1line.db3’}

Name to save the database as
Returns:

mcmc : pymc object

The function returns a PyMC object that contains the traces of

  • The free parameters in the model,
  • The straight line component x0, y0, Rx, Ry
  • model_a_1 and model_b_1 are the calculated slope and y-intercept for each step in the chain
  • The outlier component x0_3, y0_3, Rx_3, Ry_3, frac_3
  • (frac_3 is the fraction of points that belong to the outlier component)
  • True x values x

To access these chains,:

mcmc.trace('x0_3')[:]

Notes

  • The lower and upper limits on the range and loc parameters should be set far away from the expected parameter range and central values and decided

on before looking at the data. * If the input parameters result in zero probability, None is returned.

cbasspy.fit_mean_var module

intrinsic_scatter(obs, obs_err, mean_lower, mean_upper, sd_lower, sd_upper, nsamples=2000, ntune=1000, accept_prob=0.9, jeffreys_prior=False)

Function for fitting the mean and additional intrinsic scattter for a dataset.

Returns:

trace: pymc3.MultiTrace

Output trace from your fitting - pymc3.MultiTrace

cbasspy.fitexey module

This module may be incorrect, it definitely needs checking.

fitexey(x_obs, y_obs, x_sigma, y_sigma, niter=500000.0, nburn=10000.0, nthin=20, ntune=100.0, dbname='./model_db.pickle', seed=1234)

Fit a straight line to data with errors in x and y. y = ax+b Returns traces of the slope, a, and intercept, b. This method is unbiased - Gull 1989. x_obs: x values y_obs: y values x_sigma: errorbar on x y_sigma: errorbar on y niter: Number of steps in MCMC chain nburn: Number of steps to burn nthin: Factor by which to thin the chain ntune: Number of steps between tuning of step size (only during burn in) dbname: Not currently used seed: Random seed

cbasspy.masks module

produce_dec_mask(NSIDE, rotTo, cutMin, cutMax)

Return a declination mask, masking region between cutMin and cutMax Pixels are set to 0 within the mask and 1 outside.

Parameters:

NSIDE : int

Nside of mask, e.g. NSIDE=512

rotTo : str, {‘C’,’G’,’E’}

Coordinate system to rotate mask to.

cutMin : float

Lower declination bound, e.g. cutMin=-90.

cutMax : float

Upper declination bound, e.g. cutMax=0.
Returns:

fullR : array_like

HEALPix map.
produce_frac_mask(inMap, frac, account_for_masked=False)

Return a mask excluding the brightest fraction of the inMap. The masked pixels are set to 0, unmasked pixels to 1.

Parameters:

inMap : array_like

HEALPix sky map

frac : float

Fraction of the sky to mask.

account_for_masked : boolean, False

frac of unmasked pixels if True, else frac of all pixels
Returns:

mask : array_like

HEALPix map.
produce_plane_mask(NSIDE, rotTo, cutOff)

Return a mask of the galactic plane. Pixels are set to 0 if too close to the plane and 1 if far enough away.

Parameters:

NSIDE : int

Nside of mask, e.g. NSIDE=512

rotTo : str, {‘C’,’G’,’E’}

Coordinate system to rotate mask to.

cutOff : float

cutOff in degrees from galactic plane.
Returns:

fullR : array_like

HEALPix map.

cbasspy.noise module

fitNoiseModelToTodFast(freq, I_N, niter=5000.0, verbose=True)

Fits a white + red noise model to a periodogram using a normal approximation to the posterior. Draws posterior from samples drawn from the normal approx.

Parameters:

freq : array_like

Frequencies [Hz]

I_N : array_like

Periodogram

niter : int

Number of steps in chain
Returns:

Mean[alpha] : float

Mean[sigma_r] : float

Mean[sigma_w] : float

Cov[alpha,alpha] : float

Cov[alpha,sigma_r] : float

Cov[alpha,sigma_w] : float

Cov[sigma_r,sigma_r] : float

Cov[sigma_r,sigma_w] : float

Cov[sigma_w,sigma_w] : float

See also

fitNoiseModelToTodFitsFile
Does full MCMC
fitNoiseModelToTodFitsFile(freq, I_N, niter=50000.0, nburn=25000.0, thinFactor=1, ntune=1000.0, verbose=False)

Fits a white + red noise model to a periodogram using the full Whittle Pseudo-Likelihood. Method described in Luke Jew’s Thesis, Chapter 3.

Parameters:

Data parameters

^^^^^^^^^^^^^^^

freq : array_like

Frequency in Hz

I_N : array_like

Periodogram

MCMC parameters

^^^^^^^^^^^^^^^

niter : int

Number of steps in chain

nburn : int

Number of initial steps to burn

thinFactor : int

Thinning factor

ntune : int

How often to tune step size
Returns:

Mean[alpha] : float

Max posterior alpha : float

Mean[sigma_r] : float

Max posterior sigma_r : float

Mean[sigma_w] : float

Max posterior sigma_w : float

Mean[f_knee] : float

Max posterior f_knee : float

Var[f_knee] : float

Cov[alpha,alpha] : float

Cov[alpha,sigma_r] : float

Cov[alpha,sigma_w] : float

Cov[sigma_r,sigma_r] : float

Cov[sigma_r,sigma_w] : float

Cov[sigma_w,sigma_w] : float

chi2 : float

dof : int

See also

ln_noise_model
Noise model that is fit to the data.
cbasspy.const
Fundamental constants
cbasspy.angles
Angle conversions to subtract sky signal.
leastsq_fitNoiseModelToTodFitsFile(hdulist, chanName, model_sky_map, sample_freq=100.0, minChunkLength=100, savetxt=False, savetxtStub='example')

Least squares fitting method. Needs to average lots of chunks together to estimate P_xx with normal error bars

p[0] = alpha p[1] = f_knee p[2] = sigma_w
Parameters:

hdulist : pyfits object

TOD Fits file

chanName : {‘I1’,’I2’,’Q1’,’Q2’,’U1’,’U2’}

Which channel to fit to

model_sky_map : array_like

HEALPix model sky map

sample_freq : float, optional {const.sample_freq.value}

Sample frequency of TOD

minChunkLength : int, optional {100}

Minimum chunk length to fit to

savetxt = {False,True}

Whether to save outputs

savetxtStub : str, optional {‘example’}

Name to attach to output files.
Returns:

p : array_like

Best fitting parameters, alpha, f_knee, sigma_w

mesg : str

Fitting string

flag : int
ln_noise_model(p, f, f_0=0.1)

Returns the natural logarithm of a white + 1/f noise model power spectrum.

ln_P_xx = ln_noise_model(p,f,f_0=0.1)
= np.log(p[2]**2 + p[1]**2*np.power(f/f_0,-p[0])))
Parameters:

p : array_like

Noise parameters shape = [3,]. p[0] = alpha, p[1] = sigma_r, p[2] = sigma_w

f : array_like

Frequency in Hz
Returns:

ln_P_xx : array_like

Noise power spectrum.
minimizeParams(f, ln_P_xx, p_guess)

Parameters needed to minimise the rms value of the difference between ln_P_xx and the power spectrum computed by ln_noise_model. residsFunc is used to calculate the difference.

Parameters:

f : array_like

Frequency in Hz

ln_P_xx : array_like

Natural logarithm of power spectrum.

p_guess : array_like

The starting values of the parameters, shape = [3,]. p_guess[0] = alpha p_guess[1] = sigma_r p_guess[2] = sigma_w

See also

residsFunc
The residuals function.
old_ln_noise_model(p, f)

Returns the natural logarithm of a white + 1/f noise model power spectrum.

Please note, this function parametrises the noise model with f_knee - not sigma_r ln_P_xx = ln_noise_model(p,f)

= np.log(p[2]**2 + p[1]**2*np.power(f,-p[0])))
Parameters:

p : array_like

Noise parameters shape = [3,]. p[0] = alpha, p[1] = f_knee, p[2] = sigma_w

f : array_like

Frequeny in Hz.
residsFunc(p, f, ln_P_xx)

Returns the difference between ln_P_xx and the power spectrum computed by ln_noise_model(p,f) with params p.

p[0] = alpha p[1] = sigma_r p[2] = sigma_w
Parameters:

p : array_like

Noise parameters

f : array_like

Frequency in Hz

ln_P_xx : array_like

Natural logarithm of power spectrum to compare to.
Returns:

diff : array_like

Difference

See also

ln_noise_model
Function called

cbasspy.noise_fit_12 module

calcSpectralContaminantChisq(f, I_N, ps, freqLower, freqUpper)
noise_fit_12(outputFileName, fileListName, model_sky_name, convert2Kelvin=1.0, doPol=1.0, niter=20000, nburn=10000, thinFactor=1, ntune=500, cutName='NIGHT', spikeThreshold=0.035, mpi_on=True)

Calculates noise stats for I1/2,Q1/2 and U1/2 Previous versions averaged these to I, Q and U. noise_fit_12(outputFileName, fileListName, model_sky_name, convert2Kelvin=1., doPol =1.,

niter=int(20e3), nburn=int(10e3), thinFactor=int(1), ntune=int(5e2), cutName=’NIGHT’, spikeThreshold=0.035 mpi_on=True):
Parameters:

outputFileName : str

Name of db3 file to create and save results in. If outputFileName already exists then a table called noiseStats is created. If noiseStats already exists then it is overwritten

fileListName : str

Name of txt file that lists the TOD fits files

model_sky_name : str

Name of healpix map .fits file that is used to generate an expected sky signal. MAP MUST BE IN GALACTIC COORDINATES

convert2Kelvin : float, (1.)

Scaling factor by which the TOD and map are multiplied by

doPol : int, {1,0}

if 0 only I, if 1 IQU

niter : int

number of steps to take in MCMC chains

nburn : int

number of steps to burn

thinFactor : int

Factor by which to thin the chains

ntune : int

how frequently to tune the step proposal distribution during tuning

mpi_on : boolean, {(True),False}

Whether this should be run on multiple cores.

cutName : str

Name of extension in fits file, default ‘NIGHT’

spikeThreshold : float, (0.035)

If model sky changes by >spikeThreshold between 2 samples then don’t use that chunk. Current default is 0.035, it used to be 0.025
Returns:

nothing :

cbasspy.overwriteFitsNoiseStats module

ASrSw2Fk(alpha, sigma_r, sigma_w, f_0=0.1)

Take alpha, sigma_r, sigma_w and f_0 (default 0.1Hz) Return f_knee

findFitsExtNo(hdulist, cutName)

extNo = findFitsExtNo(hdulist,cutName) Find the header extension with name cutName. Returns extNo, hdulist[extNo] is the table of interest.

overwriteFitsNoiseStats_fromSql(fileListName, inSqlFileName, cutName='NIGHT')

Function to overwrite the noise statistics in the fits files listed in fileList with the noise properties in the SQL database.

Parameters:

fileListName : str

Name of the .txt file containing the file list

inSqlFileName : str

Name of the SQL databse

cutName : str, ‘NIGHT’

Which cut extension to overwrite
Returns:

filesWithoutNoiseStats : array_like

A list of files where the default noise properties were written.

filesWithNoiseStats : array_like

A list of files where a match was found in the database.
wavg_12_calc(inSqlFileName, outSqlFileName)

Script to calculate the weighted average noise statistics for each fits file. Takes a sqlite database as input. Outputs results in a second sqlite database.

weightedNoiseStatsMean(aa, ar, aw, rr, rw, ww, alpha_mu, sigma_r_mu, sigma_w_mu, error=False)

weightedMean [,covariance] = weightedNoiseStatsMean(aa,ar,aw,rr,rw,ww,alpha_mu,sigma_r_mu,sigma_w_mu,error=False) Calculate the weighted average noise properties, accounting for parameter covariance

aa,ar,aw,rr,rw,ww: arrays of (co)variance of alpha, sigma_r and sigma_w alpha_mu: array of alpha values sigma_r_mu: array of sigma_r values sigma_w_mu: array of sigma_w values error: If True just return weighted average, if False also return covariance.

cbasspy.pymcCustomPriors module

BoundedNegOneOverX_like(value, lower, upper)
BoundedOneOverX_like(value, lower, upper)
p(X)~1/X when lower<X<upper
0 else

Returns log probability

JeffCuPrior_like(value, sigs, nus, nu0, lower, upper, beta)
JeffEmPrior_like(value, sigs, nus, lower, upper, Te)
JeffNupPrior_like(value, sigs, nus, f0, Y, lower, upper)
JeffSiPrior_like(value, sigs, nus, nu0, lower, upper, G)

Jeffreys prior for spectral indices when likelihood is L~1/sigma * exp(-(X-A*G*(nu/nu0)^beta)/2sigma^2) Where A is an amplitude and G is a fixed function of frequency e.g. for power law with no curvature. G = 1.0 e.g. for a power law with curvature, G = exp(0.5*C*ln(nu_i/nu_0)) e.g. for modified black body (i.e. thermal dust) G = (nu_i/nu_0) * [exp(h*nu_0/k*T)-1]/[exp(h*nu_i/k*T)-1]

Args: sigs = array of uncertanties nus = array of frequencies nu0 = reference frequency lower = lower bound on spectral index upper = upper bound on spectral index value = proposed value of spectral index G = either float or an array of floats.

Returns: ln(probability)

JeffSpectralPrior_like(value, sigs, nus, nu0, lower, upper, Y, X, factor)
JeffTdPrior_like(value, sigs, nus, nu0, lower, upper, Y)
TruncatedExpX_like(value, lower, upper)
TruncatedX_like(value, lower, upper)

cbasspy.scalePolCovsBy2 module

scalePolCovsBy2(cov_name, output_str='_polCovCor')

Scale the Q and U covariances by 2. Descart thinks that Q1 and Q2 (U1 and U2) are independent - they are not. If you make a map using both 1 and 2 channels then you will need to apply this scaling.

Parameters:

cov_name : str

Name of covariance fits file

output_str : str, ‘_polCovCor’

String to attached to output file name. The output file name takes the form:

'{0}{1}.fits'.format(cov_name[0:-5],output_str)
Returns:

cov_map : array_like

Corrected healpix map.

cbasspy.simulate module

groundFromTemplate(az, templateFileName)

Simulates the ground signal at a given azimuth and elevation by taking a sample ground template.

Parameters:

az : array_like

azimuths to sample at [degrees]

templateFileName : str

Template filename (and directory).
Returns:

signal : array_like

Simulated ground TOD

Notes

LJ - 05/07/18 AJB - 05/07/18

overwrite_tod_fits_file(fitsFile, whichField, values)

Overwrites the field [whichField] in extension 1 of the TOD Fits file [fitsFile] with the values in [values].

Parameters:

fitsFile : str

File name

whichField : str

Field name

values : float, array_like

values to overwrite field with
polSkyFromMap(RA, DEC, MJD, LAT, LON, Q_sim_map, U_sim_map)

Simulate stokes Q and U signals from input Q and U maps ~~in Celestial coordinates~~.

Future work Support input maps in other coordinate systems (e.g. Galactic) Support angles in degrees

Parameters:

RA : array_like

Right ascension in radians

DEC : array_like

Declination in radians

MJD : array_like

MJD

LAT : float

Latitude of telescope, in radians

LON : float

Longitude of telescope, in radians

Q_sim_map : array_like

Simulated Q Healpix sky map, in Celestial coordinates

U_sim_map : array_like

Simulated U Healpix sky map, in Celestial coordinates
Returns:

Q_a : array_like

Rotated Q signal in telescope coordinates

U_a : array_like

Rotated U signal in telescope coordinates
skyFromMap(glat, glon, sky_sim_map)

Simulate sky-signal TOD from the HEALPix map sky_sim_map (in Galactic coordiantes) at coordinates glat and glon (in radians).

Parameters:

glat : array_like

array of Galactic latitudes [radians]

glon : array_like

array of Galactic longitudes [radians]

sky_sim_map : array_like

HEALPix sky map.
Returns:

signal : array_like

Simulated sky signal
whiteAndRedNoiseSim(sigma_w, f_knee, alpha, sample_freq, N, return_pdf=False)

nosie (,f,P) = whiteAndRedNoiseSim(sigma_w,f_knee,alpha,sample_freq,N,return_pdf=False) Simulate red and white noise in TOD.

Parameters:

sigma_w : float

white noise level in K

f_knee : float

knee frequency in Hz

alpha : float

power law

sample_freq : float

for C-BASS this is 100 Hz

N : int

length of noise to simulate
Returns:

noise : array_like

Simulated noise signal

cbasspy.smart_ud_grade module

alm_ud_grade(map_in, nside_out)

Upgrade or downgrade map in alm space.

Parameters:

map_in : array_like

input map

nside_out : int

output resolution
Returns:

map_out : array_like

up/down graded map
smart_ud_grade(map_in, cov_in, nside_out, pess=False)

Upgrade or downgrade the resolution of a healpix map

When degrading the resolution, the value of the output superpixel is set to the inverse variance weighted mean of the input children pixels. When upgrading the resolution, the value of the children pixels are set to that of the superpixel.

When degrading the resolution, the input inverse variances are summed to give the superpixel. When upgrading the resolution, the input inverse variance is divided by the number of pixels.

The default healpy.ud_grade calculates the average of children pixels. Multiplying the average by the number of daughter pixels in a superpixel returns their sum. ud_grade has the option to multiply the outputs by the number of children, setting power=-2.

Parameters:

map_in : array_like

input map

cov_in : array_like

input covariance

nside_out : int

desired resolution of output map
Returns:

map_out : array_like

up/down graded map

cov_out : array_like

up/down graded cov