waveform_util

Header file for waveform specific utilites

Functions

double match(std::complex<double> *data1, std::complex<double> *data2, double *SN, double *frequencies, int length)

template<class T>
void create_coherent_GW_detection(std::string *detectors, int detector_N, T **frequencies, int *lengths, bool reuse_WF, gen_params_base<T> *gen_params, std::string generation_method, std::complex<T> **responses)

Parameters:

reuse_WF – If using the same exact frequencies for each detector, we can save a lot of time by only computing the waveform once
gen_params – tc in gen_params refers to the time of coalescence, relative to the initial time of the data stream, for the FIRST detector, the other detectors are shifted by the appropriate amount
responses – [out] Responses for the source at each detector, same order as detectors parameter — should be pre allocated shape [detector_N][lengths[i]]

template<class T> void create_single_GW_detection(std::complex<T> *response, std::string detector, T *frequencies, int length, gen_params_base<T> *gen_params, std::string generation_method)

template<class T>
void create_coherent_GW_detection_reuse_WF(std::string *detectors, int detector_N, T *frequencies, int lengths, gen_params_base<T> *gen_params, std::string generation_method, std::complex<T> **responses)

Parameters:

gen_params – tc in gen_params refers to the time of coalescence, relative to the initial time of the data stream, for the FIRST detector, the other detectors are shifted by the appropriate amount
responses – [out] Responses for the source at each detector, same order as detectors parameter — should be pre allocated shape [detector_N][length]

double data_snr(double *frequencies, int length, std::complex<double> *data, std::complex<double> *response, double *psd)

double data_snr_maximized_extrinsic(double *frequencies, int length, std::complex<double> *data, double *psd, std::string detector, std::string generation_method, gen_params *param)

Utility to calculate the snr of a fourier transformed data stream while maximizing over the coalescence parameters phic and tc.

The gen_params structure holds the parameters for the template to be used (the maximimum likelihood parameters)

Parameters:

frequencies – Frequencies used by data
length – length of the data
data – input data in the fourier domain
psd – PSD for the detector that created the data
detector – Name of the detector —See noise_util for options
generation_method – Generation method for the template — See waveform_generation.cpp for options
param – gen_params structure for the template

double data_snr_maximized_extrinsic(double *frequencies, int length, double *data_real, double *data_imag, double *psd, std::string detector, std::string generation_method, gen_params *param)

Light wrapper for the data_snr_maximized_extrinsic method.

Splits the data into real and imaginary, so all the arguments are C-safe

Parameters:

frequencies – Frequencies used by data
length – length of the data
data_real – input data in the fourier domain — real part
data_imag – input data in the fourier domain — imaginary part
psd – PSD for the detector that created the data
detector – Name of the detector —See noise_util for options
generation_method – Generation method for the template — See waveform_generation.cpp for options
param – gen_params structure for the template

double calculate_snr(std::string sensitivity_curve, std::complex<double> *waveform, double *frequencies, int length, std::string integration_method = "SIMPSONS", double *weights = NULL, bool log10_freq = false)

Caclulates the snr given a detector and waveform (complex) and frequencies.

This function computes the un-normalized snr: \sqrt( ( H | H ) )

Parameters:

sensitivity_curve – detector name - must match the string of populate_noise precisely
waveform – complex waveform
frequencies – double array of frequencies that the waveform is evaluated at
length – length of the above two arrays

double calculate_snr_internal(double *psd, std::complex<double> *waveform, double *frequencies, int length, std::string integration_method = "SIMPSONS", double *weights = NULL, bool log10_freq = false)

double calculate_snr_internal(double *psd, std::complex<double> *waveform, const Quadrature *QuadMethod)

double calculate_snr(std::string sensitivity_curve, std::string detector, std::string generation_method, gen_params_base<double> *params, double *frequencies, int length, std::string integration_method = "SIMPSONS", double *weights = NULL, bool log10_freq = false)

Routine to calculate the SNR of a template with GSL quadrature integration.

Sometimes, this is faster than the `grid’ style integration

Supports sky-averaged templates, but this should only be used with non-precessing waveforms

int calculate_snr_gsl(double *snr, std::string sensitivity_curve, std::string detector, std::string generation_method, gen_params_base<double> *params, double f_min, double f_max, double relative_error)

Routine to calculate the SNR of a template with GSL quadrature integration.

Sometimes, this is faster than the `grid’ style integration

Supports sky-averaged templates, but this should only be used with non-precessing waveforms

Parameters:

snr – [out] SNR
sensitivity_curve – Noise curve
detector – Detector to compute response — can be empty is SA
generation_method – Generation method
params – Source Parameters
f_min – Lower frequency bound
f_max – Upper frequency bound
relative_error – Relative error threshold

int calculate_snr_gsl(double *snr, std::string sensitivity_curve, std::string detector, std::string generation_method, gen_params_base<double> *params, double f_min, double f_max, double relative_error, gsl_integration_workspace *w, int np)

Parameters:

sensitivity_curve – Noise curve
detector – Detector to compute response — can be empty is SA
generation_method – Generation method
params – Source Parameters
f_min – Lower frequency bound
f_max – Upper frequency bound
relative_error – Relative error threshold
w – User-allocated gsl_integration_workspace
np – Size of gsl_integration_workspace allocation

double integrand_snr_SA_subroutine(double f, void *subroutine_params)

Internal function to calculate the SNR integrand for sky-averaged waveforms.

NOTE: Only works for GR spin aligned

double integrand_snr_subroutine(double f, void *subroutine_params): Internal function to calculate the SNR integrand for full waveforms.

template<class T> int fourier_detector_response_horizon(T *frequencies, int length, waveform_polarizations<T>*, std::complex<T> *detector_response, T theta, T phi, std::string detector)

Parameters:

frequencies – array of frequencies corresponding to waveform
length – length of frequency/waveform arrays
detector_response – [out] detector response
theta – polar angle (rad) theta in detector frame
phi – azimuthal angle (rad) phi in detector frame
detector – detector - list of supported detectors in noise_util

template<class T> int fourier_detector_response_horizon(T *frequencies, int length, waveform_polarizations<T>*, std::complex<T> *detector_response, T theta, T phi, T psi, std::string detector)

Parameters:

frequencies – array of frequencies corresponding to waveform
length – length of frequency/waveform arrays
detector_response – [out] detector response
theta – polar angle (rad) theta in detector frame
phi – azimuthal angle (rad) phi in detector frame
psi – polarization angle (rad) phi in detector frame
detector – detector - list of supported detectors in noise_util

template<class T> int time_detector_response_horizon(T *times, int length, std::complex<T> *response, std::string detector, std::string generation_method, gen_params_base<T> *parameters)

Parameters:

times – double array of frequencies for the waveform to be evaluated at
length – integer length of all the arrays
response – [out] complex array for the output plus polarization waveform
generation_method – String that corresponds to the generation method - MUST BE SPELLED EXACTLY
parameters – structure containing all the source parameters

template<class T> int fourier_detector_response_horizon(T *frequencies, int length, std::complex<T> *response, std::string detector, std::string generation_method, gen_params_base<T> *parameters)

Function to produce the detector response caused by impinging gravitational waves from a quasi-circular binary.

By using the structure parameter, the function is allowed to be more flexible in using different method of waveform generation - not all methods use the same parameters

This puts the responsibility on the user to pass the necessary parameters

Detector options include classic interferometers like LIGO/VIRGO (coming soon: ET and LISA)

This is a wrapper that combines generation with response functions: if producing mulitple responses for one waveform (ie stacking Hanford, Livingston, and VIRGO), it will be considerably more efficient to calculate the waveform once, then combine each response manually

Parameters:

frequencies – double array of frequencies for the waveform to be evaluated at
length – integer length of all the arrays
response – [out] complex array for the output plus polarization waveform
generation_method – String that corresponds to the generation method - MUST BE SPELLED EXACTLY
parameters – structure containing all the source parameters

template<class T> int fourier_detector_response_equatorial(T *frequencies, int length, waveform_polarizations<T>*, std::complex<T> *detector_response, T ra, T dec, T psi, double gmst, T *times, T LISA_alpha0, T LISA_phi0, T theta_j_ecl, T phi_j_ecl, std::string detector)

Parameters:

frequencies – array of frequencies corresponding to waveform
length – length of frequency/waveform arrays
detector_response – [out] detector response
ra – Right Ascension in rad
dec – Declination in rad
psi – polarization angle (rad)
gmst – greenwich mean sidereal time
detector – detector - list of supported detectors in noise_util

template<class T> int time_detector_response_equatorial(T *times, int length, std::complex<T> *response, std::string detector, std::string generation_method, gen_params_base<T> *parameters)

Parameters:

times – double array of frequencies for the waveform to be evaluated at
length – integer length of all the arrays
response – [out] complex array for the output plus polarization waveform
generation_method – String that corresponds to the generation method - MUST BE SPELLED EXACTLY
parameters – structure containing all the source parameters

template<class T> int fourier_detector_response_equatorial(T *frequencies, int length, std::complex<T> *response, std::string detector, std::string generation_method, gen_params_base<T> *parameters)

Function to produce the detector response caused by impinging gravitational waves from a quasi-circular binary for equatorial coordinates for TERRESTIAL detectors, where the earth’s rotation during the signal is minor.

By using the structure parameter, the function is allowed to be more flexible in using different method of waveform generation - not all methods use the same parameters

This puts the responsibility on the user to pass the necessary parameters

Detector options include classic interferometers like LIGO/VIRGO (coming soon: ET and LISA)

This is a wrapper that combines generation with response functions: if producing mulitple responses for one waveform (ie stacking Hanford, Livingston, and VIRGO), it will be considerably more efficient to calculate the waveform once, then combine each response manually

Parameters:

frequencies – double array of frequencies for the waveform to be evaluated at
length – integer length of all the arrays
response – [out] complex array for the output plus polarization waveform
generation_method – String that corresponds to the generation method - MUST BE SPELLED EXACTLY
parameters – structure containing all the source parameters

template<class T> int fourier_detector_response_equatorial(T *frequencies, int length, std::complex<T> *response, std::string detector, std::string generation_method, gen_params_base<T> *parameters, T *times)

Function to produce the detector response caused by impinging gravitational waves from a quasi-circular binary for equatorial coordinates.

By using the structure parameter, the function is allowed to be more flexible in using different method of waveform generation - not all methods use the same parameters

This puts the responsibility on the user to pass the necessary parameters

Detector options include classic interferometers like LIGO/VIRGO (coming soon: ET and LISA)

This is a wrapper that combines generation with response functions: if producing mulitple responses for one waveform (ie stacking Hanford, Livingston, and VIRGO), it will be considerably more efficient to calculate the waveform once, then combine each response manually

Parameters:

frequencies – double array of frequencies for the waveform to be evaluated at
length – integer length of all the arrays
response – [out] complex array for the output plus polarization waveform
generation_method – String that corresponds to the generation method - MUST BE SPELLED EXACTLY
parameters – structure containing all the source parameters

template<class T> int time_detector_response(T *times, int length, std::complex<T> *response, std::string detector, std::string generation_method, gen_params_base<T> *parameters)

template<class T> int fourier_detector_response(T *frequencies, int length, std::complex<T> *response, std::string detector, std::string generation_method, gen_params_base<T> *parameters, T *times = NULL): Wrapper to handle all detector_response calls — horizon and equatorial.

int boundary_number(std::string method)

Utility to inform the fisher routine how many logical boundaries should be expected.

The automatic derivative code requires a new tape for each logical branch of the program, so each waveform_generation method needs to add the number of branches here

void time_phase_corrected_autodiff(double *times, int length, double *frequencies, gen_params_base<double> *params, std::string generation_method, bool correct_time, int *tapes_in = NULL, int order = 1)

Mapping from phase to time AUTODIFFERENTIATION using ADOL-C.

Currently, only IMRPhenomD is used for the time-phase relationship, as this is simpler and captures almost all of the physics. The generation_method parameter does NOTHING right now.

This is NOT autodiff safe. ADOL-C does not support wrapping a section of code as active twice. To find the derivative of this function, you must use the hessian function of ADOL-C

Made for use with detectors like LISA

Uses the deprecated postmerger calculation for IMRPhenomD — this is ADOL-C friendly

Using https://arxiv.org/abs/1809.04799 t = (1/2PI) d phi/ d f

correct_time is currently not supported

For IMRPhenomPv2, the phase has to be wrapped, because arctan is taken of the waveform because of the euler rotations. This might make the numerical derivative unpredictable

if the order is set to 2, returns (1/2PI) d^2 phi / d f^2

template<class T> void time_phase_corrected(T *times, int length, T *frequencies, gen_params_base<T> *params, std::string generation_method, bool correct_time, int order = 1)

Mapping from phase to time NUMERICAL.

Currently, only IMRPhenomD is used for the time-phase relationship, as this is simpler and captures almost all of the physics. The generation_method parameter does NOTHING right now.

Made for use with detectors like LISA

Using https://arxiv.org/abs/1809.04799 t = (1/2PI) d phi/ d f

correct_time is currently not supported

For IMRPhenomPv2, the phase has to be wrapped, because arctan is taken of the waveform because of the euler rotations. This might make the numerical derivative unpredictable

Just uses a second order numerical derivative for now

Second order is only available for 1 frequency at a time — something’s funky with this

int fourier_detector_amplitude_phase(double *frequencies, int length, double *amplitude, double *phase, std::string detector, std::string generation_method, gen_params *parameters)

Calculates the amplitude (magnitude) and phase (argument) of the response of a given detector.

This is for general waveforms, and will work for precessing waveforms

Not as fast as non-precessing, but that can’t be helped. MUST include plus/cross polarizations

template<class T> void transform_orientation_coords(gen_params_base<T> *parameters, std::string generation_method, std::string detector)

template<class T> void map_extrinsic_angles(gen_params_base<T> *params)

void assign_freq_boundaries(double *freq_boundaries, double *intermediate_freqs, int boundary_num, gen_params_base<double> *input_params, std::string generation_method)

For ADOL-C, assigns the frequencies at which logical breaks occur in the waveform.

ADOL-C cannot handle logically branches, so a tape must be made for each branch separately. This function returns the locations of these branches

void integration_bounds(gen_params_base<double> *params, std::string generation_method, std::string detector, std::string sensitivity_curve, double fmin, double fmax, double signal_to_noise, double tol, double *integration_bounds, bool autodiff)

Utility to find the integration bounds for Fisher matrices for increasing speed of Fisher evaluation.

Numerically finds the frequencies at which the Fisher should be evaluated at.

Uses the bisection search algorithm for the cases where the waveform enters/leaves the band at SNR>1

Uses a 100 pt grid search (logarithmically spaced) if the signal has SNR<1 when entering and leaving

integrand_bounds[0] ~ frequency at which |h|/(sqrt S) ~signal_to_noise +/- tol

integrand_bounds[1] ~ frequency at which |h|/(sqrt S) ~signal_to_noise +/- tol

If autodiff set to true, ADOL-C routines are used. This is openmp safe, not necessarily thread safe.

Parameters:

params – Parameters of the waveform
detector – Detector to use for the response function
sensitivity_curve – Sensitivity curve to use (must be one of the analytic curves in the detector_utilitiy file
fmin – minimum frequency to use (specific to the detector)
fmax – max frequency to use (specific to the detector)
signal_to_noise – Target ratio of |h|/ sqrt(S) (typically ~0.1)
tol – This is a numerical algorithm, so the tolerance must be specified
integration_bounds – [out] bounds fo the integral shape — [2] — (fmin,fmax)

int observation_bounds(double sampling_freq, double integration_time, std::string detector, std::string sensitivity_curve, std::string generation_method, gen_params_base<double> *params, double *freq_bounds, bool autodiff)

Determines the integration bounds for the log likelihood or fisher given some observation time, sampling frequency, detector, and sensitivity curve.

Sensitivity curve has to be one of the options in detector_util analytic options

The current scheme is to use the frequency bounds determined by the SNR if the binary spends less than the integration time in band. If the merger spends more time in band than the integration time, the frequencies are determined to be (f_integration_time, f_high_band)

Accounts for the integration time, unlike the integration_bounds routine

returns 0 if successful, returns 1 if bounds could not be found due to roundoff, and returns 2 if entirely unsuccessful

autodiff set to true will use ADOL-C for time-phase relationships. This is more accurate, but slower and only threadsafe with OpenMP

Parameters:

sampling_freq – Frequency at which the detector operates
integration_time – Time of observation in seconds
detector – Detector to use for the response function
sensitivity_curve – Sensitivity curve to use — must match analytic choices in detector_util
generation_method – method to use for the waveform generation
params – parameters of the source
freq_bounds – [out] Output bounds

int Tbm_to_freq(gen_params_base<double> *params, std::string generation_method, double Tbm, double *freq, double tol, bool autodiff, int max_iteration, bool relative_time)

Convenience function to Calculate the time before merger using numerical methods.

Tbm should be positive

tol refers to the tolerance of the bisection search. If using numerical derivatives, tolerance is not guaranteed because of errors introduced in the phase-time conversion

Also, tolerance is associated with Tbm. Closer to merger, this can correspond to much larger errors on frequency because the relationship is so steep between time and frequency near merger.

Uses numerical if autodiff set to false— omp safe and thread safe

Numerical accuracy is spotty. The second derivative isn’t always reliable

Relative time:

true -- Tbm is interpreted as the time before tc, and is shifted accordingly

false -- Tbm is interpreted as some absolute time, and is not shifted. Assumed to be in the same time coordinates as tc. ie, if Tbm set to 0, freq will be the frequency at a time tc before tc

Return values — See util.cpp, newton_raphson_method_1d

Parameters:

params – Generation parameters of the source
generation_method – Generation method for the waveform
Tbm – target time before merger — in seconds or years (years if YEARS==true)
freq – Frequency at the input time before merger
tol – Tolerance for the scheme
autodiff – Use autodiff routines instead of numerical derivatives
relative_time – True if shifting time by tc (time before merger) and false if not shifting time (absolute time frequency)

void Tbm_subroutine(double f, double *t, double *tp, void *param)

template<class T> void postmerger_params(gen_params_base<T> *params, std::string generation_method, T *fpeak, T *fdamp, T *fRD)

Calculates the postmerger parameters for a parameter set.

calculates fpeak, fdamp, and fRD

void threshold_times(gen_params_base<double> *params, std::string generation_method, double T_obs, double T_wait, double f_lower, double f_upper, std::string SN, double SNR_thresh, double *threshold_times_out, double tolerance)

Utility for calculating the threshold times before merger that result in an SNR>SNR_thresh

See arXiv 1902.00021

Binary must merge within time T_wait

SNR is calculated with frequencies [f(t_mer),f(t_mer-T_obs)] or [f(t_mer),0] depending on whether the binary has merged or not

Assumes sky average — Only supports PhenomD for now

If no time before merger satisfies the requirements, both are set to -1

Parameters:

generation_method – Generation method to use for the waveform
T_obs – Observation time — also specifies the frequency spacing (\delta f = 1./T_obs)
T_wait – Wait time — Maximum time for binaries to coalesce
f_lower – Lower bound of search
f_upper – upper bound of search
SN – Noise curve name
SNR_thresh – Threshold SNR
threshold_times_out – [out] Output frequencies
tolerance – Percent tolerance on SNR search

void threshold_times(gen_params_base<double> *params, std::string generation_method, double T_obs, double T_wait, double *freqs, double *SN, int length, double SNR_thresh, double *threshold_times_out, double tolerance)

Utility for calculating the threshold times before merger that result in an SNR>SNR_thresh.

See arXiv 1902.00021

Binary must merge within time T_wait

SNR is calculated with frequencies [f(t_mer),f(t_mer-T_obs)] or [f(t_mer),0] depending on whether the binary has merged or not

Assumes sky average — Only supports PhenomD for now — No angular dependence used ( only uses plus polarization — assumes iota = psi = 0 )

Assumes this is for multiband — ie stellar mass BHs — Only uses pn approximation of time frequency relation

If no time before merger satisfies the requirements, both are set to -1

Only supports tensor polarizations !!!

Parameters:

generation_method – Generation method to use for the waveform
T_obs – Observation time — also specifies the frequency spacing (\delta f = 1./T_obs)
T_wait – Wait time — Maximum time for binaries to coalesce
freqs – Maximum frequency array
SN – Noise curve array, should be prepopulated from f_lower to f_upper with spacing 1./T_obs
length – Length of maximum frequency array
SNR_thresh – Threshold SNR
threshold_times_out – [out] Output frequencies
tolerance – Percent tolerance on SNR search

double integrand_threshold_subroutine(double f, void *subroutine_params)

double snr_threshold_subroutine(double fmin, double fmax, double rel_err, gen_params_base<double> *params, std::string generation_method, std::string SN, gsl_integration_workspace *w, int np)

int threshold_times_gsl(gen_params_base<double> *params, std::string generation_method, double T_obs, double T_wait, double fmin, double fmax, std::string SN, double SNR_thresh, double *threshold_times_out, double *T_obs_SNR, double tolerance, gsl_integration_workspace *w, int np)

Utility for calculating the threshold times before merger that result in an SNR>SNR_thresh —GSL quad integration implementation.

See arXiv 1902.00021

Binary must merge within time T_wait

SNR is calculated with frequencies [f(t_mer),f(t_mer-T_obs)] or [f(t_mer),0] depending on whether the binary has merged or not

Assumes sky average — Only supports PhenomD for now — No angular dependence used ( only uses plus polarization — assumes iota = psi = 0 )

If no time before merger satisfies the requirements, both are set to -1

ALL temporal quantities in seconds or Hz

Return values:

-- success

11-- Failure: SNR was 0 in lower bound

12-- Failure: SNR was 0 in upper bound

-- partial success: Closest values output, but roundoff error prevented the routine from reaching the desired accuracy

-- partial success: numerical inversion for time -> frequency had errors 

-- partial success: numerical inversion for time -> frequency had failed, and a PN approximate had to be used

-- max iterations reached

Parameters:

generation_method – Generation method to use for the waveform
T_obs – Observation time — also specifies the frequency spacing (\delta f = 1./T_obs)
T_wait – Wait time — Maximum time for binaries to coalesce
fmin – Maximum frequency array
fmax – Maximum frequency array
SN – Noise curve array, should be prepopulated from f_lower to f_upper with spacing 1./T_obs
SNR_thresh – Threshold SNR
threshold_times_out – [out] Output times
tolerance – Percent tolerance on SNR search

int threshold_times_full_gsl(gen_params_base<double> *params, std::string generation_method, double T_obs, double T_wait, double fmin, double fmax, std::string SN, double SNR_thresh, double *threshold_times_out, double tolerance, gsl_integration_workspace *w, int np)

NOT FINISHED DO NOT USE Utility for calculating the threshold times before merger that result in an SNR>SNR_thresh —GSL quad integration implementation

NOT FINISHED — DO NOT USE

See arXiv 1902.00021

Binary must merge within time T_wait

SNR is calculated with frequencies [f(t_mer),f(t_mer-T_obs)] or [f(t_mer),0] depending on whether the binary has merged or not

Assumes sky average — Only supports PhenomD for now — No angular dependence used ( only uses plus polarization — assumes iota = psi = 0 )

If no time before merger satisfies the requirements, both are set to -1

ALL temporal quantities in seconds or Hz

Return values:

-- success

11-- Failure: SNR was 0 in lower bound

12-- Failure: SNR was 0 in upper bound

-- partial success: Closest values output, but roundoff error prevented the routine from reaching the desired accuracy

-- partial success: numerical inversion for time -> frequency had errors 

-- partial success: numerical inversion for time -> frequency had failed, and a PN approximate had to be used

-- max iterations reached

Parameters:

generation_method – Generation method to use for the waveform
T_obs – Observation time — also specifies the frequency spacing (\delta f = 1./T_obs)
T_wait – Wait time — Maximum time for binaries to coalesce
fmin – Maximum frequency array
fmax – Maximum frequency array
SN – Noise curve array, should be prepopulated from f_lower to f_upper with spacing 1./T_obs
SNR_thresh – Threshold SNR
threshold_times_out – [out] Output times
tolerance – Percent tolerance on SNR search

double TaylorF2ReducedSpinChirpTime(const double fStart, const double m1, const double m2, const double s1z, const double s2z, const int PNO)