add some utils and TDP function

hasasia/sensitivity.py

@@ -758,7 +758,11 @@ def to_pickle(self, filepath):
             pickle.dump(self, fout)
 
     def fidx(self,f):
-        return np.argmin(np.abs(self.freqs-f))
+        """Get the indices of a frequencies in the frequency array."""
+        if isinstance(f, int) or isinstance(f, float):
+            f = np.array([f])
+            f = np.asarray(f)
+        return np.array([np.argmin(abs(ff-self.freqs)) for ff in f])
 
     @property
     def S_eff(self):

hasasia/skymap.py

@@ -363,6 +363,24 @@ def detection_prob(self, F_thresh, snr=None, Npsrs=None, ave=None, prob_kwargs={
         Use None for the Fe and the number of pulsars for the Fp stat.
         '''
         return 1 - self.false_dismissal_prob(F_thresh, snr, Npsrs, ave, prob_kwargs)
+
+    def total_detection_probability(self, F_thresh, snr=None, Npsrs=None, ave=None, prob_kwargs={'h0':None, 'fidx': None}):
+        '''
+        See equation ## in Rosado.
+        Can be interpretted as the probability of detecting a single source
+        in *any* frequency bin.
+        
+        '''
+        return 1. - np.prod(
+            [self.false_dismissal_prob(
+                F_thresh,
+                snr=snr,
+                Npsrs=Npsrs,
+                ave=ave,
+                prob_kwargs={'h0':h0,'fidx':fidx})
+             for h0, fidx in zip(prob_kwargs['h0'],prob_kwargs['fidx'])
+            ],
+            axis=0)
 
     def _fdp_angle_averaged(self, F_thresh, h0, fidx):
         '''

hasasia/utils.py

@@ -93,6 +93,28 @@ def strain_and_chirp_mass_to_luminosity_distance(h, M_c, f0):
             * 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 char_strain_to_strain_amp(hc, fc, df):
+    r'''
+    Calculate the strain amplitude of single sources given
+    their characteristic strains.
+
+    Parameters
+    ----------
+    hc : array_like
+        Characteristic strain of the single sources.
+    fc : array_like
+        Observed orbital frequency bin centers.
+    df : array_like
+        Observed orbital frequency bin widths.
+
+    Returns
+    -------
+    hs : 
+        Strain amplitude of the single sources.
+
+    '''
+    return hc * np.sqrt(df/fc)
+
 def theta_phi_to_SkyCoord(theta, phi):
     r'''
     Takes in a celestial longitude and lattitude and returns an `astropy.SkyCoord` object.
@@ -140,4 +162,29 @@ def skycoord_to_Jname(skycoord):
     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]
+    return 'J' + coord_pieces[0] + coord_pieces[1] + sign + coord_pieces[2] +coord_pieces[3]
+
+def distance_on_sphere(lat1, lon1, lat2, lon2):
+    ''''
+    Returns the distance between two points on the surface of the unit sphere.
+    
+    Parameters
+    ==========
+    lat1, lat2 - float
+        lattitude of first and second point respectively
+    lon1, lon2 - float
+        longitude of first and second point respectively
+    Returns
+    =======
+    distance - float
+        distance between two points
+    '''
+    dlat = lat2 - lat1
+    dlon = lon2 - lon1
+
+    # Haversine formula for distance
+    a = np.sin(dlat / 2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2)**2
+    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))
+    distance = c
+
+    return distance