Commit 6d86b4ad authored by Jonathan Blackman's avatar Jonathan Blackman

Refactor NRSur7dq2 and fix for waveform review

-There are now two approximants:
    NRSur7dq2_HM uses all modes available (ell <= 4)
    NRSur7dq2_L2only uses only the (2, +/-2) modes in the co-orbital frame
-Use lalsimulation lalinference_o2 conventions in ChooseTDWaveform
 As there is not a consistent convention, we follow SEOBNRv3.
-Fix various issues with GPS time, fmin, fref
-Add StartFrequency function giving smallest allowed fmin and fref
-Add documentation
-Fix small numeric error in waveform evaluation
parent cd767089
......@@ -48,6 +48,7 @@ extern "C" {
* @defgroup LALSimIMRSEOBNRROM_c LALSimIMRSEOBNRvxROMXXX.c
* @defgroup LALSimIMRSEOBNRv2ChirpTime_c LALSimIMRSEOBNRv2ChirpTime.c
* @defgroup LALSimIMRPSpinInspiralRD_c LALSimIMRPSpinInspiralRD.c
* @defgroup LALSimNRSur7dq2_c LALSimIMRNRSur7dq2.c
* @}
*
* @addtogroup LALSimIMR_h
......@@ -247,6 +248,55 @@ int XLALSimInspiralNRWaveformGetHplusHcross(
const char *NRDataFile /**< Location of NR HDF file */
);
/* in module LALSimIMRNRSur7dq2.c */
double XLALSimInspiralNRSur7dq2StartFrequency(
REAL8 m1, /**< mass of companion 1 (kg) */
REAL8 m2, /**< mass of companion 2 (kg) */
REAL8 s1x, /**< initial value of S1x */
REAL8 s1y, /**< initial value of S1y */
REAL8 s1z, /**< initial value of S1z */
REAL8 s2x, /**< initial value of S2x */
REAL8 s2y, /**< initial value of S2y */
REAL8 s2z /**< initial value of S2z */
);
int XLALSimInspiralNRSur7dq2Polarizations(
REAL8TimeSeries **hplus, /**< OUTPUT h_+ vector */
REAL8TimeSeries **hcross, /**< OUTPUT h_x vector */
REAL8 phiRef, /**< orbital phase at reference pt. */
REAL8 inclination, /**< inclination angle */
REAL8 deltaT, /**< sampling interval (s) */
REAL8 m1, /**< mass of companion 1 (kg) */
REAL8 m2, /**< mass of companion 2 (kg) */
REAL8 distnace, /**< distance of source (m) */
REAL8 fMin, /**< start GW frequency (Hz) */
REAL8 fRef, /**< reference GW frequency (Hz) */
REAL8 s1x, /**< reference value of S1x */
REAL8 s1y, /**< reference value of S1y */
REAL8 s1z, /**< reference value of S1z */
REAL8 s2x, /**< reference value of S2x */
REAL8 s2y, /**< reference value of S2y */
REAL8 s2z, /**< reference value of S2z */
bool quadrupole_only /**< If true, only use (2, +/-2) modes in the coorbital frame */
);
SphHarmTimeSeries *XLALSimInspiralNRSur7dq2Modes(
REAL8 phiRef, /**< orbital phase at reference pt. */
REAL8 deltaT, /**< sampling interval (s) */
REAL8 m1, /**< mass of companion 1 (kg) */
REAL8 m2, /**< mass of companion 2 (kg) */
REAL8 S1x, /**< x-component of the dimensionless spin of object 1 */
REAL8 S1y, /**< y-component of the dimensionless spin of object 1 */
REAL8 S1z, /**< z-component of the dimensionless spin of object 1 */
REAL8 S2x, /**< x-component of the dimensionless spin of object 2 */
REAL8 S2y, /**< y-component of the dimensionless spin of object 2 */
REAL8 S2z, /**< z-component of the dimensionless spin of object 2 */
REAL8 fMin, /**< start GW frequency (Hz) */
REAL8 fRef, /**< reference GW frequency (Hz) */
REAL8 distnace /**< distance of source (m) */
);
/* in module LALSimNRTunedTides.c */
double XLALSimNRTunedTidesComputeKappa2T(
REAL8 m1_SI, /**< Mass of companion 1 (kg) */
......
This diff is collapsed.
/* * Copyright (C) 2017 Jonathan Blackman
* NRSur7dq2 NR surrogate model.
* Paper: https://arxiv.org/abs/1705.07089
* Based on the python implementation found at:
* https://www.black-holes.org/data/surrogates/index.html
* which uses the same NRSur7dq2.hdf5 data file.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with with program; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
* MA 02111-1307 USA
*/
#include "LALSimNRSurrogateUtilities.c"
/****************************** Constants ***********************************/
// These are to rescale the mass ratio fit range from [0.99, 2.01] to [-1, 1].
static const double NRSUR7DQ2_Q_FIT_OFFSET = -2.9411764705882355; // = (0.5*(2.01 - 0.99)) * (2 / (2.01 - 0.99))
static const double NRSUR7DQ2_Q_FIT_SLOPE = 1.9607843137254901; // = 2 / (2.01 - 0.99)
// Allow a tiny bit of extrapolation in mass ratio and spin
static const double NRSUR7DQ2_Q_MIN = 0.999;
static const double NRSUR7DQ2_Q_MAX = 2.001;
static const double NRSUR7DQ2_CHI_MAX = 0.801;
static const int NRSUR7DQ2_LMAX = 4;
static const double NRSUR7DQ2_START_TIME = -4500.0;
// Surrogate model data, in LAL_DATA_PATH. File available in lalsuite-extra or at
// https://www.black-holes.org/surrogates
static const char NRSUR7DQ2_DATAFILE[] = "NRSur7dq2.h5";
/***********************************************************************************/
/****************************** Type definitions ***********************************/
/***********************************************************************************/
/**
* Data used in a single scalar NRSur7dq2 fit
*/
typedef struct tagFitData {
gsl_matrix_long *basisFunctionOrders; /**< matrix of (n_coefs x 7) basis function orders
giving the polynomial order in f(q), chiA components, and chiB components. */
gsl_vector *coefs; /**< coefficient vector of length n_coefs */
int n_coefs; /**< Number of coefficients in the fit */
} FitData;
/**
* Data used in a single vector NRSur7dq2 fit
*/
typedef struct tagVectorFitData {
gsl_matrix_long *basisFunctionOrders; /**< matrix of (n_coefs x 7) basis function orders
giving the polynomial order in f(q), chiA components, and chiB components. */
gsl_vector *coefs; /**< coefficient vector of length n_coefs */
gsl_vector_long *componentIndices; /**< Each fit coefficient applies to a single component
of the vector; this gives the component indices. */
int n_coefs; /**< Number of coefficients in the fit */
int vec_dim; /**< Dimension of the vector */
} VectorFitData;
/**
* NRSur7dq2 data for a single dynamics node
*/
typedef struct tagDynamicsNodeFitData {
FitData *omega_data; /**< A fit to the orbital angular frequency */
VectorFitData *omega_copr_data; /**< A 2d vector fit for the x and y components of
Omega^{coorb}(t) in arxiv 1705.07089 */
VectorFitData *chiA_dot_data; /**< A 3d vector fit for the coorbital components of the
time derivative of chiA taken in the coprecessing frame */
VectorFitData *chiB_dot_data; /**< A 3d vector fit for the coorbital components of the
time derivative of chiB taken in the coprecessing frame */
} DynamicsNodeFitData;
/**
* NRSur7dq2 data for a single waveform data piece.
* Waveform data pieces are real time-dependent components of the waveform in the coorbital frame.
* See the beginning of section IV of https://arxiv.org/abs/1705.07089
*/
typedef struct tagWaveformDataPiece {
int n_nodes; /**< Number of empirical nodes */
FitData **fit_data; /**< FitData at each empirical node */
gsl_matrix *empirical_interpolant_basis; /**< The empirical interpolation matrix */
gsl_vector_long *empirical_node_indices; /**< The empirical node indices */
} WaveformDataPiece;
/**
* All WaveformDataPieces needed to evaluate all modes with a fixed value of ell.
* For m=0 modes we model the real and imaginary parts separately.
* For m>0, we model the real and imaginary parts of
* X_pm^{ell, m} = frac{1}{2}( h^{ell, m} pm h^{ell, -m} )
* where this h is in the coorbital frame.
* (see https://arxiv.org/abs/1705.07089)
*/
typedef struct tagWaveformFixedEllModeData {
int ell; /**< The fixed value of ell */
WaveformDataPiece *m0_real_data; /**< The real (ell, 0) mode data piece */
WaveformDataPiece *m0_imag_data; /**< The imag (ell, 0) mode data piece */
WaveformDataPiece **X_real_plus_data; /**< One Re[X_+] for each 1 <= m <= ell */
WaveformDataPiece **X_real_minus_data; /**< One Re[X_-] for each 1 <= m <= ell */
WaveformDataPiece **X_imag_plus_data; /**< One Im[X_+] for each 1 <= m <= ell */
WaveformDataPiece **X_imag_minus_data; /**< One Im[X_-] for each 1 <= m <= ell */
} WaveformFixedEllModeData;
/**
* All data needed by the full NRSur7dq2 model
*/
typedef struct tagNRSur7dq2Data {
UINT4 setup; /**< Indicates if this has been initialized */
int LMax; /**< Maximum ell mode that will ever be evaluated */
gsl_vector *t_ds; /** Vector of the dynamics surrogate node times, not including half times */
gsl_vector *t_ds_half_times; /**< t_1/2, t_3/2, and t_5/2 used to start up integration*/
gsl_vector *t_coorb; /**< Vector of the coorbital surrogate output times. */
DynamicsNodeFitData **ds_node_data; /** A DynamicsNodeFitData for each time in t_ds.*/
DynamicsNodeFitData **ds_half_node_data; /** A DynamicsNodeFitData for each time in t_ds_half_times. */
WaveformFixedEllModeData **coorbital_mode_data; /** One for each 2 <= ell <= LMax */
} NRSur7dq2Data;
/***********************************************************************************/
/****************************** Function declarations*******************************/
/***********************************************************************************/
static void NRSur7dq2_Init_LALDATA(void);
static int NRSur7dq2_Init(NRSur7dq2Data *data, LALH5File *file);
static void NRSur7dq2_LoadDynamicsNode(DynamicsNodeFitData **ds_node_data, LALH5File *sub, int i);
static void NRSur7dq2_LoadCoorbitalEllModes(WaveformFixedEllModeData **coorbital_mode_data, LALH5File *file, int i);
static void NRSur7dq2_LoadWaveformDataPiece(LALH5File *sub, WaveformDataPiece **data, bool invert_sign);
static bool NRSur7dq2_IsSetup(void);
static double ipow(double base, int exponent); // integer powers
static double NRSur7dq2_eval_fit(FitData *data, double *x);
static void NRSur7dq2_eval_vector_fit(
double *res, // Result
VectorFitData *data, // Data for fit
double *x // size 7, giving mass ratio q, and dimensionless spin components
);
static void NRSur7dq2_normalize_y(
double chiANorm,
double chiBNorm,
double *y
); // Helper to keep spin magnitude constant
static void NRSur7dq2_normalize_results(
double normA,
double normB,
gsl_vector **quat,
gsl_vector **chiA,
gsl_vector **chiB
);
static void NRSur7dq2_rotate_spins(gsl_vector **chiA, gsl_vector **chiB, gsl_vector *phi);
static void NRSur7dq2_ds_fit_x(
double *x, //Result, length 7
double q, // Mass ratio
double *y // [q0, qx, qy, qz, phi, chiAx, chiAy, chiAz, chiBx, chiBy, chiBz]
);
static void NRSur7dq2_assemble_dydt(
double *dydt, // Result, length 11
double *y, // ODE solution at the current time step
double *Omega_coorb_xy, // a form of time derivative of the coprecessing frame
double omega, // orbital angular frequency in the coprecessing frame
double *chiA_dot, // chiA time derivative
double *chiB_dot // chiB time derivative
);
static double cubic_interp(double xout, double *x, double *y);
static gsl_vector *spline_array_interp(gsl_vector *xout, gsl_vector *x, gsl_vector *y);
static double NRSur7dq2_get_omega(size_t node_index, double q, double *y0);
static double NRSur7dq2_get_t_ref(double omega_ref, double q, double *chiA0, double *chiB0, double *q_ref, double phi_ref);
static void NRSur7dq2_get_time_deriv_from_index(
double *dydt, // Output: dy/dt evaluated at the ODE time node with index i0. Must have space for 11 entries.
int i0, // Time node index. i0=-1, -2, and -3 are used for time nodes 1/2, 3/2, and 5/2 respectively.
double q, // Mass ratio
double *y // Current ODE state: [q0, qx, qy, qz, orbphase, chiAx, chiAy, chiAz, chiBx, chiBy, chiBz]
);
static void NRSur7dq2_get_time_deriv(
double *dtdy, // Output: dy/dt evaluated at time t. Must have space for 11 entries.
double t, // Time at which the ODE should be evaluated.
double q, // Mass ratio
double *y // Current ODE state
);
static int NRSur7dq2_initialize_at_dynamics_node(
double *dynamics_data, // ODE output
double t_ref, // reference time. t_ds[i0] will be close to t_ref.
double q, // mass ratio
double *chiA0, // chiA at t_ref.
double *chiB0, // chiB at t_ref.
double phi_ref, // orbital phase at t_ref.
double *q_ref, // quaternion at t_ref.
double normA, // |chiA|
double normB // |chiB|
);
static void NRSur7dq2_initialize_RK4_with_half_nodes(
double *dynamics_data, // ODE output
double *time_steps, // Output: first three time steps. Should have size 3.
double *dydt0, // Output: dydt at node 0. Should have size 11.
double *dydt1, // Output: dydt at node 1. Should have size 11.
double *dydt2, // Output: dydt at node 2. Should have size 11.
double *dydt3, // Output: dydt at node 3. Should have size 11.
double normA, // |chiA|
double normB, // |chiB|
double q // mass ratio
);
static int NRSur7dq2_initialize_RK4(
double *dynamics_data, // ODE output
double *time_steps, // Output: first three time steps. Should have size 3.
double *dydt0, // Output: dydt at node i0 + 0. Should have size 11.
double *dydt1, // Output: dydt at node i0 + 1. Should have size 11.
double *dydt2, // Output: dydt at node i0 + 2. Should have size 11.
double *dydt3, // Output: dydt at node i0 + 3. Should have size 11.
double normA, // |chiA|
double normB, // |chiB|
double q, // mass ratio
int i0 // the node that is already initialized
);
static void NRSur7dq2_integrate_AB4(
double *dynamics_data, // ODE output
double *time_steps, // The first three time steps beginning at i_start.
double *dydt0, // dydt at node i_start
double *dydt1, // dydt at node i_start + 1
double *dydt2, // dydt at node i_start + 2
double *dydt3, // dydt at node i_start + 3
double normA, // |chiA|
double normB, // |chiB|
double q, // mass ratio
int i_start // nodes i_start through i_start+3 are already initialized
);
static int NRSur7dq2_IntegrateDynamics(
double *dynamics_data, // Output: length (n * 11), where entries 11*m <= i < 11*(m+1) are
// [q0, qx, qy, qz, varphi, chiAx, chiAy, chiAz, chiBx, chiBy, chiBz]
double q, // Mass ratio mA / mB
double *chiA0, // chiA at the reference point
double *chiB0, // chiB at the reference point
double omega_ref, // orbital angular frequency at the reference point
double phi_ref, // orbital phase at the reference point
double *q_ref // coprecessing quaterion at the reference point
);
static void NRSur7dq2_eval_data_piece(
gsl_vector *result, // Output: Should have already been assigned space
double q, // Mass ratio
gsl_vector **chiA, // 3 gsl_vector *s, one for each (coorbital) component
gsl_vector **chiB, // similar to chiA
WaveformDataPiece *data // The data piece to evaluate
);
static void NRSur7dq2_core(
MultiModalWaveform **h, // Output. Dimensionless waveform modes sampled on t_coorb
double q, // Mass ratio mA / mB
double *chiA0, // chiA at the reference point
double *chiB0, // chiB at the reference point
double omega_ref, // orbital angular frequency at the reference point
double phi_ref, // orbital phase at the reference point
double *q_ref, // coprecessing quaterion at the reference point
bool quadrupole_only // If true, only use (2, +/-2) modes in the coorbital frame
);
......@@ -52,6 +52,8 @@ UNUSED static int read_matrix(const char dir[], const char fname[], gsl_matrix *
UNUSED static int CheckVectorFromHDF5(LALH5File *file, const char name[], const double *v, size_t n);
UNUSED static int ReadHDF5RealVectorDataset(LALH5File *file, const char *name, gsl_vector **data);
UNUSED static int ReadHDF5RealMatrixDataset(LALH5File *file, const char *name, gsl_matrix **data);
UNUSED static int ReadHDF5LongVectorDataset(LALH5File *file, const char *name, gsl_vector_long **data);
UNUSED static int ReadHDF5LongMatrixDataset(LALH5File *file, const char *name, gsl_matrix_long **data);
UNUSED static void PrintInfoStringAttribute(LALH5File *file, const char attribute[]);
UNUSED static int ROM_check_version_number(LALH5File *file, INT4 version_major_in, INT4 version_minor_in, INT4 version_micro_in);
#endif
......@@ -105,6 +107,7 @@ static void err_handler(const char *reason, const char *file, int line, int gsl_
XLALPrintError("gsl: %s:%d: %s - %d\n", file, line, reason, gsl_errno);
}
// Helper functions to read gsl_vector and gsl_matrix data with error checking
static int read_vector(const char dir[], const char fname[], gsl_vector *v) {
size_t size = strlen(dir) + strlen(fname) + 2;
......@@ -262,6 +265,111 @@ static int ReadHDF5RealMatrixDataset(LALH5File *file, const char *name, gsl_matr
return 0;
}
static int ReadHDF5LongVectorDataset(LALH5File *file, const char *name, gsl_vector_long **data) {
LALH5Dataset *dset;
UINT4Vector *dimLength;
size_t n;
if (file == NULL || name == NULL || data == NULL)
XLAL_ERROR(XLAL_EFAULT);
dset = XLALH5DatasetRead(file, name);
if (dset == NULL)
XLAL_ERROR(XLAL_EFUNC);
if (XLALH5DatasetQueryType(dset) != LAL_I8_TYPE_CODE) {
XLALH5DatasetFree(dset);
XLAL_ERROR(XLAL_ETYPE, "Dataset `%s' is wrong type", name);
}
dimLength = XLALH5DatasetQueryDims(dset);
if (dimLength == NULL) {
XLALH5DatasetFree(dset);
XLAL_ERROR(XLAL_EFUNC);
}
if (dimLength->length != 1) {
XLALH5DatasetFree(dset);
XLAL_ERROR(XLAL_EDIMS, "Dataset `%s' must be 1-dimensional", name);
}
n = dimLength->data[0];
XLALDestroyUINT4Vector(dimLength);
if (*data == NULL) {
*data = gsl_vector_long_alloc(n);
if (*data == NULL) {
XLALH5DatasetFree(dset);
XLAL_ERROR(XLAL_ENOMEM, "gsl_vector_long_alloc(%zu) failed", n);
}
}
else if ((*data)->size != n) {
XLALH5DatasetFree(dset);
XLAL_ERROR(XLAL_EINVAL, "Expected gsl_vector `%s' of size %zu", name, n);
}
// Now read the data
if (XLALH5DatasetQueryData((*data)->data, dset) < 0) {
XLALH5DatasetFree(dset);
XLAL_ERROR(XLAL_EFUNC);
}
XLALH5DatasetFree(dset);
return 0;
}
static int ReadHDF5LongMatrixDataset(LALH5File *file, const char *name, gsl_matrix_long **data) {
LALH5Dataset *dset;
UINT4Vector *dimLength;
size_t n1, n2;
if (file == NULL || name == NULL || data == NULL)
XLAL_ERROR(XLAL_EFAULT);
dset = XLALH5DatasetRead(file, name);
if (dset == NULL)
XLAL_ERROR(XLAL_EFUNC);
if (XLALH5DatasetQueryType(dset) != LAL_I8_TYPE_CODE) {
XLALH5DatasetFree(dset);
XLAL_ERROR(XLAL_ETYPE, "Dataset `%s' is wrong type", name);
}
dimLength = XLALH5DatasetQueryDims(dset);
if (dimLength == NULL) {
XLALH5DatasetFree(dset);
XLAL_ERROR(XLAL_EFUNC);
}
if (dimLength->length != 2) {
XLALH5DatasetFree(dset);
XLAL_ERROR(XLAL_EDIMS, "Dataset `%s' must be 2-dimensional", name);
}
n1 = dimLength->data[0];
n2 = dimLength->data[1];
XLALDestroyUINT4Vector(dimLength);
if (*data == NULL) {
*data = gsl_matrix_long_alloc(n1, n2);
if (*data == NULL) {
XLALH5DatasetFree(dset);
XLAL_ERROR(XLAL_ENOMEM, "gsl_matrix_long_alloc(%zu, %zu) failed", n1, n2);
}
}
else if ((*data)->size1 != n1 || (*data)->size2 != n2) {
XLALH5DatasetFree(dset);
XLAL_ERROR(XLAL_EINVAL, "Expected gsl_matrix_long `%s' of size %zu x %zu", name, n1, n2);
}
// Now read the data
if (XLALH5DatasetQueryData((*data)->data, dset) < 0) {
XLALH5DatasetFree(dset);
XLAL_ERROR(XLAL_EFUNC);
}
XLALH5DatasetFree(dset);
return 0;
}
static void PrintInfoStringAttribute(LALH5File *file, const char attribute[]) {
LALH5Generic gfile = {.file = file};
int len = XLALH5AttributeQueryStringValue(NULL, 0, gfile, attribute) + 1;
......
......@@ -153,6 +153,8 @@ static const char *lalSimulationApproximantNames[] = {
INITIALIZE_NAME(SpinTaylorT4Fourier),
INITIALIZE_NAME(SpinTaylorT2Fourier),
INITIALIZE_NAME(SpinDominatedWf),
INITIALIZE_NAME(NRSur7dq2_HM),
INITIALIZE_NAME(NRSur7dq2_L2only),
INITIALIZE_NAME(NR_hdf5),
};
#undef INITIALIZE_NAME
......@@ -871,6 +873,20 @@ int XLALSimInspiralChooseTDWaveform(
S2x, S2y, S2z, numrel_data_path);
XLALFree(numrel_data_path);
break;
case NRSur7dq2_HM:
/* Waveform-specific sanity checks */
/* Call the waveform driver routine */
ret = XLALSimInspiralNRSur7dq2Polarizations(hplus, hcross,
phiRef, i, deltaT, m1, m2, r, f_min, f_ref,
S1x, S1y, S1z, S2x, S2y, S2z, false);
break;
case NRSur7dq2_L2only:
/* Waveform-specific sanity checks */
/* Call the waveform driver routine */
ret = XLALSimInspiralNRSur7dq2Polarizations(hplus, hcross,
phiRef, i, deltaT, m1, m2, r, f_min, f_ref,
S1x, S1y, S1z, S2x, S2y, S2z, true);
break;
default:
......@@ -4056,6 +4072,8 @@ int XLALSimInspiralImplementedTDApproximants(
case SEOBNRv4:
case SEOBNRv4_opt:
case NR_hdf5:
case NRSur7dq2_HM:
case NRSur7dq2_L2only:
return 1;
default:
......@@ -4488,6 +4506,8 @@ int XLALSimInspiralGetSpinSupportFromApproximant(Approximant approx){
case SpinDominatedWf:
case SEOBNRv3:
case NR_hdf5:
case NRSur7dq2_HM:
case NRSur7dq2_L2only:
spin_support=LAL_SIM_INSPIRAL_PRECESSINGSPIN;
break;
case SpinTaylorF2:
......@@ -4610,6 +4630,8 @@ int XLALSimInspiralApproximantAcceptTestGRParams(Approximant approx){
case SpinDominatedWf:
case NumApproximants:
case NR_hdf5:
case NRSur7dq2_HM:
case NRSur7dq2_L2only:
testGR_accept=LAL_SIM_INSPIRAL_NO_TESTGR_PARAMS;
break;
case TaylorF2:
......
......@@ -366,6 +366,8 @@ typedef enum tagApproximant {
SpinDominatedWf, /**< Time domain, inspiral only, 1 spin, precessing waveform, Tapai et al, arXiv: 1209.1722
* @remarks Implemented in lalsimulation (time domain). */
NR_hdf5, /**< Time domain, NR waveform from HDF file. From INSERT LINKS HERE */
NRSur7dq2_HM, /**< Time domain, fully precessing NR surrogate model with up to ell=4 modes, arxiv: 1705.07089 */
NRSur7dq2_L2only, /**< NRSur7dq2_HM but using only the (2, +/-2) modes in the coorbital frame */
SEOBNRv4_ROM_NRTidal, /**< Low-mass double-spin frequency domain reduced order model of spin-aligned EOBNR model SEOBNRv4 [Bohe et al, arXiv:1611.03703] with tidal phase corrections [Dietrich et al, arXiv:1706.02969]
* @remarks Implemented in lalsimulation (frequency domain). */
IMRPhenomD_NRTidal, /**< Uses arxiv:1706.02969 to upgrad IMRPhenomD to a tidal approximant
......
This diff is collapsed.
/* Copyright (C) 2018 Jonathan Blackman
* Utility functions that are useful for evaluating NR surrogates.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with with program; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
* MA 02111-1307 USA
*/
#include <math.h>
#include <lal/XLALError.h>
#include <lal/Units.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif
/**
* Structure for a multi-modal waveform. Can store any set of modes.
*/
typedef struct tagMultiModalWaveform {
REAL8 phiRef; /**< orbital phase at reference pt. */
int n_modes; /**< Number of modes */
int n_times; /**< Number of time samples in each mode */
int *lvals; /**< The ell values of each mode (n_modes entries) */
int *mvals; /**< The m values of each mode (n_modes entries) */
gsl_vector **modes_real_part; /**< The real part of each mode - n_modes vectors of length n_times */
gsl_vector **modes_imag_part; /**< The imaginary part of each mode - n_modes vectors of length n_times */
} MultiModalWaveform;
/**
* Helper structure for computing WignerDMatrices, which require x(t)^n for many values of n.
*/
typedef struct tagRealPowers {
int LMax; /**< Maximum L; powers computed depend on LMax. */
int n_entries; /**< Number of entries in **powers. */
int n_times; /**< Length of each entry. */
gsl_vector **powers; /**< The data; x^n. */
} RealPowers;
/**
* Helper structure for computing WignerDMatrices, which require z(t)^n for many values of n.
*/
typedef struct tagComplexPowers {
int LMax; /**< Maximum L; powers computed depend on LMax. */
int n_entries; /**< Number of entries in **powers. */
int n_times; /**< Length of each entry. */