util
General utilities (functions and structures) independent of modelling method
Defines
-
GWAT_LN2
natural log of 2
-
GWAT_G_SI
Gravitational constant in SI
-
GWAT_TWOPI
2*PI
-
GWAT_MSUN_SI
! Nominal solar mass, kg
-
GWAT_MTSUN_SI
Geometrized nominal solar mass, s !
!
GWAT_MTSUN_SI = GWAT_GMSUN_SI / c^3
-
GWAT_MRSUN_SI
Geometrized nominal solar mass, m !
!
GWAT_MRSUN_SI = GWAT_GMSUN_SI / c^3
-
PBSTR
-
PBWIDTH
-
MAX_TOL_ATAN
Functions
-
double gsl_LU_lndet(double **matrix, int dim)
-
int mvn_sample(int samples, double *mean, double **cov, int dim, double **output)
- Parameters:
output – [out] Size [samples][dim]
-
int mvn_sample(int samples, double *mean, double **cov, int dim, gsl_rng *r, double **output)
Samples from a multivariate-gaussian with covariance cov and mean mean, with dimension dim and puts the output in output Taken from blogpost: https://juanitorduz.github.io/multivariate_normal/.
- Parameters:
output – [out] Size [samples][dim]
-
void vector_union(std::vector<double> A, std::vector<double> B, std::vector<double> *C)
Finds union of two lists which have been pre-sorted. Returns the union, also sorted.
Allocates memory for C, needs to be deallocated
-
void matrix_multiply(double **A, double **B, double **C, int dim1, int dim2, int dim3)
Matrix multiplication.
Takes A — [dim1,dim2], B — [dim2,dim3], and produces C — [dim1,dim3]
Arrays should be pre-allocated
-
template<class T, class U>
void transform_parameters(gen_params_base<T> *param_in, gen_params_base<U> *param_out) Simple utility to copy the members of param_in to param_out, for whatever types those are.
This is sometimes required to jump from double params to adouble params
Memory MUST be allocated for betappe and bppe (if present), and this needs to be deallocated
-
template<class T, class U>
void transform_parameters(gen_params_base<T> *param_in, gen_params_base<U> **param_out) TESTING Simple utility to copy the members of param_in to param_out, for whatever types those are.
This is sometimes required to jump from double params to adouble params
Memory MUST be allocated for betappe and bppe (if present), and this needs to be deallocated
-
int newton_raphson_method_1d(void (*f_on_fprime)(double x, double *func, double *func_prime, void *param), double initial_guess, double tolerance, int max_iterations, void *parameters, double *solution)
Newton Raphson method for a function of one dimension.
f_on_fprime is a function pointer that should assign the value result = f(x) / f’(x)
param are just an extra needed parameters
return values:
0 -- success 1 -- Bad value for derivative 2 -- max iterations reached
-
template<class T>
T A0_from_DL(T chirpmass, T DL, bool sky_average) Transforms between chirpmass and DL to overall amplitude factor A0.
All quantities in seconds
-
template<class T>
T DL_from_A0(T chirpmass, T A0, bool sky_average) Transforms between amplitude factor A0 and chirpmass to DL.
All quantities in seconds
-
void initiate_LumD_Z_interp(gsl_interp_accel **Z_DL_accel_ptr, gsl_spline **Z_DL_spline_ptr)
Function that uses the GSL libraries to interpolate pre-calculated Z-D_L data.
Initiates the requried functions — GSL interpolation requires allocating memory before hand
-
void free_LumD_Z_interp(gsl_interp_accel **Z_DL_accel_ptr, gsl_spline **Z_DL_spline_ptr)
Frees the allocated interpolation function.
-
adouble Z_from_DL_interp(adouble DL, gsl_interp_accel *Z_DL_accel_ptr, gsl_spline *Z_DL_spline_ptr)
Function that returns Z from a given luminosity Distance — only Planck15
adouble version for ADOL-C calculations
-
double Z_from_DL_interp(double DL, gsl_interp_accel *Z_DL_accel_ptr, gsl_spline *Z_DL_spline_ptr)
Function that returns Z from a given luminosity Distance — only Planck15
-
double Z_from_DL(double DL, std::string cosmology)
Calculates the redshift given the luminosity distance.
Based on Astropy.cosmology calculations — see python script in the ./data folder of the project — numerically calculated given astropy.cosmology’s definitions (http://docs.astropy.org/en/stable/cosmology/) and used scipy.optimize to fit to a power series, stepping in half powers of DL. These coefficients are then output to a header file (D_Z_config.h) which are used here to calculate redshift. Custom cosmologies etc can easily be acheived by editing the python script D_Z_config.py, the c++ functions do not need modification. They use whatever data is available in the header file.
5 cosmological models are available (this argument must be spelled exactly, although case insensitive):
PLANCK15, PLANCK13, WMAP9, WMAP7, WMAP5
-
double DL_from_Z(double Z, std::string cosmology)
Calculates the luminosity distance given the redshift.
Based on Astropy.cosmology calculations — see python script in the ./data folder of the project — numerically calculated given astropy.cosmology’s definitions (http://docs.astropy.org/en/stable/cosmology/) and used scipy.optimize to fit to a power series, stepping in half powers of Z. These coefficients are then output to a header file (D_Z_config.h) which are used here to calculate distance. Custom cosmologies etc can easily be acheived by editing the python script D_Z_config.py, the c++ functions do not need modification. They use whatever data is available in the header file. If the functional form of the fitting function changes, these functions DO need to change.
5 cosmological models are available (this argument must be spelled exactly):
PLANCK15, PLANCK13, WMAP9, WMAP7, WMAP5
-
double cosmology_interpolation_function(double x, double *coeffs, int interp_degree)
Custom interpolation function used in the cosmology calculations.
Power series in half power increments of x, up to 11/2. powers of x
-
int cosmology_lookup(std::string cosmology)
Helper function for mapping cosmology name to an internal index.
-
template<class T>
T fcontact(T M_detector, T DL, std::string cosmology) Rough threshold for when to truncate waveforms for neutron stars.
24 km separation
- Parameters:
M_detector – Detector frame total mass in seconds
DL – Luminsosity distance in seconds
cosmology – Cosmology to use for getting the redshift from dl
-
double powerlaw_from_uniform(double x0, double x1, double power, double uniform_random_number)
Converts a uniformly random random number into a powerlaw distribution.
-
double gsl_maxwell_boltzmann_distribution(double sigma, gsl_rng *r)
Calculates a random number from the Maxwell-Boltzmann distribution using 3 gaussian random numbers.
Maxwell-Boltzmann describes a gaussian distributed vector in 3 space, so the magnitude is the RMS of the components of that vector
Assumed mean is at 0
Sigma would be equivalent to sqrt(kT/m) in the real maxwell boltzmann distribution
-
template<class T>
void list_intersect_ptrs(T **A, int lenA, T **B, int lenB, T **C, int *lenC) SHOULDN’T BE USED FOR DOUBLESCustom list-interesction implementation for sorted lists.
Uses pointers for efficiency — NOT COPIED — used std library for that
Takes A and B and puts the intersection in C — C = A /\ B
Takes the length of C and puts it in lenC
C should have allocated memory for at least as large as the smallest list
WILL NOT work in-place, A, B, and C must be separate
-
template<class T>
void list_intersect(T *A, int lenA, T *B, int lenB, T **C, int *lenC) Custom list-interesction implementation for sorted lists.
Uses pointers for efficiency — NOT COPIED — used std library for that
Takes A and B and puts the intersection in C — C = A /\ B
Takes the length of C and puts it in lenC
C should have allocated memory for at least as large as the smallest list
WILL NOT work in-place, A, B, and C must be separate
-
template<class T>
void list_intersect_value(T *A, int lenA, T *B, int lenB, T *C, int *lenC) Custom list-interesction implementation for sorted lists.
Uses pointers for efficiency — NOT COPIED — used std library for that
Takes A and B and puts the intersection in C — C = A /\ B
Takes the length of C and puts it in lenC
C should have allocated memory for at least as large as the smallest list
WILL NOT work in-place, A, B, and C must be separate
-
void rm_fisher_dim(double **input, int full_dim, double **output, int reduced_dim, int *removed_dims)
Removes a dimension from a matrix (made with fishers in mind, but this is general)
Note: the removed_dims list is indexed from 0
-
int gsl_cholesky_matrix_invert(double **input, double **inverse, int dim)
-
int normalized_gsl_cholesky_matrix_invert(double **input, double **inverse, int dim)
Normalize the Fisher matrix before inversion to try and tame singularity issues:
arXiv:1108.1826v2 Equation 60-61
-
adouble Z_from_DL(adouble DL, std::string cosmology)
Calculates the redshift given the luminosity distance adouble version for ADOL-C implementation.
-
adouble DL_from_Z(adouble Z, std::string cosmology)
Calculates the luminosity distance given the redshift adouble version for ADOL-C implementation.
-
adouble cosmology_interpolation_function(adouble x, double *coeffs, int interp_degree)
Custom interpolation function used in the cosmology calculations adouble version for ADOL-C.
-
double std_omega(double RA, double std_RA, double std_DEC, double cov_RA_DEC)
map the error on RA and DEC to error on solid angle \Omega
All quantities in rad
From arXiv:0906.4269
-
void printProgress(double percentage)
routine to print the progress of a process to the terminal as a progress bar
Call everytime you want the progress printed
-
void allocate_FFTW_mem_forward(fftw_outline *plan, int length)
Allocate memory for FFTW3 methods used in a lot of inner products input is a locally defined structure that houses all the pertinent data.
-
void allocate_FFTW_mem_reverse(fftw_outline *plan, int length)
Allocate memory for FFTW3 methods used in a lot of inner products —INVERSE input is a locally defined structure that houses all the pertinent data.
-
void deallocate_FFTW_mem(fftw_outline *plan)
deallocates the memory used for FFTW routines
-
double **allocate_2D_array(int dim1, int dim2)
Utility to malloc 2D array.
-
int **allocate_2D_array_int(int dim1, int dim2)
-
void deallocate_2D_array(double **array, int dim1, int dim2)
Utility to free malloc’d 2D array.
-
void deallocate_2D_array(int **array, int dim1, int dim2)
-
double ***allocate_3D_array(int dim1, int dim2, int dim3)
Utility to malloc 3D array.
-
int ***allocate_3D_array_int(int dim1, int dim2, int dim3)
-
void deallocate_3D_array(double ***array, int dim1, int dim2, int dim3)
Utility to free malloc’d 2D array.
-
void deallocate_3D_array(int ***array, int dim1, int dim2, int dim3)
Utility to free malloc’d 2D array.
-
void hann_window(double *window, int length)
-
void tukey_window(double *window, int length, double alpha)
Tukey window function for FFTs.
As defined by https://en.wikipedia.org/wiki/Window_function
-
template<class T>
void terr_pol_iota_from_equat_sph(T RA, T DEC, T thetaj, T phij, T *pol, T *iota) Routine to transform from the equatorial coordinate system spherical polar description of the total angular momentum to the detector specific polarization angle and inclination angle (of the total angular momentum)
For the terrestrial network, this calculates the polarization angle for a detector at the center of Earth in equatorial coordinates, which is the polarization angle used in the detector_response_equatorial functions
Polarization angle is defined as:
tan \psi = ( J.z - (J.N)(z.N) ) / (N.(Jxz))
Where z is the axis of rotation of earth (equatorial z), and N is the line of sight to the source (direction of propagation is -N)
The inclination angle is between the direction of propagation and the total angular momentum
- Parameters:
RA – Right ascension in rad
DEC – Declination in rad
thetaj – spherical polar angle of the total angular momentum in equatorial coordinates
phij – spherical azimuthal angle of the total angular momentum in equatorial coordinates
pol – polarization defined by tan \psi = ( J.z - (J.N)(z.N) ) / (N.(Jxz))
iota – Inclination angle of the TOTAL angular momentum and the direction of propagation -N
-
template<class T>
void ecl_from_eq(T theta_eq, T phi_eq, T *theta_ecl, T *phi_ecl) transform spherical angles from equatorial to ecliptic
Rotation about the vernal equinox (x-hat in both coordinate systems) by the axial tilt or obliquity
From wikipedia
Compared against astropy and agrees
- Parameters:
theta_eq – Equatorial spherical polar angle
phi_eq – Equatorial spherical azimuthal angle
theta_ecl – Ecliptic spherical polar angle
phi_ecl – Ecliptic spherical Azimuthal angle
-
template<class T>
void equatorial_from_SF(T *SFvec, T thetal, T phil, T thetas, T phis, T iota, T phi_ref, T *EQvec) Utility to map source frame vectors to equatorial frame vectors.
Needs the equatorial vectors for the line of sight and the spin angular momentum
Needs the inclination angle and the reference phase
Works by constructing a third vector from the cross product of the two others, in both frames. This defines a family of 3 vectors that can uniquely determine a rotation matrix from source frame to equatorial, which was done analytically in mathematica ( see the nb anglenb.nb)
K = LxN
A = {LSF,NSF, KSF}
B={LEQ,NEQ,KEQ}
B = R.A
R = B.A^-1
Verified with specific cases and with dot products with vectors before and after rotation
- Parameters:
SFvec – Input, source frame vector, as defined by LAL (L = z_hat) in cartesian
thetal – Polar angle in equatorial coordinates of the orbital angular momentum at the reference frequency
phil – Azimuthal angle in equatorial coordinates of the orbital angular momentum at the reference frequency
thetas – Polar angle in equatorial coordinates of the source (not direction of propagation)
phis – Azimuthal angle in equatorial coordinates of the source (not direction of propagation)
iota – Inclination angle between the orbiatl angular momentum and the direction of propagation at the reference frequency
phi_ref – Reference frequency for the waveform (defines the LAL frame ) ( orbital phase) at the reference frequency
EQvec – [out] Out vector in equatorial coordinates in cartesian
-
double calculate_eta(double mass1, double mass2)
Calculates the symmetric mass ration from the two component masses.
-
adouble calculate_eta(adouble mass1, adouble mass2)
-
double calculate_chirpmass(double mass1, double mass2)
Calculates the chirp mass from the two component masses.
The output units are whatever units the input masses are
-
adouble calculate_chirpmass(adouble mass1, adouble mass2)
-
double calculate_mass1(double chirpmass, double eta)
Calculates the larger mass given a chirp mass and symmetric mass ratio.
Units of the output match the units of the input chirp mass
-
adouble calculate_mass1(adouble chirpmass, adouble eta)
-
double calculate_mass2(double chirpmass, double eta)
Calculates the smaller mass given a chirp mass and symmetric mass ratio.
Units of the output match the units of the input chirp mass
-
adouble calculate_mass2(adouble chirpmass, adouble eta)
-
double calculate_mass1_Mcq(double chirpmass, double q)
Calculates the larger mass given a chirp mass and asymmetric mass ratio.
Units of the output match the units of the input chirp mass
-
adouble calculate_mass1_Mcq(adouble chirpmass, adouble q)
-
double calculate_mass2_Mcq(double chirpmass, double q)
Calculates the smaller mass given a chirp mass and asymmetric mass ratio.
Units of the output match the units of the input chirp mass
-
adouble calculate_mass2_Mcq(adouble chirpmass, adouble q)
-
template<class T>
void celestial_horizon_transform(T RA, T DEC, double gps_time, T LONG, T LAT, T *phi, T *theta) Utility to transform from celestial coord RA and DEC to local horizon coord for detector response functions.
Outputs are the spherical polar angles defined by North as 0 degrees azimuth and the normal to the earth as 0 degree polar
- Parameters:
RA – Right acsension (rad)
DEC – Declination (rad)
gps_time – GPS time
LONG – Longitude (rad)
LAT – Latitude (rad)
phi – [out] horizon azimuthal angle (rad)
theta – [out] horizon polar angle (rad)
-
template<class T>
T gps_to_GMST(T gps_time) Utility to transform from gps time to GMST https://aa.usno.navy.mil/faq/docs/GAST.php.
-
void decimal_to_HMS(double decimal, int *hour, int *min, double *second)
Converts decimal into hours, minutes and seconds.
I just got tired of keeping all these steps straight..
-
template<class T>
void transform_cart_sph(T *cartvec, T *sphvec) utility to transform a vector from cartesian to spherical (radian)
order:
cart: x, y, z
spherical: r, polar, azimuthal
Tested in all octants
-
template<class T>
void transform_sph_cart(T *sphvec, T *cartvec) utility to transform a vector from spherical (radian) to cartesian
Range: \theta \el [0,\pi]
\phi \el [0,2 \pi]
order:
cart: x, y, z
spherical: r, polar, azimuthal
Tested in all octants
-
template<class T>
void unwrap_array(T *in, T *out, int len) Unwrap angles from arctan.
Stolen from stack exchange.. https://stackoverflow.com/questions/15634400/continous-angles-in-c-eq-unwrap-function-in-matlab
Modified
-
template<class T>
T trapezoidal_sum_uniform(double delta_x, int length, T *integrand) Trapezoidal sum rule to approximate discrete integral - Uniform spacing.
This version is faster than the general version, as it has half the function calls
Something may be wrong with this function - had an overall offset for real data that was fixed by using the simpsons rule - not sure if this was because of a boost in accuracy or because something is off with the trapezoidal sum
-
template<class T>
T trapezoidal_sum(double *delta_x, int length, T *integrand) Trapezoidal sum rule to approximate discrete integral - Non-Uniform spacing.
This version is slower than the uniform version, but will handle non-uniform spacing
-
template<class T>
T simpsons_sum(double delta_x, int length, T *integrand) Simpsons sum rule to approximate discrete integral - Uniform spacing.
More accurate than the trapezoidal rule, but must be uniform
-
long factorial(long num)
Local function to calculate a factorial.
-
double pow_int(double base, int power)
Local power function, specifically for integer powers.
Much faster than the std version, because this is only for integer powers
-
adouble pow_int(adouble base, int power)
-
template<class T>
std::complex<T> XLALSpinWeightedSphericalHarmonic(T theta, T phi, int s, int l, int m) Shamelessly stolen from LALsuite
- Parameters:
theta – polar angle (rad)
phi – azimuthal angle (rad)
s – spin weight
l – mode number l
m – mode number m
-
double cbrt_internal(double base)
Fucntion that just returns the cuberoot.
-
adouble cbrt_internal(adouble base)
Fucntion that just returns the cuberoot ADOL-C doesn’t have the cbrt function (which is faster), so have to use the power function.
Variables
-
const double gamma_E = 0.5772156649015328606065120900824024310421
Euler number
-
const double c = 299792458.
Speed of light m/s
-
const double G = GWAT_G_SI * (1.988409902147041637325262574352366540e30)
Gravitational constant in m**3/(s**2 SolMass)
-
const double MSOL_SEC = 4.925491025543575903411922162094833998e-6
G/c**3 seconds per solar mass
-
const double GWAT_MPC_SI = 3.085677581491367278913937957796471611e22
Meters in a megaparsec
-
const double MPC_SEC = GWAT_MPC_SI / c
consts.kpc.to(‘m’)*1000/c Mpc in sec
-
const double T_year = 31557600.
1 year in seconds — ie seconds/year
-
const double T_day = 24 * 3600.
1 day in seconds
-
const double AXIAL_TILT = 0.409092627749
Earth’s Axial tilt in radian
-
const double AU_SEC = 499.005
1AU in seconds
-
const double h_planck = 4.135667696e-15
Planck’s constant (unbarred) in ev sec
-
const double DOUBLE_COMP_THRESH = 1e-10
-
struct fftw_outline
- #include <util.h>
-
template<class T>
struct sph_harm - #include <util.h>
Container for set of spin weighted spherical harmonics (all spins are -2)
-
template<class T>
class gen_params_base - #include <util.h>
Structure for interfacing with the libraries.
Structure to interface with the libraries - Units are in solar masses and mpc contains the generation parameters including source parameters and theory parameters
NOTE not all the members of this structure need to be assigned for usage. In fact, some are reduntant. It’s up to the user to determine what fields require an assignment. (Sorry)
Public Functions
-
inline void print_properties()
Public Members
-
bool tidal_love = true
Boolean flag indicating binary love relations should be used
-
bool tidal_love_error = false
Boolean flag indicating error marginalization over residual EoS dependence in binary love relations should be performed.
-
bool equatorial_orientation = false
boolean flag indicating equatorial orientation coordinates should be used
-
bool horizon_coord = false
Boolean flag indicating local, horizon coordinates should be used
-
double gmst
Greenwich Mean Sidereal time (for detector orientation - start of data — IN RADIANS NOT HOURS
-
bool NSflag1 = false
BOOL flag for early termination of NS binaries
-
bool NSflag2 = false
-
bool dep_postmerger = false
Flag to force the use of deprecated postmerger calculations — ADOLC friendly
-
bool shift_time = true
Shift time detemines if times are shifted so coalescence is more accurately
-
bool shift_phase = true
Shift time detemines if phic or phiRef is used
-
bool sky_average = false
-
bool precess_reduced_flag = false
-
T theta_j_ecl
Polar angle in ecliptic coordinates
Azimuthal angle in ecliptic coordinates for the total angular momentum — internal, and should not be specified
-
int Nmod_beta = 0
-
int Nmod_alpha = 0
-
int Nmod_sigma = 0
-
int Nmod_phi = 0
-
int *betai = NULL
-
int *alphai = NULL
-
int *sigmai = NULL
-
int *phii = NULL
-
double *bppe = NULL
ppE b parameter (power of the frequency) - vector for multiple modifications
-
double *appe = NULL
-
int Nmod = 0
Number of phase modificatinos
-
int PNorder = 35
-
bool EA_region1 = false
-
bool alpha_param = true
-
bool include_l1 = false
-
gsl_spline *Z_DL_spline_ptr = NULL
-
gsl_interp_accel *Z_DL_accel_ptr = NULL
-
inline void print_properties()
-
class gen_params : public gen_params_base<double>
- #include <util.h>
convience wrapper for the gen_params_base class
If using the code in the intended way, this is all the user should ever have to use. Just allows the user to drop the template parameter
Also implemented for backwards compatibility with previous versions of the code
-
template<class T>
struct useful_powers - #include <util.h>
To speed up calculations within the for loops, we pre-calculate reoccuring powers of M*F and Pi, since the pow() function is prohibatively slow.
Powers of PI are initialized once, and powers of MF need to be calculated once per for loop (if in the inspiral portion).
use the functions precalc_powers_ins_amp, precalc_powers_ins_phase, precalc_powers_pi to initialize
-
template<class T>
struct source_parameters - #include <util.h>
For internal data transfers.
Structure to facililate parameter tranfers - All dimensionful quantities are in seconds
Public Functions
-
void populate_source_parameters(gen_params_base<T> *param_in)
Builds the structure that shuttles source parameters between functions -updated version to incorporate structure argument.
Populates the structure that is passed to all generation methods - contains all relavent source parameters
Template type of source parameters and gen_parameters must match
Public Members
-
bool shift_phase = true
Shift time detemines if phic or phiRef is used
-
bool dep_postmerger = false
Flag to force the use of deprecated postmerger calculations — ADOLC friendly
-
bool NSflag1
-
bool NSflag2
-
int PNorder = 35
-
double *bppe
-
double *appe
-
int Nmod
Number of modifications to phase
-
bool sky_average
-
bool shift_time = true
Boolean — shift time to 0 before shifting to tc or not
-
gsl_spline *Z_DL_spline_ptr = NULL
-
gsl_interp_accel *Z_DL_accel_ptr = NULL
-
int Nmod_beta = 0
-
int Nmod_alpha = 0
-
int Nmod_sigma = 0
-
int Nmod_phi = 0
-
int *betai
-
int *alphai
-
int *sigmai
-
int *phii
-
bool include_l1 = false
-
bool EA_nan_error_message = false
-
bool alpha_param = true
-
bool EA_region1 = false
Public Static Functions
-
static source_parameters<T> populate_source_parameters_old(T mass1, T mass2, T Luminosity_Distance, T *spin1, T *spin2, T phi_c, T t_c, bool sky_average)
Builds the structure that shuttles source parameters between functions- outdated in favor of structure argument.
Populates the structure that is passed to all generation methods - contains all relavent source parameters
- Parameters:
mass1 – mass of the larger body - in Solar Masses
mass2 – mass of the smaller body - in Solar Masses
Luminosity_Distance – Luminosity Distance in Mpc
spin2 – spin vector of the larger body {sx,sy,sz}
phi_c – spin vector of the smaller body {sx,sy,sz}
t_c – coalescence phase
sky_average – coalescence time
-
void populate_source_parameters(gen_params_base<T> *param_in)