# -*- coding: utf-8 -*-
import numpy as np
from enterprise import constants as const
from enterprise.signals import (deterministic_signals, parameter, signal_base,
utils)
[docs]def fdm_block(Tmin, Tmax, amp_prior='log-uniform', name='fdm',
amp_lower=-18, amp_upper=-11,
freq_lower=-9, freq_upper=-7,
use_fixed_freq=False, fixed_freq=-8):
"""
Returns deterministic fuzzy dark matter model:
1. FDM parameterized by frequency, phase,
and amplitude (mass and DM energy density).
:param Tmin:
Min time to search, probably first TOA (MJD).
:param Tmax:
Max time to search, probably last TOA (MJD).
:param amp_prior:
Prior on log10_A.
:param logmin:
log of minimum FDM amplitude for prior (log10)
:param logmax:
log of maximum FDM amplitude for prior (log10)
:param name:
Name of FDM signal.
:param amp_upper, amp_lower, freq_upper, freq_lower:
The log-space bounds on the amplitude and frequency priors.
:param use_fixed_freq:
Whether to do a fixed-frequency run and not search over the frequency.
:param fixed_freq:
The frequency value to do a fixed-frequency run with.
"""
# BWM parameters
amp_name = '{}_log10_A'.format(name)
log10_A_fdm = parameter.Uniform(amp_lower, amp_upper)(amp_name)
if use_fixed_freq is True:
log10_f_fdm = fixed_freq
if use_fixed_freq is False:
freq_name = '{}_log10_f'.format(name)
log10_f_fdm = parameter.Uniform(freq_lower, freq_upper)(freq_name)
phase_e_name = '{}_phase_e'.format(name)
phase_e_fdm = parameter.Uniform(0, 2*np.pi)(phase_e_name)
phase_p = parameter.Uniform(0, 2*np.pi)
fdm_wf = fdm_delay(log10_A=log10_A_fdm, log10_f=log10_f_fdm,
phase_e=phase_e_fdm, phase_p=phase_p)
fdm = deterministic_signals.Deterministic(fdm_wf, name=name)
return fdm
[docs]def cw_block_circ(amp_prior='log-uniform', dist_prior=None,
skyloc=None, log10_fgw=None,
psrTerm=False, tref=0, name='cw'):
"""
Returns deterministic, cirular orbit continuous GW model:
:param amp_prior:
Prior on log10_h. Default is "log-uniform."
Use "uniform" for upper limits, or "None" to search over
log10_dist instead.
:param dist_prior:
Prior on log10_dist. Default is "None," meaning that the
search is over log10_h instead of log10_dist. Use "log-uniform"
to search over log10_h with a log-uniform prior.
:param skyloc:
Fixed sky location of CW signal search as [cos(theta), phi].
Search over sky location if ``None`` given.
:param log10_fgw:
Fixed log10 GW frequency of CW signal search.
Search over GW frequency if ``None`` given.
:param ecc:
Fixed log10 distance to SMBHB search.
Search over distance or strain if ``None`` given.
:param psrTerm:
Boolean for whether to include the pulsar term. Default is False.
:param name:
Name of CW signal.
"""
if dist_prior is None:
log10_dist = None
if amp_prior == 'uniform':
log10_h = parameter.LinearExp(-18.0, -11.0)('{}_log10_h'.format(name))
elif amp_prior == 'log-uniform':
log10_h = parameter.Uniform(-18.0, -11.0)('{}_log10_h'.format(name))
elif dist_prior == 'log-uniform':
log10_dist = parameter.Uniform(-2.0, 4.0)('{}_log10_dL'.format(name))
log10_h = None
# chirp mass [Msol]
log10_Mc = parameter.Uniform(6.0, 10.0)('{}_log10_Mc'.format(name))
# GW frequency [Hz]
if log10_fgw is None:
log10_fgw = parameter.Uniform(-9.0, -7.0)('{}_log10_fgw'.format(name))
else:
log10_fgw = parameter.Constant(log10_fgw)('{}_log10_fgw'.format(name))
# orbital inclination angle [radians]
cosinc = parameter.Uniform(-1.0, 1.0)('{}_cosinc'.format(name))
# initial GW phase [radians]
phase0 = parameter.Uniform(0.0, 2*np.pi)('{}_phase0'.format(name))
# polarization
psi_name = '{}_psi'.format(name)
psi = parameter.Uniform(0, np.pi)(psi_name)
# sky location
costh_name = '{}_costheta'.format(name)
phi_name = '{}_phi'.format(name)
if skyloc is None:
costh = parameter.Uniform(-1, 1)(costh_name)
phi = parameter.Uniform(0, 2*np.pi)(phi_name)
else:
costh = parameter.Constant(skyloc[0])(costh_name)
phi = parameter.Constant(skyloc[1])(phi_name)
if psrTerm:
# orbital phase
p_phase = parameter.Uniform(0, np.pi)
p_dist = parameter.Normal(0, 1)
else:
p_phase = None
p_dist = 0
# continuous wave signal
wf = cw_delay(cos_gwtheta=costh, gwphi=phi, cos_inc=cosinc,
log10_mc=log10_Mc, log10_fgw=log10_fgw,
log10_h=log10_h, log10_dist=log10_dist,
phase0=phase0, psi=psi,
psrTerm=True, p_dist=p_dist, p_phase=p_phase,
phase_approx=True, check=False,
tref=tref)
cw = CWSignal(wf, ecc=False, psrTerm=psrTerm)
return cw
[docs]def cw_block_ecc(amp_prior='log-uniform', skyloc=None, log10_F=None,
ecc=None, psrTerm=False, tref=0, name='cw'):
"""
Returns deterministic, eccentric orbit continuous GW model:
:param amp_prior:
Prior on log10_h and log10_Mc/log10_dL. Default is "log-uniform" with
log10_Mc and log10_dL searched over. Use "uniform" for upper limits,
log10_h searched over.
:param skyloc:
Fixed sky location of CW signal search as [cos(theta), phi].
Search over sky location if ``None`` given.
:param log10_F:
Fixed log-10 orbital frequency of CW signal search.
Search over orbital frequency if ``None`` given.
:param ecc:
Fixed eccentricity of SMBHB search.
Search over eccentricity if ``None`` given.
:param psrTerm:
Boolean for whether to include the pulsar term. Default is False.
:param name:
Name of CW signal.
"""
if amp_prior == 'uniform':
log10_h = parameter.LinearExp(-18.0, -11.0)('{}_log10_h'.format(name))
elif amp_prior == 'log-uniform':
log10_h = None
# chirp mass [Msol]
log10_Mc = parameter.Uniform(6.0, 10.0)('{}_log10_Mc'.format(name))
# luminosity distance [Mpc]
log10_dL = parameter.Uniform(-2.0, 4.0)('{}_log10_dL'.format(name))
# orbital frequency [Hz]
if log10_F is None:
log10_Forb = parameter.Uniform(-9.0, -7.0)('{}_log10_Forb'.format(name))
else:
log10_Forb = parameter.Constant(log10_F)('{}_log10_Forb'.format(name))
# orbital inclination angle [radians]
cosinc = parameter.Uniform(-1.0, 1.0)('{}_cosinc'.format(name))
# periapsis position angle [radians]
gamma_0 = parameter.Uniform(0.0, np.pi)('{}_gamma0'.format(name))
# Earth-term eccentricity
if ecc is None:
e_0 = parameter.Uniform(0.0, 0.99)('{}_e0'.format(name))
else:
e_0 = parameter.Constant(ecc)('{}_e0'.format(name))
# initial mean anomaly [radians]
l_0 = parameter.Uniform(0.0, 2.0*np.pi)('{}_l0'.format(name))
# mass ratio = M_2/M_1
q = parameter.Constant(1.0)('{}_q'.format(name))
# polarization
pol_name = '{}_pol'.format(name)
pol = parameter.Uniform(0, np.pi)(pol_name)
# sky location
costh_name = '{}_costheta'.format(name)
phi_name = '{}_phi'.format(name)
if skyloc is None:
costh = parameter.Uniform(-1, 1)(costh_name)
phi = parameter.Uniform(0, 2*np.pi)(phi_name)
else:
costh = parameter.Constant(skyloc[0])(costh_name)
phi = parameter.Constant(skyloc[1])(phi_name)
# continuous wave signal
wf = compute_eccentric_residuals(cos_gwtheta=costh, gwphi=phi,
log10_mc=log10_Mc, log10_dist=log10_dL,
log10_h=log10_h, log10_F=log10_Forb,
cos_inc=cosinc, psi=pol, gamma0=gamma_0,
e0=e_0, l0=l_0, q=q, nmax=400,
pdist=None, pphase=None, pgam=None,
tref=tref, check=False)
cw = CWSignal(wf, ecc=True, psrTerm=psrTerm)
return cw
@signal_base.function
def cw_delay(toas, pos, pdist,
cos_gwtheta=0, gwphi=0, cos_inc=0,
log10_mc=9, log10_fgw=-8, log10_dist=None, log10_h=None,
phase0=0, psi=0,
psrTerm=False, p_dist=1, p_phase=None,
evolve=False, phase_approx=False, check=False,
tref=0):
"""
Function to create GW incuced residuals from a SMBMB as
defined in Ellis et. al 2012,2013.
:param toas:
Pular toas in seconds
:param pos:
Unit vector from the Earth to the pulsar
:param pdist:
Pulsar distance (mean and uncertainty) [kpc]
:param cos_gwtheta:
Cosine of Polar angle of GW source in celestial coords [radians]
:param gwphi:
Azimuthal angle of GW source in celestial coords [radians]
:param cos_inc:
cosine of Inclination of GW source [radians]
:param log10_mc:
log10 of Chirp mass of SMBMB [solar masses]
:param log10_fgw:
log10 of Frequency of GW (twice the orbital frequency) [Hz]
:param log10_dist:
log10 of Luminosity distance to SMBMB [Mpc],
used to compute strain, if not None
:param log10_h:
log10 of GW strain,
used to compute distance, if not None
:param phase0:
Initial GW phase of source [radians]
:param psi:
Polarization angle of GW source [radians]
:param psrTerm:
Option to include pulsar term [boolean]
:param p_dist:
Pulsar distance parameter
:param p_phase:
Use pulsar phase to determine distance [radian]
:param evolve:
Option to include/exclude full evolution [boolean]
:param phase_approx:
Option to include/exclude phase evolution across observation time
[boolean]
:param check:
Check if frequency evolves significantly over obs. time [boolean]
:param tref:
Reference time for phase and frequency [s]
:return: Vector of induced residuals
"""
# convert units to time
mc = 10**log10_mc * const.Tsun
fgw = 10**log10_fgw
gwtheta = np.arccos(cos_gwtheta)
inc = np.arccos(cos_inc)
p_dist = (pdist[0] + pdist[1]*p_dist)*const.kpc/const.c
if log10_h is None and log10_dist is None:
raise ValueError("one of log10_dist or log10_h must be non-None")
elif log10_h is not None and log10_dist is not None:
raise ValueError("only one of log10_dist or log10_h can be non-None")
elif log10_h is None:
dist = 10**log10_dist * const.Mpc / const.c
else:
dist = 2 * mc**(5/3) * (np.pi*fgw)**(2/3) / 10**log10_h
if check:
# check that frequency is not evolving significantly over obs. time
fstart = fgw * (1 - 256/5 * mc**(5/3) * fgw**(8/3) * toas[0])**(-3/8)
fend = fgw * (1 - 256/5 * mc**(5/3) * fgw**(8/3) * toas[-1])**(-3/8)
df = fend - fstart
# observation time
Tobs = toas.max()-toas.min()
fbin = 1/Tobs
if np.abs(df) > fbin:
print('WARNING: Frequency is evolving over more than one '
'frequency bin.')
print('f0 = {0}, f1 = {1}, df = {2}, fbin = {3}'.format(fstart, fend, df, fbin))
return np.ones(len(toas)) * np.nan
# get antenna pattern funcs and cosMu
# write function to get pos from theta,phi
fplus, fcross, cosMu = utils.create_gw_antenna_pattern(pos, gwtheta, gwphi)
# get pulsar time
toas -= tref
if p_dist > 0:
tp = toas-p_dist*(1-cosMu)
else:
tp = toas
# orbital frequency
w0 = np.pi * fgw
phase0 /= 2 # convert GW to orbital phase
# omegadot = 96/5 * mc**(5/3) * w0**(11/3) # Not currently used in code
# evolution
if evolve:
# calculate time dependent frequency at earth and pulsar
omega = w0 * (1 - 256/5 * mc**(5/3) * w0**(8/3) * toas)**(-3/8)
omega_p = w0 * (1 - 256/5 * mc**(5/3) * w0**(8/3) * tp)**(-3/8)
if p_dist > 0:
omega_p0 = w0 * (1 + 256/5
* mc**(5/3) * w0**(8/3) * p_dist*(1-cosMu))**(-3/8)
else:
omega_p0 = w0
# calculate time dependent phase
phase = phase0 + 1/32/mc**(5/3) * (w0**(-5/3) - omega**(-5/3))
if p_phase is None:
phase_p = phase0 + 1/32/mc**(5/3) * (w0**(-5/3) - omega_p**(-5/3))
else:
phase_p = (phase0 + p_phase
+ 1/32*mc**(-5/3) * (omega_p0**(-5/3) - omega_p**(-5/3)))
elif phase_approx:
# monochromatic
omega = w0
if p_dist > 0:
omega_p = w0 * (1 + 256/5
* mc**(5/3) * w0**(8/3) * p_dist*(1-cosMu))**(-3/8)
else:
omega_p = w0
# phases
phase = phase0 + omega * toas
if p_phase is not None:
phase_p = phase0 + p_phase + omega_p*toas
else:
phase_p = (phase0 + omega_p*toas
+ 1/32/mc**(5/3) * (w0**(-5/3) - omega_p**(-5/3)))
# no evolution
else:
# monochromatic
omega = np.pi*fgw
omega_p = omega
# phases
phase = phase0 + omega * toas
phase_p = phase0 + omega * tp
# define time dependent coefficients
At = -0.5*np.sin(2*phase)*(3+np.cos(2*inc))
Bt = 2*np.cos(2*phase)*np.cos(inc)
At_p = -0.5*np.sin(2*phase_p)*(3+np.cos(2*inc))
Bt_p = 2*np.cos(2*phase_p)*np.cos(inc)
# now define time dependent amplitudes
alpha = mc**(5./3.)/(dist*omega**(1./3.))
alpha_p = mc**(5./3.)/(dist*omega_p**(1./3.))
# define rplus and rcross
rplus = alpha*(-At*np.cos(2*psi)+Bt*np.sin(2*psi))
rcross = alpha*(At*np.sin(2*psi)+Bt*np.cos(2*psi))
rplus_p = alpha_p*(-At_p*np.cos(2*psi)+Bt_p*np.sin(2*psi))
rcross_p = alpha_p*(At_p*np.sin(2*psi)+Bt_p*np.cos(2*psi))
# residuals
if psrTerm:
res = fplus*(rplus_p-rplus)+fcross*(rcross_p-rcross)
else:
res = -fplus*rplus - fcross*rcross
return res
@signal_base.function
def bwm_delay(toas, pos, log10_h=-14.0, cos_gwtheta=0.0, gwphi=0.0, gwpol=0.0, t0=55000,
antenna_pattern_fn=None):
"""
Function that calculates the earth-term gravitational-wave
burst-with-memory signal, as described in:
Seto et al, van haasteren and Levin, phsirkov et al, Cordes and Jenet.
This version uses the F+/Fx polarization modes, as verified with the
Continuous Wave and Anisotropy papers.
:param toas: Time-of-arrival measurements [s]
:param pos: Unit vector from Earth to pulsar
:param log10_h: log10 of GW strain
:param cos_gwtheta: Cosine of GW polar angle
:param gwphi: GW azimuthal polar angle [rad]
:param gwpol: GW polarization angle
:param t0: Burst central time [day]
:param antenna_pattern_fn:
User defined function that takes `pos`, `gwtheta`, `gwphi` as
arguments and returns (fplus, fcross)
:return: the waveform as induced timing residuals (seconds)
"""
# convert
h = 10 ** log10_h
gwtheta = np.arccos(cos_gwtheta)
t0 *= const.day
# antenna patterns
if antenna_pattern_fn is None:
apc = utils.create_gw_antenna_pattern(pos, gwtheta, gwphi)
else:
apc = antenna_pattern_fn(pos, gwtheta, gwphi)
# grab fplus, fcross
fp, fc = apc[0], apc[1]
# combined polarization
pol = np.cos(2 * gwpol) * fp + np.sin(2 * gwpol) * fc
# Return the time-series for the pulsar
return pol * h * np.heaviside(toas - t0, 0.5) * (toas - t0)
@signal_base.function
def bwm_sglpsr_delay(toas, sign, log10_A=-15, t0=55000):
"""
Function that calculates the earth-term gravitational-wave
burst-with-memory signal for an optimally oriented source in a single pulsar
:param toas: Time-of-arrival measurements [s]
:param log10_A: log10 of the amplitude of the ramp (delta_f/f)
:param t0: Burst central time [day]
:return: the waveform as induced timing residuals (seconds)
"""
A = 10 ** log10_A
t0 *= const.day
# Return the time-series for the pulsar
def heaviside(x):
return 0.5 * (np.sign(x) + 1)
# return 0 #Fix the return to 0 in order to test what the heck is wrong with red noise detection in bwm
return A * np.sign(sign) * heaviside(toas - t0) * (toas - t0)
@signal_base.function
def compute_eccentric_residuals(toas, theta, phi, cos_gwtheta, gwphi,
log10_mc, log10_dist, log10_h, log10_F, cos_inc,
psi, gamma0, e0, l0, q, nmax=400, pdist=1.0,
pphase=None, pgam=None, psrTerm=False,
tref=0, check=False):
"""
Simulate GW from eccentric SMBHB. Waveform models from
Taylor et al. (2015) and Barack and Cutler (2004).
WARNING: This residual waveform is only accurate if the
GW frequency is not significantly evolving over the
observation time of the pulsar.
:param toa: pulsar observation times
:param theta: polar coordinate of pulsar
:param phi: azimuthal coordinate of pulsar
:param gwtheta: Polar angle of GW source in celestial coords [radians]
:param gwphi: Azimuthal angle of GW source in celestial coords [radians]
:param log10_mc: Base-10 lof of chirp mass of SMBMB [solar masses]
:param log10_dist: Base-10 uminosity distance to SMBMB [Mpc]
:param log10_F: base-10 orbital frequency of SMBHB [Hz]
:param inc: Inclination of GW source [radians]
:param psi: Polarization of GW source [radians]
:param gamma0: Initial angle of periastron [radians]
:param e0: Initial eccentricity of SMBHB
:param l0: Initial mean anomoly [radians]
:param q: Mass ratio of SMBHB
:param nmax: Number of harmonics to use in waveform decomposition
:param pdist: Pulsar distance [kpc]
:param pphase: Pulsar phase [rad]
:param pgam: Pulsar angle of periastron [rad]
:param psrTerm: Option to include pulsar term [boolean]
:param tref: Fidicuial time at which initial parameters are referenced [s]
:param check: Check if frequency evolves significantly over obs. time
:returns: Vector of induced residuals
"""
# convert from sampling
F = 10.0**log10_F
mc = 10.0**log10_mc
dist = 10.0**log10_dist
if log10_h is not None:
h0 = 10.0**log10_h
else:
h0 = None
inc = np.arccos(cos_inc)
gwtheta = np.arccos(cos_gwtheta)
# define variable for later use
cosgwtheta, cosgwphi = np.cos(gwtheta), np.cos(gwphi)
singwtheta, singwphi = np.sin(gwtheta), np.sin(gwphi)
sin2psi, cos2psi = np.sin(2*psi), np.cos(2*psi)
# unit vectors to GW source
m = np.array([singwphi, -cosgwphi, 0.0])
n = np.array([-cosgwtheta*cosgwphi, -cosgwtheta*singwphi, singwtheta])
omhat = np.array([-singwtheta*cosgwphi, -singwtheta*singwphi, -cosgwtheta])
# pulsar position vector
phat = np.array([np.sin(theta)*np.cos(phi), np.sin(theta)*np.sin(phi),
np.cos(theta)])
fplus = 0.5 * (np.dot(m, phat)**2 - np.dot(n, phat)**2) / (1+np.dot(omhat, phat))
fcross = (np.dot(m, phat)*np.dot(n, phat)) / (1 + np.dot(omhat, phat))
cosMu = -np.dot(omhat, phat)
# get values from pulsar object
toas = toas.copy() - tref
if check:
# check that frequency is not evolving significantly over obs. time
y = utils.solve_coupled_ecc_solution(F, e0, gamma0, l0, mc, q,
np.array([0.0, toas.max()]))
# initial and final values over observation time
Fc0, ec0, gc0, phic0 = y[0, :]
Fc1, ec1, gc1, phic1 = y[-1, :]
# observation time
Tobs = 1/(toas.max()-toas.min())
if np.abs(Fc0-Fc1) > 1/Tobs:
print('WARNING: Frequency is evolving over more than one frequency bin.')
print('F0 = {0}, F1 = {1}, delta f = {2}'.format(Fc0, Fc1, 1/Tobs))
return np.ones(len(toas)) * np.nan
# get gammadot for earth term
gammadot = utils.get_gammadot(F, mc, q, e0)
# get number of harmonics to use
if not isinstance(nmax, int):
if e0 < 0.999 and e0 > 0.001:
nharm = int(nmax(e0))
elif e0 < 0.001:
nharm = 2
else:
nharm = int(nmax(0.999))
else:
nharm = nmax
# no more than 100 harmonics
nharm = min(nharm, 100)
##### earth term #####
splus, scross = utils.calculate_splus_scross(nmax=nharm, mc=mc, dl=dist,
h0=h0, F=F, e=e0, t=toas.copy(),
l0=l0, gamma=gamma0,
gammadot=gammadot, inc=inc)
##### pulsar term #####
if psrTerm:
# pulsar distance
pd = pdist
# convert units
pd *= const.kpc / const.c
# get pulsar time
tp = toas.copy() - pd * (1-cosMu)
# solve coupled system of equations to get pulsar term values
y = utils.solve_coupled_ecc_solution(F, e0, gamma0, l0, mc,
q, np.array([0.0, tp.min()]))
# get pulsar term values
if np.any(y):
Fp, ep, gp, phip = y[-1, :]
# get gammadot at pulsar term
gammadotp = utils.get_gammadot(Fp, mc, q, ep)
# get phase at pulsar
if pphase is None:
lp = phip
else:
lp = pphase
# get angle of periastron at pulsar
if pgam is None:
gp = gp
else:
gp = pgam
# get number of harmonics to use
if not isinstance(nmax, int):
if e0 < 0.999 and e0 > 0.001:
nharm = int(nmax(e0))
elif e0 < 0.001:
nharm = 2
else:
nharm = int(nmax(0.999))
else:
nharm = nmax
# no more than 1000 harmonics
nharm = min(nharm, 100)
splusp, scrossp = utils.calculate_splus_scross(nmax=nharm, mc=mc,
dl=dist, h0=h0,
F=Fp, e=ep,
t=toas.copy(),
l0=lp, gamma=gp,
gammadot=gammadotp,
inc=inc)
rr = (fplus*cos2psi - fcross*sin2psi) * (splusp - splus) + \
(fplus*sin2psi + fcross*cos2psi) * (scrossp - scross)
else:
rr = np.ones(len(toas)) * np.nan
else:
rr = - (fplus*cos2psi - fcross*sin2psi) * splus - \
(fplus*sin2psi + fcross*cos2psi) * scross
return rr
[docs]def CWSignal(cw_wf, ecc=False, psrTerm=False, name='cw'):
BaseClass = deterministic_signals.Deterministic(cw_wf, name=name)
class CWSignal(BaseClass):
def __init__(self, psr):
super(CWSignal, self).__init__(psr)
self._wf[''].add_kwarg(psrTerm=psrTerm)
if ecc:
pgam = parameter.Uniform(0, 2*np.pi)('_'.join([psr.name,
'pgam',
name]))
self._params['pgam'] = pgam
self._wf['']._params['pgam'] = pgam
return CWSignal
@signal_base.function
def generalized_gwpol_psd(f, log10_A_tt=-15, log10_A_st=-15,
log10_A_vl=-15, log10_A_sl=-15,
kappa=10/3, p_dist=1.0):
"""
PSD for a generalized mixture of scalar+vector dipole radiation
and tensorial quadrupole radiation from SMBHBs.
"""
df = np.diff(np.concatenate((np.array([0]), f[::2])))
euler_e = 0.5772156649
pdist = p_dist * const.kpc / const.c
orf_aa_tt = (2/3) * np.ones(len(f))
orf_aa_st = (2/3) * np.ones(len(f))
orf_aa_vl = 2*np.log(4*np.pi*f*pdist) - 14/3 + 2*euler_e
orf_aa_sl = np.pi**2*f*pdist/4 - \
np.log(4*np.pi*f*pdist) + 37/24 - euler_e
prefactor = (1 + kappa**2) / (1 + kappa**2 * (f / const.fyr)**(-2/3))
gwpol_amps = 10**(2*np.array([log10_A_tt, log10_A_st,
log10_A_vl, log10_A_sl]))
gwpol_factors = np.array([orf_aa_tt*gwpol_amps[0],
orf_aa_st*gwpol_amps[1],
orf_aa_vl*gwpol_amps[2],
orf_aa_sl*gwpol_amps[3]])
S_psd = prefactor * (gwpol_factors[0, :] * (f / const.fyr)**(-4/3) +
np.sum(gwpol_factors[1:, :], axis=0) *
(f / const.fyr)**(-2)) / \
(8*np.pi**2*f**3)
return S_psd * np.repeat(df, 2)
@signal_base.function
def fdm_delay(toas, log10_A, log10_f, phase_e, phase_p):
"""
Function that calculates the earth-term gravitational-wave
fuzzy dark matter signal, as described in:
Kato et al. (2020).
:param toas: Time-of-arrival measurements [s]
:param log10_A: log10 of GW strain
:param log10_f: log10 of GW frequency
:param phase_e: The Earth-term phase of the GW
:param phase_p: The Pulsar-term phase of the GW
:return: the waveform as induced timing residuals (seconds)
"""
# convert
A = 10 ** log10_A
f = 10 ** log10_f
# Return the time-series for the pulsar
return - A / (2 * np.pi * f) * (np.sin(2 * np.pi * f * toas + phase_e) - np.sin(2 * np.pi * f * toas + phase_p))