Merge branch 'lalinference_o2_nrsur-dev-refactor' into 'lalinference_o2'

Refactor NRSur7dq2 and fix for waveform review

See merge request !343

lalsimulation/src/LALSimIMR.h

@@ -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) */

lalsimulation/src/LALSimIMRNRSur7dq2.h

+/* *  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
+);

lalsimulation/src/LALSimIMRSEOBNRROMUtilities.c

@@ -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;

lalsimulation/src/LALSimInspiral.c

@@ -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:

lalsimulation/src/LALSimInspiral.h

@@ -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

lalsimulation/src/LALSimNRSurrogateUtilities.h

+/*  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. */