From 0ca51fa73b51b12adf45bcc5e7b94db0dcff5003 Mon Sep 17 00:00:00 2001
From: "gregory.ashton" <gregory.ashton@ligo.org>
Date: Thu, 4 May 2017 17:08:55 +0200
Subject: [PATCH] Adds method to calculate frequency range of signal
---
pyfstat/__init__.py | 1 +
pyfstat/injection_helper_functions.py | 82 +++++++++++++++++++++++++++
2 files changed, 83 insertions(+)
create mode 100644 pyfstat/injection_helper_functions.py
diff --git a/pyfstat/__init__.py b/pyfstat/__init__.py
index 3ada8dc..d54f078 100644
--- a/pyfstat/__init__.py
+++ b/pyfstat/__init__.py
@@ -3,4 +3,5 @@ from __future__ import division as _division
from .core import BaseSearchClass, ComputeFstat, Writer, SemiCoherentSearch, SemiCoherentGlitchSearch
from .mcmc_based_searches import MCMCSearch, MCMCGlitchSearch, MCMCSemiCoherentSearch, MCMCFollowUpSearch, MCMCTransientSearch
from .grid_based_searches import GridSearch, GridUniformPriorSearch, GridGlitchSearch, FrequencySlidingWindow
+from .injection_helper_functions import get_frequency_range_of_signal
diff --git a/pyfstat/injection_helper_functions.py b/pyfstat/injection_helper_functions.py
new file mode 100644
index 0000000..04d854f
--- /dev/null
+++ b/pyfstat/injection_helper_functions.py
@@ -0,0 +1,82 @@
+""" Code used with permissision from Sylvia Zhu to calculate the range in
+ frequency space that a signal occupies due to spindown and Doppler
+ modulations
+"""
+
+import numpy as np
+from astropy import units as u
+from astropy.coordinates import SkyCoord
+from astropy.time import Time
+
+# 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
+# equinox, or March 20.
+# Each day the Earth moves 2*pi/365 radians around its orbit.
+
+
+def _eqToEcl(alpha, delta):
+ source = SkyCoord(alpha*u.radian, delta*u.radian, frame='gcrs')
+ out = source.transform_to('geocentrictrueecliptic')
+ return np.array([out.lon.radian, out.lat.radian])
+
+
+def _eclToEq(lon, lat):
+ source = SkyCoord(lon*u.radian, lat*u.radian,
+ frame='geocentrictrueecliptic')
+ out = source.transform_to('gcrs')
+ return np.array([out.ra.radian, out.dec.radian])
+
+
+def _calcDopplerWings(
+ s_freq, s_alpha, s_delta, lonStart, lonStop, numTimes=100):
+ e_longitudes = np.linspace(lonStart, lonStop, numTimes)
+ v_over_c = 1e-4
+ s_lon, s_lat = _eqToEcl(s_alpha, s_delta)
+
+ vertical = s_lat
+ horizontals = s_lon - e_longitudes
+
+ dopplerShifts = s_freq * np.sin(horizontals) * np.cos(vertical) * v_over_c
+ return np.amin(dopplerShifts), np.amax(dopplerShifts)
+
+
+def _calcSpindownWings(freq, fdot, minStartTime, maxStartTime):
+ timespan = maxStartTime - minStartTime
+ return 0.5 * timespan * np.abs(fdot) * np.array([-1, 1])
+
+
+def get_frequency_range_of_signal(F0, F1, Alpha, Delta, minStartTime,
+ maxStartTime):
+ """ Calculate the frequency range that a signal will occupy
+
+ Parameters
+ ----------
+ F0, F1, Alpha, Delta: float
+ Frequency, derivative, and sky position for the signal (all angles in
+ radians)
+ minStartTime, maxStartTime: float
+ GPS time of the start and end of the data span
+
+ Returns
+ -------
+ [Fmin, Fmax]: array
+ The minimum and maximum frequency span
+ """
+ YEAR_IN_DAYS = 365.25
+ tEquinox = 79
+
+ minStartTime_t = Time(minStartTime, format='gps').to_datetime().timetuple()
+ maxStartTime_t = Time(minStartTime, 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
+
+ dopplerWings = _calcDopplerWings(F0, Alpha, Delta, lonStart, lonStop)
+ spindownWings = _calcSpindownWings(F0, F1, minStartTime, maxStartTime)
+ return np.array([F0, F0]) + dopplerWings + spindownWings
--
GitLab