Adds compute pstar functionality

pyfstat/grid_based_searches.py

@@ -501,7 +501,7 @@ class DMoff_NO_SPIN(GridSearch):
         self.c = 2.998e8
         self.SIDEREAL_DAY = 23*60*60 + 56*60 + 4.0916
         self.TERRESTRIAL_DAY = 86400.
-        a0 = self.Re/self.c*np.cos(self.par['Delta'])
+        a0 = self.Re/self.c  # *np.cos(self.par['Delta'])
         self.m0 = np.max([4, int(np.ceil(2*np.pi*self.par['F0']*a0))])
         logging.info('m0 = {}'.format(self.m0))
 
@@ -516,7 +516,8 @@ class DMoff_NO_SPIN(GridSearch):
         self.SSBprec = 2
         self.out_file = '{}/{}_gridFS_SSBPREC2.txt'.format(
             self.outdir, self.label)
-        self.F0s = [self.par['F0']+j/self.SIDEREAL_DAY for j in range(-4, 5)]
+        self.F0s = [self.par['F0']+j/self.SIDEREAL_DAY
+                    for j in range(-self.m0, self.m0+1)]
         self.run()
         twoF_SUM = np.sum(self.data[:, -1])
 

pyfstat/helper_functions.py

@@ -9,6 +9,7 @@ import logging
 import inspect
 import peakutils
 from functools import wraps
+from scipy.stats.distributions import ncx2
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -151,3 +152,34 @@ def get_comb_values(F0, frequencies, twoF, period, N=4):
     comb_idxs = [np.argmin(np.abs(frequencies-F0-F)) for F in comb_frequencies]
     return comb_frequencies, twoF[comb_idxs], freq_err*np.ones(len(comb_idxs))
 
+
+def compute_P_twoFstarcheck(twoFstarcheck, twoFcheck, M0, plot=False):
+    """ Returns the unnormalised pdf of twoFstarcheck given twoFcheck """
+    upper = 4+twoFstarcheck + 0.5*(2*(4*M0+2*twoFcheck))
+    rho2starcheck = np.linspace(1e-1, upper, 500)
+    integrand = (ncx2.pdf(twoFstarcheck, 4*M0, rho2starcheck)
+                 * ncx2.pdf(twoFcheck, 4, rho2starcheck))
+    if plot:
+        fig, ax = plt.subplots()
+        ax.plot(rho2starcheck, integrand)
+        fig.savefig('test')
+    return np.trapz(integrand, rho2starcheck)
+
+
+def compute_pstar(twoFcheck_obs, twoFstarcheck_obs, m0, plot=False):
+    M0 = 2*m0 + 1
+    upper = 4+twoFcheck_obs + (2*(4*M0+2*twoFcheck_obs))
+    twoFstarcheck_vals = np.linspace(1e-1, upper, 500)
+    P_twoFstarcheck = np.array(
+        [compute_P_twoFstarcheck(twoFstarcheck, twoFcheck_obs, M0)
+         for twoFstarcheck in twoFstarcheck_vals])
+    C = np.trapz(P_twoFstarcheck, twoFstarcheck_vals)
+    idx = np.argmin(np.abs(twoFstarcheck_vals - twoFstarcheck_obs))
+    if plot:
+        fig, ax = plt.subplots()
+        ax.plot(twoFstarcheck_vals, P_twoFstarcheck)
+        ax.fill_between(twoFstarcheck_vals[:idx+1], 0, P_twoFstarcheck[:idx+1])
+        ax.axvline(twoFstarcheck_vals[idx])
+        fig.savefig('test')
+    pstar_l = np.trapz(P_twoFstarcheck[:idx+1]/C, twoFstarcheck_vals[:idx+1])
+    return 2*np.min([pstar_l, 1-pstar_l])