diff --git a/pyfstat/__init__.py b/pyfstat/__init__.py index 3ada8dc9d78a31783b4807b0bd3933a9a83cd94b..d54f07893640c17d08af342373ccf90a74074328 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 0000000000000000000000000000000000000000..04d854fcfd96f72067f30cb8db4ce48041de5416 --- /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