Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

include chromatic noise processes #20

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
58 changes: 58 additions & 0 deletions hasasia/sensitivity.py
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@
'get_Tspan',
'get_TspanIJ',
'corr_from_psd',
'corr_from_psd_chrom',
'quantize_fast',
'red_noise_powerlaw',
'Agwb_from_Seff_plaw',
Expand Down Expand Up @@ -1074,6 +1075,63 @@ def corr_from_psd(freqs, psd, toas, fast=True):
integrand = psd*np.cos(2*np.pi*freqs*tm[:,:,np.newaxis])#df*
return np.trapz(integrand, axis=2, x=freqs)#np.sum(integrand,axis=2)#

def corr_from_psd_chrom(freqs, psd, toas, obs_freqs, chr_idx, fref=1400., fast=True):
"""
Calculates the correlation matrix over a set of TOAs for a given power
spectral density.

Parameters
----------

freqs : array
Array of freqs over which the psd is given.

psd : array
Power spectral density to use in calculation of correlation matrix.

toas : array
Pulsar times-of-arrival to use in correlation matrix.

obs_freqs : array
observation frequency of each ToA

chr_idx : float
Spectral index of the powerlaw amplitude. 2 for DM noise, 4 for Chrom.

fref: float, optional
reference frequency for amplitude powerlaw. Usually set to 1400 MHz

fast : bool, optional
Fast mode uses a matix inner product, while the slower mode uses the
numpy.trapz function which is slower, but more accurate.

Returns
-------

corr : array
A 2-dimensional array which represents the correlation matrix for the
given set of TOAs.
"""

N_toa = len(toas)
matrix = np.ones((N_toa,N_toa))

for i in range(N_toa):
matrix[i,:] = (fref/obs_freqs[i])**chr_idx
A_matrix = matrix*matrix.transpose()

if fast:
df = np.diff(freqs)
df = np.append(df,df[-1])
tm = np.sqrt(psd*df)*np.exp(1j*2*np.pi*freqs*toas[:,np.newaxis])
integrand = np.matmul(tm, np.conjugate(tm.T))
return A_matrix*np.real(integrand)
else: #Makes much larger arrays, but uses np.trapz
t1, t2 = np.meshgrid(toas, toas, indexing='ij')
tm = np.abs(t1-t2)
integrand = psd*np.cos(2*np.pi*freqs*tm[:,:,np.newaxis])#df*
return A_matrix*np.trapz(integrand, axis=2, x=freqs)#np.sum(integrand,axis=2)#

def corr_from_psdIJ(freqs, psd, toasI, toasJ, fast=True):
"""
Calculates the correlation matrix over a set of TOAs for a given power
Expand Down