From 25b45c0b9ed7da937511ebc563ebae6184d0272f Mon Sep 17 00:00:00 2001
From: Gregory Ashton <gregory.ashton@ligo.org>
Date: Wed, 2 Aug 2017 09:57:05 +0200
Subject: [PATCH] Update of functions to calculate the covering band
- Adds function to use XLALCWSignalCoveringBand
- Fixes issues in injection_helper_function, both are checked to agree
---
pyfstat/helper_functions.py | 23 +++++++++++++++++++++++
pyfstat/injection_helper_functions.py | 12 ++++++------
2 files changed, 29 insertions(+), 6 deletions(-)
diff --git a/pyfstat/helper_functions.py b/pyfstat/helper_functions.py
index 7e2226b..d6eece2 100644
--- a/pyfstat/helper_functions.py
+++ b/pyfstat/helper_functions.py
@@ -239,3 +239,26 @@ def get_sft_array(sftfilepattern, data_duration, F0, dF0):
times = np.linspace(0, data_duration, nsfts)
return times, freqs, data
+
+
+def get_covering_band(tref, tstart, tend, uniform_prior):
+ tref = lal.LIGOTimeGPS(tref)
+ tstart = lal.LIGOTimeGPS(tstart)
+ tend = lal.LIGOTimeGPS(tend)
+ psr = lalpulsar.PulsarSpinRange()
+ for i, key in enumerate(['F0', 'F1', 'F2']):
+ if key in uniform_prior:
+ if type(uniform_prior[key]) == dict and uniform_prior[key]['type'] == 'unif':
+ l, u = uniform_prior[key]['lower'], uniform_prior[key]['upper']
+ psr.fkdot[i] = (l+u)/2.
+ psr.fkdotBand[i] = u-l
+ else:
+ psr.fkdot[i] = uniform_prior[key]
+ psr.fkdotBand[i] = 0
+ else:
+ raise ValueError(
+ 'uniform_prior should contain unif or const values of F0, F1, F2')
+ psr.refTime = tref
+ return lalpulsar.CWSignalCoveringBand(tstart, tend, psr, 0, 0, 0)
+
+
diff --git a/pyfstat/injection_helper_functions.py b/pyfstat/injection_helper_functions.py
index 972115c..3698518 100644
--- a/pyfstat/injection_helper_functions.py
+++ b/pyfstat/injection_helper_functions.py
@@ -11,6 +11,7 @@ try:
from astropy.time import Time
except ImportError:
logging.warning('Python module astropy not installed')
+import lal
# Assume Earth goes around Sun in a non-wobbling circle at constant speed;
# Still take the zero longitude to be the Earth's position during the March
@@ -34,7 +35,7 @@ def _eclToEq(lon, lat):
def _calcDopplerWings(
s_freq, s_alpha, s_delta, lonStart, lonStop, numTimes=100):
e_longitudes = np.linspace(lonStart, lonStop, numTimes)
- v_over_c = 1e-4
+ v_over_c = 2*np.pi*lal.AU_SI/lal.YRSID_SI/lal.C_SI
s_lon, s_lat = _eqToEcl(s_alpha, s_delta)
vertical = s_lat
@@ -61,23 +62,22 @@ def get_frequency_range_of_signal(F0, F1, Alpha, Delta, minStartTime,
minStartTime, maxStartTime: float
GPS time of the start and end of the data span
+ Note: assumes tref is in the middle of the data span
+
Returns
-------
[Fmin, Fmax]: array
The minimum and maximum frequency span
"""
- YEAR_IN_DAYS = 365.25
+ YEAR_IN_DAYS = lal.YRSID_SI / lal.DAYSID_SI
tEquinox = 79
minStartTime_t = Time(minStartTime, format='gps').to_datetime().timetuple()
- maxStartTime_t = Time(minStartTime, format='gps').to_datetime().timetuple()
+ maxStartTime_t = Time(maxStartTime, format='gps').to_datetime().timetuple()
tStart_days = minStartTime_t.tm_yday - tEquinox
tStop_days = maxStartTime_t.tm_yday - tEquinox
tStop_days += (maxStartTime_t.tm_year-minStartTime_t.tm_year)*YEAR_IN_DAYS
- tStart_days = 280 - tEquinox # 7 October is day 280 in a non leap year
- tStop_days = 19 + YEAR_IN_DAYS - tEquinox # the next year
-
lonStart = 2*np.pi*tStart_days/YEAR_IN_DAYS - np.pi
lonStop = 2*np.pi*tStop_days/YEAR_IN_DAYS - np.pi
--
GitLab