diff --git a/.gitignore b/.gitignore index 84229f4..748488d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,8 +1,12 @@ +#MacBook stuff +.DS_Store +.gitignore + # Byte-compiled / optimized / DLL files __pycache__/ *.py[cod] *$py.class - +**/.DS_Store # C extensions *.so diff --git a/hasasia/sensitivity.py b/hasasia/sensitivity.py index b60378f..230e093 100644 --- a/hasasia/sensitivity.py +++ b/hasasia/sensitivity.py @@ -265,6 +265,7 @@ def get_NcalInv(psr, nf=200, fmin=None, fmax=2e-7, freqs=None, NcalInv = np.linalg.inv(Ncal) #N_TOA-N_par x N_TOA-N_par TfN = np.matmul(np.conjugate(Gtilde),np.matmul(NcalInv,Gtilde.T)) / 2 + # N_freq x N_freq if return_Gtilde_Ncal: return np.real(TfN), Gtilde, Ncal elif full_matrix: @@ -315,13 +316,17 @@ class Pulsar(object): Earth-pulsar distance. Default units is kpc. """ - def __init__(self, toas, toaerrs, phi=None, theta=None, - designmatrix=None, N=None, pdist=1.0*u.kpc): + def __init__(self, toas, toaerrs, phi=None, theta=None, name=None, + designmatrix=None, N=None, pdist=1.0*u.kpc, A_rn=None, + alpha=None,): self.toas = toas self.toaerrs = toaerrs self.phi = phi self.theta = theta + self.name = name self.pdist = make_quant(pdist,'kpc') + self.A_rn = A_rn + self.alpha = alpha if N is None: self.N = np.diag(toaerrs**2) #N ==> weights @@ -334,6 +339,175 @@ def __init__(self, toas, toaerrs, phi=None, theta=None, else: self.designmatrix = designmatrix + def filter_data(self, start_time=None, end_time=None): + """ + Parameters + ========== + start_time - float + MJD at which to begin data subset. + end_time - float + MJD at which to end data subset. + + Filter data to create a time-slice of overall dataset. + Function adapted from enterprise.BasePulsar() class. + """ + if start_time is None and end_time is None: + mask = np.ones(self.toas.shape, dtype=bool) + else: + mask = np.logical_and(self.toas >= start_time * 86400, self.toas <= end_time * 86400) + + self.toas = self.toas[mask] + self.toaerrs = self.toaerrs[mask] + self.N = self.N[mask, :][:, mask] + + self.designmatrix = create_design_matrix(self.toas, RADEC=True, PROPER=True, PX=True) + #self.designmatrix = self.designmatrix[mask, :] + #dmx_mask = np.sum(self.designmatrix, axis=0) != 0.0 + #self.designmatrix = self.designmatrix[:, dmx_mask] + self._G = G_matrix(designmatrix=self.designmatrix) + + def change_cadence(self, start_time=0, end_time=1_000_000, + cadence=None, cadence_factor=4, uneven=False, + A_gwb=None, alpha_gwb=-2/3., freqs=None, + fast=True,): + """ + Parameters + ========== + start_time - float + MJD at which to begin altared cadence. + end_time - float + MJD at which to end altared cadence. + cadence - float + cadence for the modified campaign [toas/year] + cadence_facter - float + (instead of cadence) factor by which to modify the old cadence. + uneven - bool + whether or not to evenly space observation epochs + A_gwb - float + amplitude of injected gwb self-noise + alpha_gwb - float + spectral index of injected gwb self-noise. + note that this is residual space spectral index. + freqs - array + frequencies to construct the gwb noise and intrinsic noise + fast - bool + faster but slightly less accurate method to calculate noise injected in N. + + Change observing cadence in a given time range. + Recalculate pulsar noise properties. + """ + mask_before = self.toas <= start_time * 86400 + mask_after = self.toas >= end_time * 86400 + old_Ntoas = np.sum( + np.logical_and(self.toas >= start_time * 86400, + self.toas <= end_time * 86400) + ) + # store the old toas and errors + old_toas = self.toas + old_toaerrs = self.toaerrs + # calculate old cadence then modified cadence + if start_time < min(old_toas)/84600: + start_time = min(old_toas)/84600 + if end_time < min(old_toas)/84600: + print("trying to change non-existant campaign") + return 0 + duration = end_time - start_time # in MJD + old_cadence = old_Ntoas / duration * 365.25 # cad is Ntoas/year + if cadence is not None: + new_cadence = cadence + else: + new_cadence = old_cadence * cadence_factor + # create new toas and toa errors + campaign_Ntoas = int(np.floor( duration / 365.25 * new_cadence )) + campaign_toas = np.linspace(start_time, end_time, campaign_Ntoas) * 86400 + if uneven: + # FIXME check this with jeff to see what he was going for + # in sim_pta() + dt = duration / campaign_Ntoas / 8 * yr_sec + campaign_toas += np.random.uniform(-dt, dt, size=campaign_Ntoas) + self.toas = np.concatenate([old_toas[mask_before], campaign_toas, old_toas[mask_after]]) + campaign_toaerrs = np.median(old_toaerrs)*np.ones(campaign_Ntoas) + # FIXME can only use a fixed toaerr for the duration of the campaign + #self.toaerrs = np.concatenate([old_toaerrs[mask_before], campaign_toaerrs, old_toaerrs[mask_after]]) + self.toaerrs = np.ones(len(self.toas))*old_toaerrs[0] + print(f"old: {len(old_toaerrs)}, new: {len(self.toaerrs)}") + # recalculate N, designmatrix, G with new toas + N = np.diag(self.toaerrs**2) + if self.A_rn is not None: + plaw = red_noise_powerlaw(A=self.A_rn, + alpha=self.alpha, + freqs=freqs) + N += corr_from_psd(freqs=freqs, psd=plaw, toas=self.toas, fast=fast) + + if A_gwb is not None: + gwb = red_noise_powerlaw(A=A_gwb, + alpha=alpha_gwb, + freqs=freqs) + N += corr_from_psd(freqs=freqs, psd=gwb, toas=self.toas, fast=fast) + self.designmatrix = create_design_matrix(self.toas, RADEC=True, PROPER=True, PX=True) + self._G = G_matrix(designmatrix=self.designmatrix) + self.N = N + + def change_sigma(self, start_time=0, end_time=1_000_000, + new_sigma=None, sigma_factor=4, uneven=False, + A_gwb=None, alpha_gwb=-2/3., freqs=None, + fast=True,): + """ + Parameters + ========== + start_time - float + MJD at which to begin altared cadence. + end_time - float + MJD at which to end altared cadence. + new_sigma - float + unvertainty of toas for the modified campaign [microseconds] + sigma_facter - float + (instead of sigmas) factor by which to modify the campaign sigmas. + uneven - bool + whether or not to evenly space observation epochs + A_gwb - float + amplitude of injected gwb self-noise + alpha_gwb - float + spectral index of injected gwb self-noise. + note that this is residual space spectral index (alpha). + freqs - array + frequencies to construct the gwb noise and intrinsic noise + fast - bool + faster but slightly less accurate method to calculate noise injected in N. + + Change observing cadence in a given time range. + Recalculate pulsar noise properties. + """ + mask_before = self.toas <= start_time * 86400 + mask_after = self.toas >= end_time * 86400 + campaign_mask = np.logical_and(self.toas >= start_time * 86400, + self.toas <= end_time * 86400) + campaign_ntoas = np.sum(campaign_mask) + # store the old toa errors + toaerrs_campaign = self.toaerrs[campaign_mask] + # modify the campaign toa errors + if sigma_factor is not None: + toaerrs_campaign = sigma_factor * toaerrs_campaign + elif sigma_factor is None and new_sigma is not None: + toaerrs_campaign = np.ones(campaign_ntoas)*new_sigma + self.toaerrs = np.concatenate([self.toaerrs[mask_before], toaerrs_campaign, self.toaerrs[mask_after]]) + # recalculate N, designmatrix, G with new toas + N = np.diag(self.toaerrs**2) + if self.A_rn is not None: + plaw = red_noise_powerlaw(A=self.A_rn, + alpha=self.alpha, + freqs=freqs) + N += corr_from_psd(freqs=freqs, psd=plaw, toas=self.toas, fast=fast) + + if A_gwb is not None: + gwb = red_noise_powerlaw(A=A_gwb, + alpha=alpha_gwb, + freqs=freqs) + N += corr_from_psd(freqs=freqs, psd=gwb, toas=self.toas, fast=fast) + self.designmatrix = create_design_matrix(self.toas, RADEC=True, PROPER=True, PX=True) + self._G = G_matrix(designmatrix=self.designmatrix) + self.N = N + @property def G(self): """Inverse Noise Weighted Transmission Function.""" @@ -672,6 +846,98 @@ def S_effIJ(self): return self._S_effIJ + def dmu_dbeta_pl_gwb(self, A_gwb, gamma_gwb, components=20): + """Derivative of the GWB signal template (mu) w.r.t. the signal parameters (beta). + $$\mu_{\rm GWB} = \frac{A_{\rm GWB}}{12\pi^2} * \left(\frac{f}{f_{\rm ref}}\right)^{-2\alpha} \rm yr^3""" + # FIXME : what units should this be in ?? + fbins = np.linspace(1/self.Tspan, components/self.Tspan, components) + #fbins = self.freqs + # dmu_dA = A_gwb/(12*np.pi**2)*(fbins/fyr)**(-2*gamma_gwb) * \ + # yr_sec**3 * np.log(10) + # dmu_dalpha = A_gwb/(12*np.pi**2)*(fbins/fyr)**(-2*gamma_gwb) * \ + # yr_sec**3 * np.log(fbins/fyr) + dmu_dA = A_gwb/(12*np.pi**2)*(fbins/fyr)**(-1*gamma_gwb) * \ + yr_sec**3 * np.log(10) + dmu_dalpha = A_gwb/(12*np.pi**2)*(fbins/fyr)**(-1*gamma_gwb) * \ + yr_sec**3 * np.log(fbins/fyr) + + return np.array([dmu_dA, dmu_dalpha]) + + def dmu_dbeta_freespec_gwb(self, log10_gw_rhos): + """Derivative of the GWB signal template (mu) w.r.t. the signal parameters (beta). + $$\mu_{\rm GWB} = \rho_{i}^2""" + dmu_drhos = [] + for i, _ in enumerate(log10_gw_rhos): + dmu_drho = np.zeros(len(log10_gw_rhos)) + # https://www.wolframalpha.com/input?i=derivative+of+10%5Ex%5E2+wrt+x + dmu_drho[i]+=10**log10_gw_rhos**2 * np.log(10) * 2*log10_gw_rhos + dmu_drhos.append(dmu_drho) + + return np.array(dmu_drhos) + + def GWBFisherMatrix(self, + Agwb, gamma_gwb, + log10_gw_rhos, + psd_type='powerlaw', + components=20): + """ + GWB Fisher Matrix + Fisher matrix = signal derivative x noise covariance x signal derivative + """ + if psd_type == 'powerlaw': + assert(Agwb is not None) + assert(gamma_gwb is not None) + dmu_dbeta = self.dmu_dbeta_pl_gwb(Agwb, gamma_gwb, components=components) + elif psd_type == 'freespec': + assert(log10_gw_rhos is not None) + dmu_dbeta = self.dmu_dbeta_freespec_gwb(log10_gw_rhos) + else: + raise NotImplementedError(f"{psd_type} PSD not supported yet! \ + Choose \'powerlaw\' or \'freespec\'.") + fbins=np.linspace(1/self.Tspan, components/self.Tspan, components) + #fbins=self.freqs + fidxs = np.array([np.argmin(abs(self.freqs - f)) for f in fbins]) + # FIXME: check the order of opperations here + #ncalinv = get_NcalInv(psr, freqs=self.freqs[fidxs]) + print("seff: ", np.diag(self.S_eff[fidxs]).shape) + print("dmu: ", dmu_dbeta.shape) + #print("ncalinv: ", psr.NcalInv.shape) + #first = np.matmul(np.diag(psr.NcalInv[fidxs]), dmu_dbeta.T) + #first = np.matmul(np.diag(self.S_eff[fidxs]), dmu_dbeta.T) + # sigmainv = np.diag(self.S_eff[fidxs]) + # first = np.matmul(sigmainv, dmu_dbeta.T) + # print("1st matmul", first.shape) + # mathcalF = np.matmul(dmu_dbeta, first) + # print("FishMat", mathcalF.shape) + fishers = np.zeros((2,2)) + for i in range(self.S_effIJ.shape[0]): + fisher = np.matmul(dmu_dbeta, np.matmul(np.diag(self.S_effIJ[i][fidxs]), dmu_dbeta.T)) + fishers = fisher + fishers + + return fishers + + def GWBFisherUncertainty(self, + Agwb=None, gamma_gwb=None, + log10_gw_rhos=None, + psd_type='powerlaw', + components=20): + """ + GWB Fisher Matrix Uncertainty + uncertainty of the i'th param is the square root of the the + (i,i) element of the inverse Fisher matrix + """ + + mathcalF = self.GWBFisherMatrix(Agwb, gamma_gwb, + log10_gw_rhos, + psd_type=psd_type, + components=components) + + mathcalFinv = np.linalg.inv(mathcalF) + print("FishMatInv", mathcalF.shape) + #mathcalFinv = mathcalF + + return np.sqrt(np.diag(mathcalFinv)) + class DeterSensitivityCurve(SensitivityCurve): ''' @@ -1022,9 +1288,13 @@ def get_Tspan(psrs): """ Returns the total timespan from a list or arry of Pulsar objects, psrs. """ - last = np.amax([p.toas.max() for p in psrs]) - first = np.amin([p.toas.min() for p in psrs]) - return last - first + try: + last = np.amax([p.toas.max() for p in psrs]) + first = np.amin([p.toas.min() for p in psrs]) + tspan = last-first + except ValueError: + tspan = 0 + return tspan def get_TspanIJ(psr1,psr2): """ @@ -1189,6 +1459,26 @@ def red_noise_powerlaw(A, freqs, gamma=None, alpha=None): return A**2*(freqs/fyr)**(-gamma)/(12*np.pi**2) * yr_sec**3 +def psd_from_background_realization(background_hc, freqs): + r""" + Calculate the power spectral density with given background strain and frequency binning. + + Parameters + ---------- + background_hc : array + Characteristic strain (hc) of background at each frequency. + + freqs : array + Frequency bins over which the background is stored. + + Returns + ------- + S_h : array + the power spectral density of the background + """ + + return background_hc**2 / (12 * np.pi**2 * freqs[:,np.newaxis]**3) + def S_h(A, alpha, freqs): r""" Add power law red noise to the prefit residual power spectral density. diff --git a/hasasia/sim.py b/hasasia/sim.py index 2ad5d54..a6bb561 100644 --- a/hasasia/sim.py +++ b/hasasia/sim.py @@ -11,8 +11,8 @@ yr_sec = 365.25*24*3600 def sim_pta(timespan, cad, sigma, phi, theta, Npsrs=None, - A_rn=None, alpha=None, freqs=None, uneven=False, A_gwb=None, - fast=True, + A_rn=None, alpha=None, freqs=None, uneven=False, + A_gwb=None, alpha_gwb = -2/3, fast=True, psr_names = None, kwastro={'RADEC':True, 'PROPER':True, 'PX':True}): """ Make a simulated pulsar timing array. Using the available parameters, @@ -98,6 +98,8 @@ def sim_pta(timespan, cad, sigma, phi, theta, Npsrs=None, psrs = [] err_dim = pars['sigma'].ndim Timespan = np.amax(pars['timespan']) + if psr_names is None: + psr_names = [None for _ in range(Npsrs)] for ii in range(Npsrs): tspan = pars['timespan'][ii] Ntoas = int(np.floor(tspan*pars['cad'][ii])) @@ -121,14 +123,22 @@ def sim_pta(timespan, cad, sigma, phi, theta, Npsrs=None, if A_gwb is not None: gwb = red_noise_powerlaw(A=A_gwb, - alpha=-2/3., + alpha=alpha_gwb, freqs=freqs) N += corr_from_psd(freqs=freqs, psd=gwb, toas=toas, fast=fast) M = create_design_matrix(toas, **kwastro) - - p = Pulsar(toas, toaerrs, phi=pars['phi'][ii], - theta=pars['theta'][ii], designmatrix=M, N=N) + # FIXME to always loop over A_rn + if 'A_rn' in keys: + A_rn = pars['A_rn'][ii] + alpha = pars['alpha'][ii] + else: + A_rn = None + + p = Pulsar(toas, toaerrs, name=psr_names[ii], + phi=pars['phi'][ii], theta=pars['theta'][ii], + A_rn=A_rn, alpha=alpha, + designmatrix=M, N=N) psrs.append(p) return psrs diff --git a/hasasia/skymap.py b/hasasia/skymap.py index 8b42e00..0e26444 100644 --- a/hasasia/skymap.py +++ b/hasasia/skymap.py @@ -3,9 +3,11 @@ """Main module.""" import numpy as np import scipy.special as spec +import healpy as hp import astropy.units as u import astropy.constants as c from .sensitivity import DeterSensitivityCurve, resid_response, get_dt +from .utils import strain_and_chirp_mass_to_luminosity_distance __all__ = ['SkySensitivity', 'h_circ', @@ -114,8 +116,10 @@ def SNR(self, h0, iota=None, psi=None): .. _[1]: https://arxiv.org/abs/1907.04341 ''' - if iota is None or psi is None: + if iota is None and psi is None: S_eff = self.S_eff + elif psi is not None and iota is None: + raise NotImplementedError('Currently cannot marginalize over phase but not inclination.') else: S_eff = self.S_eff_full(iota, psi) return h0 * np.sqrt(self.Tspan / S_eff) @@ -144,7 +148,7 @@ def h_thresh(self, SNR=1, iota=None, psi=None): An array representing the skymap of amplitudes needed to see the given signal with the SNR threshold specified. ''' - if iota is None or psi is None: + if iota is None: S_eff = self.S_eff else: S_eff = self.S_eff_full(iota, psi) @@ -174,39 +178,49 @@ def A_gwb(self, h_div_A, SNR=1): integrand = h_div_A[:,np.newaxis]**2 / self.S_eff return SNR / np.sqrt(np.trapz(integrand,x=self.freqs,axis=0 )) - def sky_response_full(self, iota, psi): + def sky_response_full(self, iota, psi=None): r''' Calculate the signal-to-noise ratio of a source given the strain amplitude. This is based on Equation (79) from Hazboun, et al., 2019 - `[1]`_. + `[1]`_, but modified so that you calculate it for a particular inclination angle. .. math:: \rho(\hat{n})=h_0\sqrt{\frac{T_{\rm obs}}{S_{\rm eff}(f_0 ,\hat{k})}} .. _[1]: https://arxiv.org/abs/1907.04341 ''' + if iota is None and psi is not None: + raise NotImplementedError('Currently cannot marginalize over inclination but not phase.') Ap_sqr = (0.5 * (1 + np.cos(iota)**2))**2 Ac_sqr = (np.cos(iota))**2 - iota = iota if isinstance(iota, (int,float)) else np.array(iota) - psi = psi if isinstance(psi, (int,float)) else np.array(psi) - spsi = np.sin(2*np.array(psi)) - cpsi = np.cos(2*np.array(psi)) - if isinstance(Ap_sqr,np.ndarray) or isinstance(spsi,np.ndarray): - c1 = Ac_sqr[:,np.newaxis]*spsi**2 + Ap_sqr[:,np.newaxis]*cpsi**2 - c2 = (Ap_sqr[:,np.newaxis] - Ac_sqr[:,np.newaxis]) * cpsi * spsi - c3 = Ap_sqr[:,np.newaxis]*spsi**2 + Ac_sqr[:,np.newaxis]*cpsi**2 - term1 = self.Fplus[:,:,np.newaxis,np.newaxis]**2 * c1 - term2 = 2 * self.Fplus[:,:,np.newaxis,np.newaxis] * self.Fcross[:,:,np.newaxis,np.newaxis] * c2 - term3 = self.Fcross[:,:,np.newaxis,np.newaxis]**2 * c3 - elif isinstance(Ap_sqr,(int,float)) or isinstance(spsi,(int,float)): - c1 = Ac_sqr*spsi**2 + Ap_sqr*cpsi**2 - c2 = (Ap_sqr - Ac_sqr) * cpsi * spsi - c3 = Ap_sqr*spsi**2 + Ac_sqr*cpsi**2 - term1 = self.Fplus**2 * c1 - term2 = 2 * self.Fplus * self.Fcross * c2 - term3 = self.Fcross**2 * c3 - - return term1 + term2 + term3 + if psi is None: # case where we average over polarization but not inclination + # 0.5 is from averaging over polarization + # Fplus*Fcross term goes to zero + if isinstance(Ap_sqr,np.ndarray): + return (self.Fplus[:,:,np.newaxis]**2 + self.Fcross[:,:,np.newaxis]**2) * 0.5 * (Ap_sqr + Ac_sqr) + elif isinstance(Ap_sqr,(int,float)): + return (self.Fplus**2 + self.Fcross**2) * 0.5 * (Ac_sqr + Ap_sqr) + else: # case where we don't average over polarization or inclination + iota = iota if isinstance(iota, (int,float)) else np.array(iota) + psi = psi if isinstance(psi, (int,float)) else np.array(psi) + spsi = np.sin(2*np.array(psi)) + cpsi = np.cos(2*np.array(psi)) + if isinstance(Ap_sqr,np.ndarray) or isinstance(spsi,np.ndarray): + c1 = Ac_sqr[:,np.newaxis]*spsi**2 + Ap_sqr[:,np.newaxis]*cpsi**2 + c2 = (Ap_sqr[:,np.newaxis] - Ac_sqr[:,np.newaxis]) * cpsi * spsi + c3 = Ap_sqr[:,np.newaxis]*spsi**2 + Ac_sqr[:,np.newaxis]*cpsi**2 + term1 = self.Fplus[:,:,np.newaxis,np.newaxis]**2 * c1 + term2 = 2 * self.Fplus[:,:,np.newaxis,np.newaxis] * self.Fcross[:,:,np.newaxis,np.newaxis] * c2 + term3 = self.Fcross[:,:,np.newaxis,np.newaxis]**2 * c3 + elif isinstance(Ap_sqr,(int,float)) or isinstance(spsi,(int,float)): + c1 = Ac_sqr*spsi**2 + Ap_sqr*cpsi**2 + c2 = (Ap_sqr - Ac_sqr) * cpsi * spsi + c3 = Ap_sqr*spsi**2 + Ac_sqr*cpsi**2 + term1 = self.Fplus**2 * c1 + term2 = 2 * self.Fplus * self.Fcross * c2 + term3 = self.Fcross**2 * c3 + + return term1 + term2 + term3 def S_SkyI_full(self, iota, psi): """Per Pulsar Strain power sensitivity. """ @@ -221,14 +235,22 @@ def S_SkyI_full(self, iota, psi): if sky_resp.ndim == 2: self._S_SkyI_full = (RNcalInv.T[:,:,np.newaxis] * sky_resp[np.newaxis,:,:]) - if sky_resp.ndim == 4: + elif sky_resp.ndim == 3: + self._S_SkyI_full = (RNcalInv.T[:,:,np.newaxis,np.newaxis] + * sky_resp[np.newaxis,:,:,:]) + elif sky_resp.ndim == 4: self._S_SkyI_full = (RNcalInv.T[:,:,np.newaxis,np.newaxis,np.newaxis] * sky_resp[np.newaxis,:,:,:,:]) return self._S_SkyI_full def S_eff_full(self, iota, psi): - """Strain power sensitivity. """ + """ + Strain power sensitivity. + + Can calculate margninalized over polarization or not + with inclination explicit in both cases. + """ if self.pulsar_term == 'explicit': raise NotImplementedError('Currently cannot use pulsar term.') # self._S_eff_full = 1.0 / (4./5 * np.sum(self.S_SkyI, axis=1)) @@ -255,7 +277,8 @@ def S_eff(self): @property def S_SkyI(self): - """Per Pulsar Strain power sensitivity. """ + """Per Pulsar Strain power sensitivity. + (Technically, 1 over this is the per pulsar strain power sensitivity.)""" if not hasattr(self, '_S_SkyI'): t_I = self.T_I / self.Tspan RNcalInv = t_I[:,np.newaxis] / self.SnI @@ -290,6 +313,74 @@ def S_eff_mean(self): self._S_eff_mean = 1.0 / (12./5 * mean_sky) return self._S_eff_mean + @property + def dmu_dbeta_cw(self, Mc, D_l, f0, theta_cw, phi_cw, h0, + iota=None, psi=None, phi0=None): + """Derivative of the CGW signal template (mu) w.r.t. the signal parameters (beta). + $$\mu_{\rm CW} = h_+ F^+_I(\hat{k}) + h_x F^x_I(\hat{k})""" + hplus = jj + hcross = jj + h0 = h0_circ(Mc, D_l, f0) + dmu_dbeta = [] + + return np.array(dmu_dbeta) + + @property + def CWFisherMatrix(self, Mc, f0, theta_cw, phi_cw, h0, + iota=None, psi=None, phi0=None): + """Continuous GW Fisher Matrix""" + dmu_dbeta = self.dmu_dbeta_cw(Mc, f0, theta_cw, phi_cw, h0, iota=None, psi=None, phi0=None) + mathcalF = dmu_dbeta.T * self.NcalInv * dmu_dbeta + + return mathcalF + + @property + def CWFisherUncertainty(self, Mc, f0, theta_cw, phi_cw, h0, + iota=None, psi=None, phi0=None): + """GWB Fisher Matrix Uncertainty""" + mathcalF = self.CWFisherMatrix(Mc, f0, theta_cw, phi_cw, h0, + iota=None, psi=None, phi0=None) + mathcalFinv = np.linalg.inv(mathcalF) + + return np.sqrt(np.diag(mathcalFinv)) + +def calculate_detection_volume(self, f0, SNR_threshold=3.7145, M_c=1e9): + r""" + Calculates the detection volume of the PTA + at a given frequency or list of frequencies. + + Parameters + ---------- + + f0 : float + the frequency [Hz] at which to calculate detection volume + + SNR_threshold : float + the signal to noise to referene detection volume to + + M_c : float + the chirp mass [Msun] at which to reference detection volume + + Returns + ------- + + volume : float + the detection volume in Mpc^3 + + """ + NSIDE = hp.pixelfunc.npix2nside(self.S_eff.shape[1]) + dA = hp.pixelfunc.nside2pixarea(NSIDE, degrees=False) + if isinstance(f0, (int,float)): + f_idx = np.array([np.argmin(abs(self.freqs - f0))]) + elif isinstance(f0, (np.ndarray, list)): + f_idx = np.array([np.argmin(abs(self.freqs - f)) for f in f0]) + h0 = self.h_thresh(SNR=SNR_threshold) + # detection volume is is the sum of detection radius * pixel area over all pixels + volume = [dA/3.*np.sum( + strain_and_chirp_mass_to_luminosity_distance(h0[fdx], M_c, self.freqs[fdx])**3, + axis=0).value for fdx in f_idx] + return volume[0] if len(volume)==1 else volume + def h_circ(M_c, D_L, f0, Tspan, f): r""" diff --git a/hasasia/utils.py b/hasasia/utils.py index 8651e28..80c75d9 100644 --- a/hasasia/utils.py +++ b/hasasia/utils.py @@ -5,6 +5,9 @@ import scipy.special as ssp import scipy.integrate as si import scipy.stats as ss +from astropy.coordinates import SkyCoord +import astropy.units as u +import astropy.constants as c __all__ = ['create_design_matrix', 'fap', @@ -87,30 +90,115 @@ def pdf_F_signal(F, snr, Npsrs=None): N = int(4 * Npsrs) return ss.ncx2.pdf(2*F, N, snr**2) -def fdp(F0, snr, Npsrs=None, sky_ave=False): +def false_dismissal_prob(F0, snr, Npsrs=None, iota_psi_ave=False): ''' - False detection probability of the F-statistic + False dismissal probability of the F-statistic Use None for the Fe statistic and the number of pulsars for the Fp stat. ''' if Npsrs is None: N = 4 elif isinstance(Npsrs,int): N = int(4 * Npsrs) - if sky_ave: + if iota_psi_ave: return ss.chi2.cdf(2*F0, df=N, loc=snr**2) else: return ss.ncx2.cdf(2*F0, df=N, nc=snr**2) +def detection_prob(F0, snr, Npsrs=None, iota_psi_ave=False): + ''' + Detection probability of the F-statistic + Use None for the Fe and the number of pulsars for the Fp stat. + ''' + return 1 - false_dismissal_prob(F0, snr, Npsrs, iota_psi_ave) def _solve_F_given_fap(fap0=0.003, Npsrs=None): return sopt.fsolve(lambda F :fap(F, Npsrs=Npsrs)-fap0, 10) -def _solve_F_given_fdp_snr(fdp0=0.05, snr=3, Npsrs=None, sky_ave=False): +def _solve_F_given_fdp_snr(fdp0=0.05, snr=3, Npsrs=None, iota_psi_ave=False): Npsrs = 1 if Npsrs is None else Npsrs F0 = (4*Npsrs+snr**2)/2 - return sopt.fsolve(lambda F :fdp(F, snr, Npsrs=Npsrs, sky_ave=sky_ave)-fdp0, F0) + return sopt.fsolve(lambda F :false_dismissal_prob(F, snr, Npsrs=Npsrs, iota_psi_ave=iota_psi_ave)-fdp0, F0) -def _solve_snr_given_fdp_F(fdp0=0.05, F=3, Npsrs=None, sky_ave=False): +def _solve_snr_given_fdp_F(fdp0=0.05, F=3, Npsrs=None, iota_psi_ave=False): Npsrs = 1 if Npsrs is None else Npsrs snr0 = np.sqrt(2*F-4*Npsrs) - return sopt.fsolve(lambda snr :fdp(F, snr, Npsrs=Npsrs, sky_ave=sky_ave)-fdp0, snr0) + return sopt.fsolve(lambda snr :false_dismissal_prob(F, snr, Npsrs=Npsrs, iota_psi_ave=iota_psi_ave)-fdp0, snr0) + +def _solve_F0_given_SNR(snr=3, Npsrs=None): + ''' + Returns the F0 (Fe stat threshold for a specified SNR) + Use None for the Fe and the number of pulsars for the Fp stat. + ''' + Npsrs = 1 if Npsrs is None else Npsrs + return 0.5*(4.*Npsrs+snr**2.) + +def strain_and_chirp_mass_to_luminosity_distance(h, M_c, f0): + r''' + Returns the luminosity distance to a source given the strain, chirp mass, and GW frequency. + + Parameters + ---------- + h : float + The strain of the source. + M_c : float + The chirp mass of the source [Msun]. + f0 : float + The GW frequency of the source [Hz]. + + Returns + ------- + D_L : float + The luminosity distance to the source [Mpc]. + ''' + return (4*c.c / (h * u.m/u.m) + * np.power(c.G * M_c * u.Msun/c.c**3, 5/3) + * np.power(np.pi * f0 * u.Hz, 2/3)).to('Mpc') + +def theta_phi_to_SkyCoord(theta, phi): + r''' + Takes in a celestial longitude and lattitude and returns an `astropy.SkyCoord` object. + + Parameters + ---------- + phi : float, array of floats + The celestial longitude in solar system coordinates. + theta : float, array of floats + The celestial lattitude in solar system coordinates. + + Returns + ------- + skycoord - astropy.SkyCoord object + Can use this to convert to ra, dec, etc. + (e.g. SkyCoord.ra.deg) + + ''' + + return SkyCoord(phi*u.rad, ( theta - np.pi/2 )*u.rad) + +def skycoord_to_Jname(skycoord): + ''' + Takes in a SkyCoord object and returns the Jname of a pulsar with given coordinates. + + Parameters + ---------- + skycoord - astropy.SkyCoord object + Can use `theta_phi_to_SkyCoord()` to get this. + + Returns + ------- + Jname - string, array of strings + The traditional Jname of a pulsar with given coordinates + (eg. 'J1713+0747') + + ''' + coord_pieces = [ + str(skycoord.ra.hms[0]).split('.')[0], + str(skycoord.ra.hms[1]).split('.')[0], + str(abs(skycoord.dec.hms[0])).split('.')[0], + str(abs(skycoord.dec.hms[1])).split('.')[0] + ] + sign = ['-','+'][int(skycoord.dec.hms[0]>0)] + for i, piece in enumerate(coord_pieces): + if len(str(abs(int(piece)))) < 2: + coord_pieces[i] = '0' + str(piece) + return 'J' + coord_pieces[0] + coord_pieces[1] + sign + coord_pieces[2] +coord_pieces[3] \ No newline at end of file