detector_util

Header file for all detector-specific utilities

Functions

void populate_noise(double *frequencies, std::string detector, double *noise_root, int length, double integration_time = 48)

Function to populate the squareroot of the noise curve for various detectors.

If frequencies are left as NULL, standard frequency spacing is applied and the frequencies are returned, in which case the frequencies argument becomes an output array

Detector names must be spelled exactly

Detectors include:

aLIGO_analytic -- analytic approximation of the advanced LIGO sensitivity curve

Hanford_O1_fitted -- Fitted function to the O1 noise curve for Hanford

AdLIGOMidHigh -- from Emanuele Berti

LISA -- LISA sensitivity curve with out sky averaging for a single channel 

    -- does NOT include the factor of sqrt(3)/2, which is taken care of in the response function for LISA

LISA_CONF -- LISA sensitivity curve with out sky averaging for a single channel including confusion noise from the galatic white dwarf population

    -- does NOT include the factor of sqrt(3)/2, which is taken care of in the response function for LISA

LISA_SADC -- LISA sensitivity curve with sky averaging for a dual channel 

    -- does include the factor of sqrt(3)/2

    -- This DOES include the sky averaged pattern functions -- i.e. sky_average flag should NOT be used when making the waveform

LISA_SADC_CONF -- LISA sensitivity curve with sky averaging for a dual channel including confusion noise from the galatic white dwarf population

    -- does include the factor of sqrt(3)/2

    -- This DOES include the sky averaged pattern functions -- i.e. sky_average flag should NOT be used when making the waveform

AdLIGOAPlus -- AdLIGO A+ O5

AdLIGOAPlus_smoothed -- AdLIGO A+ O5 -- without spectral lines

AdLIGODesign -- Design sensitivity for LIGO O4

AdLIGODesign_smoothed -- Design sensitivity for LIGO O4 -- without spectral lines

CE1 -- Cosmic Explorer phase 1 -- circa 2035

CE2 -- Cosmic Explorer phase 2 -- circa 2045

KAGRA_opt -- KAGRA O5 05 optimistic

KAGRA_pess -- KAGRA O5 05 pessimistic

AdVIRGOPlus2_opt -- phase 2, O5, optimistic

AdVIRGOPlus2_opt_smoothed -- phase 2, O5, optimistic -- without spectral lines

AdVIRGOPlus2_pess -- phase 2, O5, pessimistic

AdVIRGOPlus2_pess_smoothed -- phase 2, O5, pessimistic -- without spectral lines

AdVIRGOPlus1 -- phase 1, O4 analog ( rescaled O5 curve according to BBH distance estimates in https://dcc.ligo.org/public/0161/P1900218/002/SummaryForObservers.pdf)

AdVIRGOPlus1_smoothed -- phase 1, O4 analog ( rescaled O5 curve according to BBH distance estimates in https://dcc.ligo.org/public/0161/P1900218/002/SummaryForObservers.pdf) -- without spectral lines

ET-D -- Einstein telescope with lines included. From ET website

ET-D_smooth -- Einstein telescope without lines. From Emanuele 

AdLIGOVoyager -- Voyager configuration of Advanced LIGO

Parameters:

frequencies – double array of frquencies (NULL)
detector – String to designate the detector noise curve to be used
noise_root – ouptput double array for the square root of the PSD of the noise of the specified detector
length – integer length of the output and input arrays
integration_time – Integration time in months (only important for LISA_conf

double aLIGO_analytic(double f)

Analytic function approximating the PSD for aLIGO.

CITE (Will?)

double LISA_analytic_SADC(double f)

Analytic function approximating the PSD for LISA sensitivity curve — this is S, not root S.

Sky averaged - dual channel

arXiv:1803.01944 — equation 13

NOTE: may need to divide by \sqrt{3}/2.. My LISA response functions already include this factor — factor 0f 3./4. (sqrt(3)/2)**2

Keeping the factor of sqrt(3)/2, since this will probably be used without the beam patter functions (that’s where the factor of sqrt(3)/2 usually appears)

double LISA_analytic(double f)

Analytic function approximating the PSD for LISA sensitivity curve — this is S, not root S and does not sky average, treats the 2 channels separately, and the geometrical factor of \sqrt{3}/2 is included in the waveform.

NON Sky averaged - single channel

arXiv:1905.08811 — equation 14-17

double LISA_POMS(double f)

\breif Optical metrology noise function

arXiv:1803.01944 — equation 10

double LISA_PACC(double f)

\breif Single test mass accelartion noise

arXiv:1803.01944 — equation 11

double LISA_SC(double f, double alpha, double beta, double kappa, double gamma, double fk)

\breif Confusion noise — this is S, not root S

arXiv:1803.01944 — 14

integration_time is the observation time in months — options are 6, 12, 24, and 48

void sort_LISA_SC_coeffs(double *alpha, double *beta, double *kappa, double *gamma, double *fk, double integration_time)

\breif LISA confusion noise coefficients

arXiv:1803.01944 — Table 1

integration_time is the observation time in months — options are 6, 12, 24, and 48

std::complex<double> Q(double theta, double phi, double iota)

Utility for the overall amplitude and phase shift for spin-aligned systems.

For spin aligned, all the extrinsic parameters have the effect of an overall amplitude modulation and phase shift

std::complex<double> Q(double theta, double phi, double iota, double psi)

Utility for the overall amplitude and phase shift for spin-aligned systems.

For spin aligned, all the extrinsic parameters have the effect of an overall amplitude modulation and phase shift

template<class T> T right_interferometer_l(T theta, T phi, T psi)

Response function of a 90 deg interferometer for longitudinal (scalar) polarization.

Theta, phi, and psi are local, horizontal coordinates relative to the detector

template<class T> T right_interferometer_b(T theta, T phi, T psi)

Response function of a 90 deg interferometer for breathing (scalar) polarization.

Theta, phi, and psi are local, horizontal coordinates relative to the detector

template<class T> T right_interferometer_y(T theta, T phi, T psi)

Response function of a 90 deg interferometer for y (vector) polarization.

Theta, phi, and psi are local, horizontal coordinates relative to the detector

template<class T> T right_interferometer_x(T theta, T phi, T psi)

Response function of a 90 deg interferometer for x (vector) polarization.

Theta, phi, and psi are local, horizontal coordinates relative to the detector

template<class T> T right_interferometer_cross(T theta, T phi)

Response function of a 90 deg interferometer for cross polarization.

Theta and phi are local, horizontal coordinates relative to the detector

template<class T> T right_interferometer_plus(T theta, T phi)

Response function of a 90 deg interferometer for plus polarization.

Theta and phi are local, horizontal coordinates relative to the detector

template<class T> void right_interferometer(det_res_pat<T> *r_pat, T theta, T phi, T psi)

Response function of a 90 deg interferometer.

Theta and phi are local, horizontal coordinates relative to the detector

With general psi

double Hanford_O1_fitted(double f)

Numerically fit PSD to the Hanford Detector’s O1.

CITE (Yunes?)

template<class T> void celestial_horizon_transform(T RA, T DEC, double gps_time, std::string detector, T *phi, T *theta)

Transform from celestial coordinates to local horizontal coords.

(RA,DEC) -> (altitude, azimuth)

Need gps_time of transformation, as the horizontal coords change in time

detector is used to specify the lat and long of the local frame

Parameters:

RA – in RAD
DEC – in RAD
phi – in RAD
theta – in RAD

void derivative_celestial_horizon_transform(double RA, double DEC, double gps_time, std::string detector, double *dphi_dRA, double *dtheta_dRA, double *dphi_dDEC, double *dtheta_dDEC)

Numerical derivative of the transformation.

Planned for use in Fisher calculations, but not currently implemented anywhere

Parameters:

RA – in RAD
DEC – in RAD

double DTOA(double theta1, double theta2, std::string detector1, std::string detector2)

calculate difference in time of arrival (DTOA) for a given source location and 2 different detectors

Parameters:

theta1 – spherical polar angle for detector 1 in RAD
theta2 – spherical polar angle for detector 2 in RAD
detector1 – name of detector one
detector2 – name of detector two

double radius_at_lat(double latitude, double elevation)

/brief Analytic approximation of the radius from the center of earth to a given location

Just the raidus as a function of angles, modelling an oblate spheroid

Parameters:

latitude – latitude in degrees
elevation – elevation in meters

template<class T> void detector_response_functions_equatorial(double D[3][3], T ra, T dec, T psi, double gmst, det_res_pat<T> *r_pat)

Calculates the response coefficients for a detector with response tensor D for a source at RA, Dec, and psi.

Taken from LALSuite

The response tensor for each of the operational detectors is precomputed in detector_util.h, but to create a new tensor, follow the outline in Anderson et al 36 PRD 63 042003 (2001) Appendix B

For terrestial detectors — psi is the polarization angle from a detector at earth’s center, aligned with equatorial coordinates

Parameters:

D – Detector Response tensor (3x3)
ra – Right ascension in rad
dec – Declination in rad
psi – polarization angle in rad
gmst – Greenwich mean sidereal time (rad)

template<class T> void detector_response_functions_equatorial(std::string detector, T ra, T dec, T psi, double gmst, det_res_pat<T> *r_pat)

Wrapping of the equatorial detector response for terrestial based detectors.

For ground based detectors, the antenna pattern functions are not functions of time.

Parameters:

detector – Detector
ra – Right ascension in rad
dec – Declination in rad
psi – polarization angle in rad
gmst – Greenwich mean sidereal time (rad)

template<class T> void detector_response_functions_equatorial(std::string detector, T ra, T dec, T psi, double gmst, T *times, int length, T LISA_alpha0, T LISA_phi0, T theta_j_ecl, T phi_j_ecl, det_res_pat<T> *r_pat)

Same as the other function, but for active and future detectors.

If terrestial detectors, it will use ra, dec, psi, and gmst

If space detector, it will use ra, dec transformed to ecliptic coord, and LISA_alpha0, LISA_phi0 and theta_j_ecl, phi_j_ecl, which are the offsets for LISA and the spherical angles for the unit vector in the total angular momemntum in the ecliptic system

Parameters:

detector – Detector
ra – Right ascension in rad
dec – Declination in rad
psi – polarization angle in rad
gmst – Greenwich mean sidereal time (rad)
times – Times at which to evaluate Fplus and Fcross, in which case the Fplus and Fcross pointers are arrays
length – Length of the arrays
LISA_alpha0 – Offset for alpha
LISA_phi0 – Offset for phi
theta_j_ecl – Ecliptic spherical polar angle of Lhat (Jhat for precessing systems)
phi_j_ecl – Ecliptic sherical azimuthal angle (Jhat for precessing systems)

template<class T> T DTOA_DETECTOR(T RA, T DEC, double GMST_rad, std::string detector1, std::string detector2)

calculate difference in time of arrival (DTOA) for a given source location and 2 different detectors

Full version, from LAL

Parameters:

RA – spherical polar angle for detector 1 in RAD
DEC – spherical polar angle for detector 2 in RAD
GMST_rad – Greenwich mean sidereal time of detection
detector1 – name of detector one
detector2 – name of detector two

template<class T> T DTOA_earth_centered_coord(T RA, T DEC, double GMST_rad, const double *loc1, const double *loc2)

calculate difference in time of arrival (DTOA) for a given source location and 2 different detectors

Parameters:

RA – spherical polar angle for detector 1 in RAD
DEC – spherical polar angle for detector 2 in RAD
GMST_rad – Greenwich mean sidereal time of detection
loc1 – Location of the first detector in Earth centered coordinates in meters
loc2 – Location of the second detector in Earth centered coordinates in meters

template<class T> T LISA_response_plus_time(T theta_s, T phi_s, T theta_j, T phi_j, T alpha_0, T phi_0, T t)

Time dependent detector response of LISA for non-precessing waveforms.

See https://arxiv.org/abs/gr-qc/0411129 or https://arxiv.org/abs/gr-qc/9703068

All the arguments are `barred’, using the notation in these two works. That is, they are relative to the solar system barycenter.

To get the second interferometer’s response, evaluate with phi_j - pi/4.

All the coordinates are in the ecliptic coordinate system

template<class T> T LISA_response_cross_time(T theta_s, T phi_s, T theta_j, T phi_j, T alpha_0, T phi_0, T t)

template<class T> T LISA_response_plus(source_parameters<T> *params, T theta_s, T phi_s, T theta_j, T phi_j, T alpha_0, T phi_0, T f)

template<class T> T LISA_response_cross(source_parameters<T> *params, T theta_s, T phi_s, T theta_j, T phi_j, T alpha_0, T phi_0, T f)

double p_single_detector(double omega, int samples)

Utility to calculate the cumulative amplitude distribution for a single detector.

P(\omega) = \int_V \Theta(\omega’(\Omega, \psi, \iota)-\omega) d\Omega d\psi dcos\iota

Integrated over the volume which \omega’ is larger than \omega

Integrates using Monte Carlo integration

Uniform sampling in \psi, cos(\iota), cos(\theta) , and \phi

Parameters:

omega – \omega = \rho/\rho_opt*
samples – number of monte carlo samples to use*

double p_N_detector(double omega, int samples, int N_detectors, std::string *detectors, int rand_seed)

Utility to calculate the cumulative amplitude distribution for three detectors.

P(\omega) = \int_V \Theta(\omega’(\Omega, \psi, \iota)-\omega) d\Omega d\psi dcos\iota

Integrated over the volume which \omega’ is larger than \omega

Integrates using Monte Carlo integration

Uniform sampling in \psi, cos(\iota), cos(DEC) , and RA

Sampled using any supported detectors. See documentation for supported detectors

Rough upper limit: HLV: 1.4 HLVK: 1.5 HLVI: 1.65 HLVKI: 1.7

Parameters:

omega – \omega = \rho/\rho_opt*
samples – number of monte carlo samples to use*
N_detectors – Number of detectors
detectors – String name of detectors to use — MUST BE SUPPORTED DETECTORS
rand_seed – Seed for random number generator

double p_single_detector_fit(double omega)

Utility to calculate the cumulative amplitude distribution for a single detector — Numerical Fit.

P(\omega) = \int_V \Theta(\omega’(\Omega, \psi, \iota)-\omega) d\Omega d\psi dcos\iota

Integrated over the volume which \omega’ is larger than \omega

see arXiv:1405.7016

Parameters:: omega – \omega = \rho/\rho_opt*

double p_triple_detector_fit(double omega)

Utility to calculate the cumulative amplitude distribution for triple detector network — Numerical Fit.

P(\omega) = \int_V \Theta(\omega’(\Omega, \psi, \iota)-\omega) d\Omega d\psi dcos\iota

Integrated over the volume which \omega’ is larger than \omega

see arXiv:1405.7016

Not accurate, apparently. Just use interpolation

Parameters:: omega – \omega = \rho/\rho_opt*

double pdet_triple_detector_fit(double rho_thresh, double rho_opt)

double p_triple_detector_interp(double omega)

Utility to calculate the cumulative amplitude distribution for triple detector network — interpolated data from https://pages.jh.edu/~eberti2/research/.

P(\omega) = \int_V \Theta(\omega’(\Omega, \psi, \iota)-\omega) d\Omega d\psi dcos\iota

Integrated over the volume which \omega’ is larger than \omega

see arXiv:1405.7016

Parameters:: omega – \omega = \rho/\rho_opt*

double p_single_detector_interp(double omega)

Utility to calculate the cumulative amplitude distribution for single detector network — interpolated data from https://pages.jh.edu/~eberti2/research/.

P(\omega) = \int_V \Theta(\omega’(\Omega, \psi, \iota)-\omega) d\Omega d\psi dcos\iota

Integrated over the volume which \omega’ is larger than \omega

see arXiv:1405.7016

Parameters:: omega – \omega = \rho/\rho_opt*

Variables

const double H_LAT = 0.81079526383

const double H_LONG = -2.08405676917

const double H_azimuth_offset = 2.199

const double H_radius = 6367299.93401105

const double H_elevation = 142.554

const double H_location[3] = {-2.16141492636e+06, -3.83469517889e+06, 4.60035022664e+06}

const double H_geometric_factor = 1.

const double L_LAT = 0.53342313506

const double L_LONG = -1.58430937078

const double L_azimuth_offset = 3.4557

const double L_radius = 6372795.50144497

const double L_elevation = -6.574

const double L_location[3] = {-7.42760447238e+04, -5.49628371971e+06, 3.22425701744e+06}

const double L_geometric_factor = 1.

const double V_LAT = 0.76151183984

const double V_LONG = 0.18333805213

const double V_azimuth_offset = 1.239

const double V_radius = 6368051.92301

const double V_elevation = 51.884

const double V_location[3] = {4.54637409900e+06, 8.42989697626e+05, 4.37857696241e+06}

const double V_geometric_factor = 1.

const double K_LAT = 0.6355068497

const double K_LONG = 2.396441015

const double K_radius = 6371060

const double K_azimuth_offset = 1.239

const double K_elevation = 414.181

const double K_location[3] = {-3777336.024, 3484898.411, 3765313.697}

const double K_geometric_factor = 1.

const double I_LAT = 0.2484185302005262

const double I_LONG = 1.3340133249409993

const double I_radius = 6376850

const double I_azimuth_offset = 1.239

const double I_elevation = 3.9624

const double I_location[3] = {1450526.82294155, 6011058.39047265, 1558018.27884102}

const double I_geometric_factor = 1.

const double CE_LAT = 0.706553

const double CE_LONG = -1.99869

const double CE_radius = 6370201

const double CE_azimuth_offset = 0

const double CE_elevation = 1190.854

const double CE_location[3] = {-2.01262e6, -4.41298e6, 4.13995e6}

const double CE_geometric_factor = 1.

const double ET1_LAT = 0.76151183984

const double ET1_LONG = 0.18333805213

const double ET1_radius = 6.368051923006691e6

const double ET1_azimuth_offset = 0.33916285222

const double ET1_elevation = 51.884

const double ET1_location[3] = {4.54637409900e+06, 8.42989697626e+05, 4.37857696241e+06}

const double ET1_geometric_factor = 0.86602540378

const double ET2_LAT = 0.76299307990

const double ET2_LONG = 0.18405858870

const double ET2_radius = 6.36802814516601e6

const double ET2_azimuth_offset = 0

const double ET2_elevation = 59.735

const double ET2_location[3] = {4.53936951685e+06, 8.45074592488e+05, 4.38540257904e+06}

const double ET2_geometric_factor = 0.86602540378

const double ET3_LAT = 0.76270463257

const double ET3_LONG = 0.18192996730

const double ET3_radius = 6.368034296190262e6

const double ET3_azimuth_offset = 0

const double ET3_elevation = 59.727

const double ET3_location[3] = {4.54240595075e+06, 8.35639650438e+05, 4.38407519902e+06}

const double ET3_geometric_factor = 0.86602540378

const double RE_polar = 6357e3

const double RE_equatorial = 6378e3

const double Hanford_D[3][3] = {{-0.392614, -0.0776134, -0.247389}, {-0.0776134, 0.319524, 0.227998}, {-0.247389, 0.227998, 0.07309}}

Response Tensor for Hanford

Calculated using arXiv:gr-qc/0008066, equation B6 with the table at the end — see mathematica script

const double Livingston_D[3][3] = {{0.411281, 0.14021, 0.247295}, {0.14021, -0.109006, -0.181616}, {0.247295, -0.181616, -0.302275}}

Response Tensor for Livingston

Calculated using arXiv:gr-qc/0008066, equation B6 with the table at the end — see mathematica script

const double Virgo_D[3][3] = {{0.243874, -0.0990838, -0.232576}, {-0.0990838, -0.447826, 0.187833}, {-0.232576, 0.187833, 0.203952}}

Response Tensor for Virgo

Calculated using arXiv:gr-qc/0008066, equation B6 with the table at the end — see mathematica script

const double Kagra_D[3][3] = {{-0.18599, 0.153167, -0.324951}, {0.153167, 0.349518, -0.170874}, {-0.324951, -0.170874, -0.163529}}

const double Indigo_D[3][3] = {{0.470823664959826499, -0.120908243530907122, 0.0279526040438164719}, {-0.120908243530907122, -0.001050025852002534, 0.115836952558600455}, {0.0279526040438164719, 0.115836952558600455, -0.469773644459096742}}

const double CE_D[3][3] = {{0.377622, -0.268333, -0.102451}, {-0.268333, -0.0883621, -0.224639}, {-0.102451, -0.224639, -0.289259}}: Response Tensor for CE placed at Great Basin desert with the yarm pointing north and the x arm pointing east

const double ET1_D[3][3] = {{0.0878588, -0.36468, -0.0208748}, {-0.36468, -0.518498, 0.475276}, {-0.0208748, 0.475276, -0.0693608}}

const double ET2_D[3][3] = {{-0.444542, 0.00279528, 0.457953}, {0.00279528, 0.401623, -0.0796882}, {0.457953, -0.0796882, -0.457081}}

const double ET3_D[3][3] = {{-0.0134683, 0.432316, -0.0687841}, {0.432316, -0.620065, -0.327299}, {-0.0687841, -0.327299, 0.133534}}

const int analytic_PSD_models_N = 6

const std::string analytic_PSD_models[6] = {"aLIGO_analytic", "Hanford_O1_fitted", "LISA_SADC", "LISA_SADC_CONF", "LISA_CONF", "LISA"}

const int interp_PSD_models_N = 23

const std::string interp_PSD_models[23] = {"AdLIGOMidHigh", "_AdLIGODesign", "AdLIGODesign", "AdLIGODesign_smoothed", "AdLIGOAPlus", "AdLIGOAPlus_smoothed", "CE1", "CE1_smoothed", "CE2", "CE2_smoothed", "AdVIRGOPlus2_opt", "AdVIRGOPlus2_opt_smoothed", "AdVIRGOPlus2_pess", "AdVIRGOPlus2_pess_smoothed", "AdVIRGOPlus1", "AdVIRGOPlus1_smoothed", "KAGRA_opt", "KAGRA_opt_smoothed", "KAGRA_pess", "KAGRA_pess_smoothed", "ET-D", "ET-D_smoothed", "AdLIGOVoyager"}

template<class T> struct det_res_pat

#include <detector_util.h>

Public Members

T *Fplus = NULL

T *Fcross = NULL

T *Fx = NULL

T *Fy = NULL

T *Fb = NULL

T *Fl = NULL

bool *active_polarizations