25 #define UNUSED __attribute__ ((unused))
39 #include <gsl/gsl_errno.h>
40 #include <gsl/gsl_bspline.h>
41 #include <gsl/gsl_blas.h>
42 #include <gsl/gsl_min.h>
43 #include <gsl/gsl_spline.h>
44 #include <gsl/gsl_complex_math.h>
45 #include <lal/Units.h>
46 #include <lal/SeqFactories.h>
47 #include <lal/LALConstants.h>
48 #include <lal/XLALError.h>
49 #include <lal/FrequencySeries.h>
51 #include <lal/StringInput.h>
52 #include <lal/Sequence.h>
53 #include <lal/LALStdio.h>
54 #include <lal/FileIO.h>
55 #include <lal/SphericalHarmonics.h>
56 #include <lal/LALSimInspiral.h>
57 #include <lal/LALSimSphHarmMode.h>
58 #include <lal/LALSimIMR.h>
59 #include <lal/H5FileIO.h>
64 #include <lal/LALConfig.h>
65 #ifdef LAL_PTHREAD_LOCK
70 #ifdef LAL_PTHREAD_LOCK
71 static pthread_once_t NRSur7dq2_is_initialized = PTHREAD_ONCE_INIT;
74 #ifdef LAL_PTHREAD_LOCK
75 static pthread_once_t NRSur7dq4_is_initialized = PTHREAD_ONCE_INIT;
156 char *dir = dirname(
path);
184 char *dir = dirname(
path);
210 UINT4 PrecessingNRSurVersion
216 XLALPrintError(
"You tried to setup a NRSurrogate model that was already initialized. Ignoring.\n");
221 gsl_vector *t_ds_with_halves = NULL;
222 ReadHDF5RealVectorDataset(
file,
"t_ds", &t_ds_with_halves);
223 gsl_vector *t_ds = gsl_vector_alloc(t_ds_with_halves->size - 3);
224 gsl_vector *t_ds_half_times = gsl_vector_alloc(3);
225 for (
i=0;
i < 3;
i++) {
226 gsl_vector_set(t_ds,
i, gsl_vector_get(t_ds_with_halves, 2*
i));
227 gsl_vector_set(t_ds_half_times,
i, gsl_vector_get(t_ds_with_halves, 2*
i + 1));
229 for (
i=3;
i < t_ds->size;
i++) {
230 gsl_vector_set(t_ds,
i, gsl_vector_get(t_ds_with_halves,
i+3));
232 gsl_vector_free(t_ds_with_halves);
234 data->t_ds_half_times = t_ds_half_times;
243 for (
i=0;
i < (t_ds->size);
i++) ds_node_data[
i] = NULL;
244 for (
i=0;
i < 3;
i++) ds_half_node_data[
i] = NULL;
248 for (
i=0;
i < (t_ds->size);
i++) {
249 if (
i < 3) {j = 2*
i;}
else {j =
i+3;}
250 snprintf(sub_name, 20,
"ds_node_%d", j);
255 snprintf(sub_name, 20,
"ds_node_%d", j+1);
261 data->ds_node_data = ds_node_data;
262 data->ds_half_node_data = ds_half_node_data;
265 gsl_vector *t_coorb = NULL;
266 ReadHDF5RealVectorDataset(
file,
"t_coorb", &t_coorb);
267 data->t_coorb = t_coorb;
271 for (
int ell_idx=0; ell_idx <
NRSUR_LMAX-1; ell_idx++) {
274 data->coorbital_mode_data = coorbital_mode_data;
292 const size_t str_size = 30;
294 UNUSED
size_t nwritten;
296 nwritten = snprintf(tmp_name, str_size,
"%s_coefs",
name);
298 (*fit_data)->coefs = NULL;
299 ReadHDF5RealVectorDataset(sub, tmp_name, &((*fit_data)->coefs));
301 nwritten = snprintf(tmp_name, str_size,
"%s_bfOrders",
name);
303 (*fit_data)->basisFunctionOrders = NULL;
304 ReadHDF5LongMatrixDataset(sub, tmp_name,
305 &((*fit_data)->basisFunctionOrders));
307 (*fit_data)->n_coefs = (*fit_data)->coefs->size;
321 const size_t str_size = 20;
323 UNUSED
size_t nwritten;
326 (*vector_fit_data)->vec_dim = size;
329 for (
size_t i=0;
i<size;
i++) {
330 nwritten = snprintf(tmp_name, str_size,
"%s_%zu",
name,
i);
334 (*vector_fit_data)->fit_data[
i] = fit_data;
346 UINT4 PrecessingNRSurVersion
348 ds_node_data[
i] =
XLALMalloc(
sizeof(*ds_node_data[
i]) );
357 if (PrecessingNRSurVersion == 0){
360 omega_copr_data->
coefs = NULL;
363 ReadHDF5RealVectorDataset(sub,
"omega_orb_coefs", &(omega_copr_data->
coefs));
364 ReadHDF5LongMatrixDataset(sub,
"omega_orb_bfOrders", &(omega_copr_data->
basisFunctionOrders));
365 ReadHDF5LongVectorDataset(sub,
"omega_orb_bVecIndices", &(omega_copr_data->
componentIndices));
366 omega_copr_data->
n_coefs = omega_copr_data->
coefs->size;
372 chiA_dot_data->
coefs = NULL;
375 ReadHDF5RealVectorDataset(sub,
"chiA_coefs", &(chiA_dot_data->
coefs));
377 ReadHDF5LongVectorDataset(sub,
"chiA_bVecIndices", &(chiA_dot_data->
componentIndices));
385 chiB_dot_data->
coefs = NULL;
394 n = dimLength->
data[0];
398 ReadHDF5RealVectorDataset(sub,
"chiB_coefs", &(chiB_dot_data->
coefs));
400 ReadHDF5LongVectorDataset(sub,
"chiB_bVecIndices", &(chiB_dot_data->
componentIndices));
408 }
else if (PrecessingNRSurVersion == 1) {
437 mode_data->
ell =
i+2;
444 snprintf(sub_name, str_size,
"hCoorb_%d_0_real",
i+2);
449 snprintf(sub_name, str_size,
"hCoorb_%d_0_imag",
i+2);
471 for (
int m=1;
m<=(
i+2);
m++) {
472 snprintf(sub_name, str_size,
"hCoorb_%d_%d_Re+",
i+2,
m);
475 snprintf(sub_name, str_size,
"hCoorb_%d_%d_Re-",
i+2,
m);
478 snprintf(sub_name, str_size,
"hCoorb_%d_%d_Im+",
i+2,
m);
481 snprintf(sub_name, str_size,
"hCoorb_%d_%d_Im-",
i+2,
m);
486 coorbital_mode_data[
i] = mode_data;
500 gsl_matrix *EI_basis = NULL;
501 ReadHDF5RealMatrixDataset(sub,
"EIBasis", &EI_basis);
503 gsl_matrix_scale(EI_basis, -1);
505 (*data)->empirical_interpolant_basis = EI_basis;
507 gsl_vector_long *node_indices = NULL;
508 ReadHDF5LongVectorDataset(sub,
"nodeIndices", &node_indices);
509 (*data)->empirical_node_indices = node_indices;
511 int n_nodes = (*data)->empirical_node_indices->size;
512 (*data)->n_nodes = n_nodes;
518 for (
int i=0;
i<n_nodes;
i++) {
520 node_data->
coefs = NULL;
522 snprintf(sub_name, str_size,
"coefs_%d",
i);
523 ReadHDF5RealVectorDataset(nodeModelers, sub_name, &(node_data->
coefs));
524 snprintf(sub_name, str_size,
"bfOrders_%d",
i);
527 (*data)->fit_data[
i] = node_data;
552 if (exponent == 0)
return 1.0;
554 while (exponent > 1) {
584 for (
i=0;
i<22;
i++) {
586 x_powers[
i] =
ipow(q_fit,
i/7);
593 for (
i=0;
i <
data->n_coefs;
i++) {
595 prod = x_powers[7 * gsl_matrix_long_get(
data->basisFunctionOrders,
i, 0)];
597 for (j=1; j<7; j++) {
598 prod *= x_powers[7 * gsl_matrix_long_get(
data->basisFunctionOrders,
i, j) + j];
600 res += gsl_vector_get(
data->coefs,
i) * prod;
621 for (
i=0;
i <
data->vec_dim;
i++) {
629 for (
i=0;
i<22;
i++) {
631 x_powers[
i] =
ipow(q_fit,
i/7);
638 for (
i=0;
i <
data->n_coefs;
i++) {
640 prod = x_powers[7 * gsl_matrix_long_get(
data->basisFunctionOrders,
i, 0)];
642 for (j=1; j<7; j++) {
643 prod *= x_powers[7 * gsl_matrix_long_get(
data->basisFunctionOrders,
i, j) + j];
645 res[gsl_vector_long_get(
data->componentIndices,
i)] += gsl_vector_get(
data->coefs,
i) * prod;
662 const REAL8 chi_wtAvg = (
q*chi1z+chi2z)/(1+
q);
663 REAL8 chiHat_val = (chi_wtAvg - 38.*eta/113.*(chi1z
664 + chi2z))/(1. - 76.*eta/113.);
665 REAL8 chi_a_val = (chi1z - chi2z)/2.;
666 *chiHat = chiHat_val;
692 fit_params[0] = log(
x[0]);
693 fit_params[1] =
x[1];
694 fit_params[2] =
x[2];
695 fit_params[3] = chiHat;
696 fit_params[4] =
x[4];
697 fit_params[5] =
x[5];
698 fit_params[6] = chi_a;
705 for (
i=0;
i<22;
i++) {
707 x_powers[
i] =
ipow(q_fit,
i/7);
709 x_powers[
i] =
ipow(fit_params[
i%7],
i/7);
714 for (
i=0;
i <
data->n_coefs;
i++) {
716 prod = x_powers[7 * gsl_matrix_long_get(
data->basisFunctionOrders,
i, 0)];
718 for (j=1; j<7; j++) {
719 prod *= x_powers[7 * gsl_matrix_long_get(
data->basisFunctionOrders,
i, j) + j];
721 res += gsl_vector_get(
data->coefs,
i) * prod;
738 for (
int i=0;
i <
data->vec_dim;
i++) {
795 for (
i=0;
i<4;
i++) {
801 for (
i=5;
i<8;
i++) {
807 for (
i=8;
i<11;
i++) {
813 for (
i=0;
i<4;
i++) {
817 for (
i=5;
i<8;
i++) {
818 y[
i] =
y[
i] * chiANorm / nA;
820 for (
i=8;
i<11;
i++) {
821 y[
i] =
y[
i] * chiBNorm / nB;
837 int n = quat[0]->size;
838 REAL8 *chiAx = chiA[0]->data;
839 REAL8 *chiAy = chiA[1]->data;
840 REAL8 *chiAz = chiA[2]->data;
841 REAL8 *chiBx = chiB[0]->data;
842 REAL8 *chiBy = chiB[1]->data;
843 REAL8 *chiBz = chiB[2]->data;
844 REAL8 *q0 = quat[0]->data;
845 REAL8 *qx = quat[1]->data;
846 REAL8 *qy = quat[2]->data;
847 REAL8 *qz = quat[3]->data;
850 for (
i=0;
i<n;
i++) {
851 nA = sqrt(chiAx[
i]*chiAx[
i] + chiAy[
i]*chiAy[
i] + chiAz[
i]*chiAz[
i]);
852 chiAx[
i] *= normA / nA;
853 chiAy[
i] *= normA / nA;
854 chiAz[
i] *= normA / nA;
859 for (
i=0;
i<n;
i++) {
860 nB = sqrt(chiBx[
i]*chiBx[
i] + chiBy[
i]*chiBy[
i] + chiBz[
i]*chiBz[
i]);
861 chiBx[
i] *= normB / nB;
862 chiBy[
i] *= normB / nB;
863 chiBz[
i] *= normB / nB;
867 for (
i=0;
i<n;
i++) {
868 nQ = sqrt(q0[
i]*q0[
i] + qx[
i]*qx[
i] + qy[
i]*qy[
i] + qz[
i]*qz[
i]);
886 int n = phi_vec->size;
888 REAL8 *phi = phi_vec->data;
890 REAL8 *chix = chiA[0]->data;
891 REAL8 *chiy = chiA[1]->data;
892 for (
i=0;
i<n;
i++) {
896 chix[
i] = tmp*cp + chiy[
i]*sp;
897 chiy[
i] = -1*tmp*sp + chiy[
i]*cp;
900 chix = chiB[0]->data;
901 chiy = chiB[1]->data;
902 for (
i=0;
i<n;
i++) {
906 chix[
i] = tmp*cp + chiy[
i]*sp;
907 chiy[
i] = -1*tmp*sp + chiy[
i]*cp;
929 x[1] =
y[5]*cp +
y[6]*sp;
930 x[2] = -1*
y[5]*sp +
y[6]*cp;
934 x[4] =
y[8]*cp +
y[9]*sp;
935 x[5] = -1*
y[8]*sp +
y[9]*cp;
948 REAL8 *Omega_coorb_xy,
953 REAL8 omega_quat_x, omega_quat_y;
961 omega_quat_x = Omega_coorb_xy[0]*cp - Omega_coorb_xy[1]*sp;
962 omega_quat_y = Omega_coorb_xy[0]*sp + Omega_coorb_xy[1]*cp;
963 dydt[0] = (-0.5)*
y[1]*omega_quat_x - 0.5*
y[2]*omega_quat_y;
964 dydt[1] = (-0.5)*
y[3]*omega_quat_y + 0.5*
y[0]*omega_quat_x;
965 dydt[2] = 0.5*
y[3]*omega_quat_x + 0.5*
y[0]*omega_quat_y;
966 dydt[3] = 0.5*
y[1]*omega_quat_y - 0.5*
y[2]*omega_quat_x;
972 dydt[5] = chiA_dot[0]*cp - chiA_dot[1]*sp;
973 dydt[6] = chiA_dot[0]*sp + chiA_dot[1]*cp;
974 dydt[7] = chiA_dot[2];
975 dydt[8] = chiB_dot[0]*cp - chiB_dot[1]*sp;
976 dydt[9] = chiB_dot[0]*sp + chiB_dot[1]*cp;
977 dydt[10] = chiB_dot[2];
990 gsl_interp_accel *acc = gsl_interp_accel_alloc();
991 gsl_spline *interpolant = gsl_spline_alloc(gsl_interp_polynomial, 4);
992 gsl_spline_init(interpolant,
x,
y, 4);
993 REAL8 res = gsl_spline_eval(interpolant, xout, acc);
994 gsl_spline_free(interpolant);
995 gsl_interp_accel_free(acc);
1010 gsl_interp_accel *acc = gsl_interp_accel_alloc();
1011 gsl_spline *spline = gsl_spline_alloc(gsl_interp_cspline,
x->size);
1012 gsl_spline_init(spline,
x->data,
y->data,
x->size);
1014 gsl_vector *res = gsl_vector_alloc(xout->size);
1016 for (
size_t i=0;
i<xout->size;
i++) {
1017 tmp = gsl_spline_eval(spline, gsl_vector_get(xout,
i), acc);
1018 gsl_vector_set(res,
i, tmp);
1020 gsl_spline_free(spline);
1021 gsl_interp_accel_free(acc);
1050 REAL8 init_orbphase,
1054 REAL8 START_TIME = 100;
1061 if (fabs(omega_ref) < 1.e-10) {
1062 XLAL_PRINT_WARNING(
"WARNING: Treating omega_ref = 0 as a flag to use t_ref = t_0 = %.2f", START_TIME);
1068 "Reference frequency omega_ref=%0.4f > 0.2, too large!\n", omega_ref);
1072 for (
i=0;
i<4;
i++) y0[
i] = init_quat[
i];
1073 y0[4] = init_orbphase;
1074 for (
i=0;
i<3;
i++) {
1080 if (omega_ref < omega_min) {
1082 "Got omega_ref or omega_low=%0.4f smaller than the minimum omega=%0.4f for this configuration!", omega_ref, omega_min);
1088 gsl_vector *t_ds = __sur_data->
t_ds;
1089 while ((omega_max <= omega_ref) && (imax < t_ds->size)) {
1091 omega_min = omega_max;
1095 if (omega_max <= omega_ref) {
1100 REAL8 t_min = gsl_vector_get(t_ds, imax-1);
1101 REAL8 t_max = gsl_vector_get(t_ds, imax);
1102 REAL8 t_ref = (t_min * (omega_max - omega_ref) + t_max * (omega_ref - omega_min)) / (omega_max - omega_min);
1129 REAL8 omega, Omega_coorb_xy[2], chiA_dot[3], chiB_dot[3];
1152 gsl_vector *t_ds = __sur_data->
t_ds;
1154 int imax = t_ds->size - 1;
1155 REAL8 t1 = gsl_vector_get(t_ds, i1);
1156 REAL8 tmax = gsl_vector_get(t_ds, imax);
1157 if (t < t1 || t > tmax) {
1161 REAL8 t2 = gsl_vector_get(t_ds, i1+1);
1166 t2 = gsl_vector_get(t_ds, i1+1);
1172 if (i0 > imax-3) i0 = imax-3;
1173 REAL8 times[4], derivs[4], dydt0[11], dydt1[11], dydt2[11], dydt3[11];
1175 for (j=0; j<4; j++) times[j] = gsl_vector_get(t_ds, i0+j);
1181 for (j=0; j<11; j++) {
1182 derivs[0] = dydt0[j];
1183 derivs[1] = dydt1[j];
1184 derivs[2] = dydt2[j];
1185 derivs[3] = dydt3[j];
1197 REAL8 *dynamics_data,
1202 REAL8 init_orbphase,
1208 int imin, imax, i0, j;
1209 REAL8 tmin, tmax, t0,
dt, y_ref[11], dydt_ref[11], *node_data;
1213 t_ds = __sur_data->
t_ds;
1215 imax = t_ds->size - 1;
1216 tmin = gsl_vector_get(t_ds, imin);
1217 tmax = gsl_vector_get(t_ds, imax);
1218 while (imax - imin > 1) {
1219 i0 = (imax + imin)/2;
1220 t0 = gsl_vector_get(t_ds, i0);
1232 if (fabs(t_ref - tmin) < fabs(tmax - t_ref)) {
1241 for (j=0; j<4; j++) y_ref[j] = init_quat[j];
1242 y_ref[4] = init_orbphase;
1243 for (j=0; j<3; j++) {
1244 y_ref[5+j] = chiA0[j];
1245 y_ref[8+j] = chiB0[j];
1252 for (j=0; j<11; j++) y_ref[j] +=
dt * dydt_ref[j];
1256 node_data = dynamics_data + i0*11;
1257 for (j=0; j<11; j++) node_data[j] = y_ref[j];
1275 REAL8 *dynamics_data,
1287 gsl_vector *t_ds = __sur_data->
t_ds;
1299 t1 = gsl_vector_get(t_ds, 0);
1300 for (
i=0;
i<3;
i++) {
1301 t2 = gsl_vector_get(t_ds,
i+1);
1302 time_steps[
i] = t2 - t1;
1307 for (
i=0;
i<3;
i++ ) {
1310 node_data = dynamics_data + 11*
i;
1314 for (j=0; j<11; j++) {
1315 y_tmp[j] = node_data[j] + 0.5*time_steps[
i]*dydt[
i][j];
1320 for (j=0; j<11; j++) {
1321 y_tmp[j] = node_data[j] + 0.5*time_steps[
i]*
k2[j];
1326 for (j=0; j<11; j++) {
1327 y_tmp[j] = node_data[j] + time_steps[
i]*
k3[j];
1332 for (j=0; j<11; j++) {
1333 y_tmp[j] = node_data[j] + (time_steps[
i]/6.0)*(dydt[
i][j] + 2*
k2[j] + 2*
k3[j] +
k4[j]);
1336 node_data = node_data + 11;
1337 for (j=0; j<11; j++) {
1338 node_data[j] = y_tmp[j];
1343 node_data = dynamics_data + 33;
1354 REAL8 *dynamics_data,
1368 gsl_vector *t_ds = __sur_data->
t_ds;
1386 t1 = gsl_vector_get(t_ds, i_start);
1387 for (
i=0;
i<3;
i++) {
1388 t2 = gsl_vector_get(t_ds, i_start+
i+1);
1389 time_steps[
i] = t2 - t1;
1394 for (
i=0;
i<3;
i++ ) {
1397 node_data = dynamics_data + 11*(i0-
i);
1401 t1 = 0.5*(gsl_vector_get(t_ds, i0-
i) + gsl_vector_get(t_ds, i0-
i-1));
1402 for (j=0; j<11; j++) {
1403 y_tmp[j] = node_data[j] - 0.5*time_steps[2-
i]*dydt[3-
i][j];
1408 for (j=0; j<11; j++) {
1409 y_tmp[j] = node_data[j] - 0.5*time_steps[2-
i]*
k2[j];
1414 for (j=0; j<11; j++) {
1415 y_tmp[j] = node_data[j] - time_steps[2-
i]*
k3[j];
1420 for (j=0; j<11; j++) {
1421 y_tmp[j] = node_data[j] - (time_steps[2-
i]/6.0)*(dydt[3-
i][j] + 2*
k2[j] + 2*
k3[j] +
k4[j]);
1424 node_data = node_data - 11;
1425 for (j=0; j<11; j++) {
1426 node_data[j] = y_tmp[j];
1430 node_data = dynamics_data + (i0 - 3)*11;
1436 t1 = gsl_vector_get(t_ds, i0);
1437 for (
i=0;
i<3;
i++) {
1438 t2 = gsl_vector_get(t_ds, i0+
i+1);
1439 time_steps[
i] = t2 - t1;
1444 for (
i=0;
i<3;
i++ ) {
1447 node_data = dynamics_data + 11*(
i+i0);
1451 t1 = 0.5*(gsl_vector_get(t_ds,
i+i0) + gsl_vector_get(t_ds,
i+i0+1));
1452 for (j=0; j<11; j++) {
1453 y_tmp[j] = node_data[j] + 0.5*time_steps[
i]*dydt[
i][j];
1458 for (j=0; j<11; j++) {
1459 y_tmp[j] = node_data[j] + 0.5*time_steps[
i]*
k2[j];
1464 for (j=0; j<11; j++) {
1465 y_tmp[j] = node_data[j] + time_steps[
i]*
k3[j];
1470 for (j=0; j<11; j++) {
1471 y_tmp[j] = node_data[j] + (time_steps[
i]/6.0)*(dydt[
i][j] + 2*
k2[j] + 2*
k3[j] +
k4[j]);
1474 node_data = node_data + 11;
1475 for (j=0; j<11; j++) {
1476 node_data[j] = y_tmp[j];
1480 node_data = dynamics_data + (i0 + 3)*11;
1496 REAL8 *dynamics_data,
1509 REAL8 dt1, dt2, dt3, dt4;
1515 gsl_vector *t_ds = __sur_data->
t_ds;
1518 dt1 = time_steps[0];
1519 dt2 = time_steps[1];
1520 dt3 = time_steps[2];
1521 for (j=0; j<11; j++) {
1529 node_data = dynamics_data + 11*(i_start + 3);
1530 for (
i=i_start+4;
i< (int)t_ds->size;
i++) {
1533 dt4 = gsl_vector_get(t_ds,
i) - gsl_vector_get(t_ds,
i-1);
1537 for (j=0; j<11; j++) ynext[j] += node_data[j];
1542 for (j=0; j<11; j++) node_data[j] = ynext[j];
1545 if (
i < (
int) t_ds->size - 1) {
1549 for (j=0; j<11; j++) {
1559 dt1 = -1 * time_steps[2];
1560 dt2 = -1 * time_steps[1];
1561 dt3 = -1 * time_steps[0];
1562 for (j=0; j<11; j++) {
1570 node_data = dynamics_data + 11*(i_start);
1571 for (
i=i_start-1;
i>=0;
i--) {
1574 dt4 = gsl_vector_get(t_ds,
i) - gsl_vector_get(t_ds,
i+1);
1578 for (j=0; j<11; j++) ynext[j] += node_data[j];
1583 for (j=0; j<11; j++) node_data[j] = ynext[j];
1590 for (j=0; j<11; j++) {
1614 gsl_vector *nodes = gsl_vector_alloc(
data->n_nodes);
1616 int i, j, node_index;
1620 for (
i=0;
i<
data->n_nodes;
i++) {
1621 node_index = gsl_vector_long_get(
data->empirical_node_indices,
i);
1622 for (j=0; j<3; j++) {
1623 x[1+j] = gsl_vector_get(chiA[j], node_index);
1624 x[4+j] = gsl_vector_get(chiB[j], node_index);
1630 gsl_blas_dgemv(CblasTrans, 1.0,
data->empirical_interpolant_basis, nodes, 0.0, result);
1632 gsl_vector_free(nodes);
1651 REAL8 *dynamics_data,
1656 REAL8 init_orbphase,
1659 UINT4 PrecessingNRSurVersion
1663 UINT4 unlim_extrap = 0;
1664 if (LALparams != NULL &&
1671 REAL8 normA = sqrt(chiA0[0]*chiA0[0] + chiA0[1]*chiA0[1] + chiA0[2]*chiA0[2]);
1672 REAL8 normB = sqrt(chiB0[0]*chiB0[0] + chiB0[1]*chiB0[1] + chiB0[2]*chiB0[2]);
1677 REAL8 Q_MAX_WARN = 1;
1678 REAL8 CHI_MAX_WARN = 0;
1680 if (PrecessingNRSurVersion == 0) {
1686 }
else if (PrecessingNRSurVersion == 1) {
1694 "Only 0 or 1 are currently allowed for PrecessingNRSurVersion\n");
1699 "Invalid mass ratio q = %0.4f < 1\n",
q);
1701 if ((
q >
Q_MAX) && (unlim_extrap == 0)) {
1703 "Too much extrapolation. Mass ratio q = %0.4f > %0.4f, the "
1704 "maximum allowed value.\n",
q,
Q_MAX);
1706 if (
q > Q_MAX_WARN) {
1708 "Extrapolating to mass ratio q = %0.4f > %0.4f, the maximum "
1709 "mass ratio used to train the surrogate.\n",
q, Q_MAX_WARN);
1711 if ((normA > CHI_MAX) && (unlim_extrap == 0)) {
1713 "Too much extrapolation. Spin magnitude |chiA| = %0.4f > %.04f "
1714 "the maximum allowed .\n", normA, CHI_MAX);
1716 if ((normB > CHI_MAX) && (unlim_extrap == 0)) {
1718 "Too much extrapolation. Spin magnitude |chiB| = %0.4f > %.04f "
1719 "the maximum allowed value.\n", normB, CHI_MAX);
1721 if (normA > CHI_MAX_WARN) {
1723 "Extrapolating to spin magnitude |chiA| = %0.4f > %0.4f, the "
1724 "maximum spin magnitude used to train the surrogate.\n",
1725 normA, CHI_MAX_WARN);
1727 if (normB > CHI_MAX_WARN) {
1729 "Extrapolating to spin magnitude |chiB| = %0.4f > %0.4f, the "
1730 "maximum spin magnitude used to train the surrogate.\n",
1731 normB, CHI_MAX_WARN);
1735 init_quat, init_orbphase, __sur_data);
1748 REAL8 time_steps[3];
1749 REAL8 dydt0[11], dydt1[11], dydt2[11], dydt3[11];
1755 i_start =
PrecessingNRSur_initialize_RK4(dynamics_data, time_steps, dydt0, dydt1, dydt2, dydt3, normA, normB,
q, i0, __sur_data);
1760 PrecessingNRSur_integrate_AB4(dynamics_data, time_steps, dydt0, dydt1, dydt2, dydt3, normA, normB,
q, i_start, __sur_data);
1780 #ifdef LAL_PTHREAD_LOCK
1792 #ifdef LAL_PTHREAD_LOCK
1805 " and NRSur7dq4 are allowed.\n");
1825 REAL8 init_orbphase,
1827 LALValue* ModeArray,
1834 if (!__sur_data->
setup) {
1841 for (
m=-ell;
m<=ell;
m++) {
1849 gsl_vector *t_ds = __sur_data->
t_ds;
1850 gsl_vector *t_coorb = __sur_data->
t_coorb;
1851 int n_ds = t_ds->size;
1852 int n_coorb = t_coorb->size;
1858 omega_ref, init_orbphase, init_quat, LALparams,
1867 gsl_vector *dynamics[11];
1868 for (j=0; j<11; j++) {
1869 dynamics[j] = gsl_vector_alloc(n_ds);
1871 for (
i=0;
i<n_ds;
i++) {
1872 for (j=0; j<11; j++) {
1873 gsl_vector_set(dynamics[j],
i, dynamics_data[11*
i + j]);
1883 gsl_vector *quat_coorb[4], *chiA_coorb[3], *chiB_coorb[3];
1896 for (j=0; j<11; j++) {
1897 gsl_vector_free(dynamics[j]);
1901 REAL8 normA = sqrt(chiA0[0]*chiA0[0] + chiA0[1]*chiA0[1] + chiA0[2]*chiA0[2]);
1902 REAL8 normB = sqrt(chiB0[0]*chiB0[0] + chiB0[1]*chiB0[1] + chiB0[2]*chiB0[2]);
1911 gsl_vector *data_piece_eval = gsl_vector_alloc(n_coorb);
1917 i0 = ell*(ell+1) - 4;
1933 for (
m=1;
m<=ell;
m++) {
1974 for (
i=0;
i<3;
i++) {
1975 gsl_vector_free(chiA_coorb[
i]);
1976 gsl_vector_free(chiB_coorb[
i]);
1977 gsl_vector_free(quat_coorb[
i]);
1979 gsl_vector_free(quat_coorb[3]);
1980 gsl_vector_free(phi_coorb);
1981 gsl_vector_free(data_piece_eval);
2010 for (
i=1;
i<4;
i++) y0[
i] = 0.0;
2020 REAL8 f_start_Hz = omega_dimless_start / (
LAL_PI * Mtot_sec);
2150 if ( ModeArray == NULL ) {
2160 &s1x, &s1y, &s1z, &s2x, &s2y, &s2z);
2165 REAL8 Mtot = massA+massB;
2168 REAL8 chiA0[3], chiB0[3];
2176 REAL8 omegaMin_dimless;
2177 REAL8 omegaRef_dimless;
2179 fMin, fRef, Mtot_sec);
2196 REAL8 init_orbphase = 0;
2201 q, chiA0, chiB0, omegaRef_dimless, init_orbphase, init_quat,
2204 if (!h_inertial_modes) {
2210 size_t length = h_inertial_modes->
n_times;
2211 gsl_vector *hplus_model_times = gsl_vector_calloc(length);
2212 gsl_vector *hcross_model_times = gsl_vector_calloc(length);
2219 bool skip_ell_modes;
2221 REAL8 prefactor = 1;
2226 skip_ell_modes =
true;
2227 for (
m=-ell;
m<=ell;
m++) {
2229 skip_ell_modes =
false;
2232 if (skip_ell_modes) {
2236 for (
m=-ell;
m<=ell;
m++) {
2248 0.5 *
LAL_PI - phiRef, -2, ell,
m);
2249 c_re = creal(swsh_coef);
2250 c_im = cimag(swsh_coef);
2258 if ((
m%2 != 0) && (labels_switched)) {
2263 for (j=0; j < (int)length; j++) {
2266 hplus_model_times->data[j] += prefactor
2267 * (c_re * hmre - c_im * hmim);
2268 hcross_model_times->data[j] -= prefactor
2269 * (c_im * hmre + c_re * hmim);
2277 gsl_vector_scale(hplus_model_times, a0);
2278 gsl_vector_scale(hcross_model_times, a0);
2281 gsl_vector *model_times = __sur_data->
t_coorb;
2283 REAL8 t0 = gsl_vector_get(model_times, 0);
2285 if (fMin >= start_freq) {
2287 init_quat, init_orbphase, __sur_data);
2288 }
else if (fMin > 0) {
2291 gsl_vector_free(hplus_model_times);
2292 gsl_vector_free(hcross_model_times);
2296 REAL8 tf = gsl_vector_get(model_times, length-1);
2297 int nt = (int) ceil((tf - t0) / deltaT_dimless);
2298 gsl_vector *output_times = gsl_vector_alloc(nt);
2299 for (j=0; j<nt; j++) {
2300 gsl_vector_set(output_times, j, t0 + deltaT_dimless*j);
2309 gsl_interp_accel *acc = gsl_interp_accel_alloc();
2310 gsl_spline *spl_hplus = gsl_spline_alloc(gsl_interp_cspline, length);
2311 gsl_spline *spl_hcross = gsl_spline_alloc(gsl_interp_cspline, length);
2312 gsl_spline_init(spl_hplus, model_times->data, hplus_model_times->data, length);
2313 gsl_spline_init(spl_hcross, model_times->data, hcross_model_times->data, length);
2314 for (j=0; j<nt; j++) {
2315 t = gsl_vector_get(output_times, j);
2316 (*hplus)->data->data[j] = gsl_spline_eval(spl_hplus, t, acc);
2317 (*hcross)->data->data[j] = gsl_spline_eval(spl_hcross, t, acc);
2321 gsl_vector_free(hplus_model_times);
2322 gsl_vector_free(hcross_model_times);
2323 gsl_vector_free(output_times);
2324 gsl_interp_accel_free(acc);
2325 gsl_spline_free(spl_hplus);
2326 gsl_spline_free(spl_hcross);
2416 if ( ModeArray == NULL ) {
2426 &S1x, &S1y, &S1z, &S2x, &S2y, &S2z);
2431 REAL8 Mtot = massA+massB;
2434 REAL8 chiA0[3], chiB0[3];
2442 REAL8 omegaMin_dimless;
2443 REAL8 omegaRef_dimless;
2445 fMin, fRef, Mtot_sec);
2462 REAL8 init_orbphase = 0;
2467 chiA0, chiB0, omegaRef_dimless, init_orbphase, init_quat,
2479 gsl_vector *model_times = __sur_data->
t_coorb;
2481 REAL8 t0 = gsl_vector_get(model_times, 0);
2483 if (fMin >= start_freq) {
2485 init_quat, init_orbphase, __sur_data);
2486 }
else if (fMin > 0) {
2492 size_t length = model_times->size;
2493 REAL8 tf = gsl_vector_get(model_times, length-1);
2494 int nt = (int) ceil((tf - t0) / deltaT_dimless);
2495 gsl_vector *output_times = gsl_vector_alloc(nt);
2497 for (j=0; j<nt; j++) {
2498 gsl_vector_set(output_times, j, t0 + deltaT_dimless*j);
2505 gsl_interp_accel *acc = gsl_interp_accel_alloc();
2510 bool skip_ell_modes;
2512 REAL8 prefactor = 1;
2516 skip_ell_modes =
true;
2517 for (
m=-ell;
m<=ell;
m++) {
2520 if (skip_ell_modes) {
2524 for (
m=-ell;
m<=ell;
m++) {
2532 if ((
m%2 != 0) && (labels_switched)) {
2540 snprintf(mode_name,
sizeof(mode_name),
"(%d, %d) mode", ell,
m);
2543 gsl_spline *spl_re = gsl_spline_alloc(gsl_interp_cspline, length);
2544 gsl_spline *spl_im = gsl_spline_alloc(gsl_interp_cspline, length);
2545 gsl_spline_init(spl_re, model_times->data,
2547 gsl_spline_init(spl_im, model_times->data,
2549 for (j=0; j<nt; j++) {
2550 t = gsl_vector_get(output_times, j);
2551 tmp_mode->
data->
data[j] = gsl_spline_eval(spl_re, t, acc);
2552 tmp_mode->
data->
data[j] += I * gsl_spline_eval(spl_im, t, acc);
2554 gsl_spline_free(spl_re);
2555 gsl_spline_free(spl_im);
2564 gsl_vector_free(output_times);
2565 gsl_interp_accel_free(acc);
2668 gsl_vector **t_dynamics,
2673 gsl_vector **orbphase,
2687 REAL8 omegaRef_dimless,
2692 REAL8 init_orbphase,
2699 if (!__sur_data->
setup) {
2707 REAL8 chiA0[3], chiB0[3];
2708 chiA0[0] = chiA0x * cos(init_orbphase) - chiA0y * sin(init_orbphase);
2709 chiA0[1] = chiA0y * cos(init_orbphase) + chiA0x * sin(init_orbphase);
2711 chiB0[0] = chiB0x * cos(init_orbphase) - chiB0y * sin(init_orbphase);
2712 chiB0[1] = chiB0y * cos(init_orbphase) + chiB0x * sin(init_orbphase);
2717 init_quat[0] = init_quat0;
2718 init_quat[1] = init_quat1;
2719 init_quat[2] = init_quat2;
2720 init_quat[3] = init_quat3;
2724 int n_ds = (__sur_data->
t_ds)->size;
2727 chiA0, chiB0, omegaRef_dimless, init_orbphase, init_quat,
2736 gsl_vector *dynamics[11];
2737 for (j=0; j<11; j++) {
2738 dynamics[j] = gsl_vector_alloc(n_ds);
2740 for (
i=0;
i<n_ds;
i++) {
2741 for (j=0; j<11; j++) {
2742 gsl_vector_set(dynamics[j],
i, dynamics_data[11*
i + j]);
2749 *t_dynamics = gsl_vector_alloc(n_ds);
2750 gsl_vector_memcpy(*t_dynamics, __sur_data->
t_ds);
2752 *quat0 = dynamics[0];
2753 *quat1 = dynamics[1];
2754 *quat2 = dynamics[2];
2755 *quat3 = dynamics[3];
2756 *orbphase = dynamics[4];
2757 *chiAx = dynamics[5];
2758 *chiAy = dynamics[6];
2759 *chiAz = dynamics[7];
2760 *chiBx = dynamics[8];
2761 *chiBy = dynamics[9];
2762 *chiBz = dynamics[10];
struct tagLALH5File LALH5File
struct tagLALH5Dataset LALH5Dataset
int XLALDictContains(const LALDict *dict, const char *key)
UINT4 XLALDictLookupUINT4Value(LALDict *dict, const char *key)
static const double Q_MAX
static PrecessingNRSurData __lalsim_NRSur7dq2_data
Global surrogate data.
static PrecessingNRSurData __lalsim_NRSur7dq4_data
static const double NRSUR7DQ4_Q_MAX
static const double NRSUR7DQ4_Q_FIT_OFFSET
static const char NRSUR7DQ4_DATAFILE[]
static const double NRSUR7DQ4_START_TIME
static const double NRSUR7DQ4_Q_FIT_SLOPE
static const double NRSUR7DQ2_Q_FIT_SLOPE
static const double NRSUR7DQ2_Q_MAX
static const double NRSUR7DQ4_Q_MAX_WARN
static const double NRSUR7DQ2_CHI_MAX_WARN
static const int NRSUR_LMAX
static const double NRSUR7DQ2_CHI_MAX
static const double NRSUR7DQ2_Q_MAX_WARN
static const double NRSUR7DQ4_CHI_MAX_WARN
static const double NRSUR7DQ2_Q_FIT_OFFSET
static const double NRSUR7DQ4_CHI_MAX
static const double NRSUR7DQ2_START_TIME
static const char NRSUR7DQ2_DATAFILE[]
Auxiliary functions for SEOBNRv1/v2 reduced order modeling codes in LALSimIMRSEOBNRv2ROMDoubleSpin....
static void NRSur_ab4_dy(REAL8 *dy, REAL8 *k1, REAL8 *k2, REAL8 *k3, REAL8 *k4, REAL8 t12, REAL8 t23, REAL8 t34, REAL8 t45, int dim)
static void MultiModalWaveform_Destroy(MultiModalWaveform *wave)
Destroy a MultiModalWaveform structure; free all associated memory.
static void TransformModesCoorbitalToInertial(MultiModalWaveform *h, MultiModalWaveform *h_coorb, gsl_vector **quat, gsl_vector *orbphase)
Transforms modes from the coorbital frame to the inertial frame.
static void MultiModalWaveform_Init(MultiModalWaveform **wave, int LMax, int n_times)
Initialize a MultiModalWaveforme.
static UNUSED int get_dimless_omega(REAL8 *omegaMin_dimless, REAL8 *omegaRef_dimless, const REAL8 fMin, const REAL8 fRef, const REAL8 Mtot_sec)
Wrapper to get dimensionless omegaMin/omegaRef from fMin/fRef.
void XLALDestroyValue(LALValue *value)
static REAL8 sum(REAL8 *arr1, REAL8 *arr2, int k)
#define XLAL_FILE_RESOLVE_PATH(fname)
UINT4Vector * XLALH5DatasetQueryDims(LALH5Dataset UNUSED *dset)
LALH5File * XLALH5GroupOpen(LALH5File UNUSED *file, const char UNUSED *name)
LALH5File * XLALH5FileOpen(const char UNUSED *path, const char UNUSED *mode)
LALH5Dataset * XLALH5DatasetRead(LALH5File UNUSED *file, const char UNUSED *name)
void * XLALCalloc(size_t m, size_t n)
void * XLALMalloc(size_t n)
#define LAL_SIM_L_MAX_MODE_ARRAY
Maximum L spherical harmonic mode that is supported in Mode Array.
Approximant
Enum that specifies the PN approximant to be used in computing the waveform.
@ NRSur7dq4
q=4 extension of NRSur7dq2, arxiv: 1905.09300
@ NRSur7dq2
Time domain, fully precessing NR surrogate model with up to ell=4 modes, arxiv: 1705....
static PrecessingNRSurData * PrecessingNRSur_LoadData(Approximant approximant)
Wrapper for Loading NRSur7dq2 and NRSur7dq4 data.
static REAL8 PrecessingNRSur_get_omega(size_t node_index, REAL8 q, REAL8 *y0, PrecessingNRSurData *__sur_data)
Computes the orbital angular frequency at a dynamics node.
static bool NRSur7dq4_IsSetup(void)
Helper function which returns whether or not the global NRSur7dq4 surrogate data has been initialized...
static PrecessingNRSurData * PrecessingNRSur_core(MultiModalWaveform **h, REAL8 q, REAL8 *chiA0, REAL8 *chiB0, REAL8 omega_ref, REAL8 init_orbphase, REAL8 *init_quat, LALValue *ModeArray, LALDict *LALparams, Approximant approximant)
static void PrecessingNRSur_normalize_y(REAL8 chiANorm, REAL8 chiBNorm, REAL8 *y)
static int PrecessingNRSur_initialize_RK4(REAL8 *dynamics_data, REAL8 *time_steps, REAL8 *dydt0, REAL8 *dydt1, REAL8 *dydt2, REAL8 *dydt3, REAL8 normA, REAL8 normB, REAL8 q, int i0, PrecessingNRSurData *__sur_data)
Initializes the AB4 ODE system from a single arbitrary dynamics node by taking 3 RK4 steps.
static int PrecessingNRSur_Init(PrecessingNRSurData *data, LALH5File *file, UINT4 PrecessingNRSurVersion)
Initialize a PrecessingNRSurData structure from an open hdf5 file.
static void NRSur7dq2_Init_LALDATA(void)
This needs to be called once, before __lalsim_NRSur7dq2_data is used.
static void PrecessingNRSur_rotate_spins(gsl_vector **chiA, gsl_vector **chiB, gsl_vector *phi_vec)
Transforms chiA and chiB from the coprecessing frame to the coorbital frame using the orbital phase.
SphHarmTimeSeries * XLALSimInspiralPrecessingNRSurModes(REAL8 deltaT, REAL8 m1, REAL8 m2, REAL8 S1x, REAL8 S1y, REAL8 S1z, REAL8 S2x, REAL8 S2y, REAL8 S2z, REAL8 fMin, REAL8 fRef, REAL8 distance, LALDict *LALparams, Approximant approximant)
This function evaluates the NRSur7dq2 or NRSur7dq4 surrogate model and returns the inertial frame mod...
static void PrecessingNRSur_LoadWaveformDataPiece(LALH5File *sub, WaveformDataPiece **data, bool invert_sign)
Loads a single NRSur coorbital waveform data piece from file into a WaveformDataPiece.
static void PrecessingNRSur_get_time_deriv_from_index(REAL8 *dydt, int i0, REAL8 q, REAL8 *y, PrecessingNRSurData *__sur_data)
Compute dydt at a given dynamics node, where y is the numerical solution to the dynamics ODE.
int XLALSimInspiralPrecessingNRSurPolarizations(REAL8TimeSeries **hplus, REAL8TimeSeries **hcross, REAL8 phiRef, REAL8 inclination, REAL8 deltaT, REAL8 m1, REAL8 m2, REAL8 distance, REAL8 fMin, REAL8 fRef, REAL8 s1x, REAL8 s1y, REAL8 s1z, REAL8 s2x, REAL8 s2y, REAL8 s2z, LALDict *LALparams, Approximant approximant)
This function evaluates the NRSur7dq2 or NRSur7dq4 surrogate model and sums over all ell <= 4 modes t...
static void PrecessingNRSur_get_time_deriv(REAL8 *dydt, REAL8 t, REAL8 q, REAL8 *y, PrecessingNRSurData *__sur_data)
Compute dydt at any time by evaluating dydt at 4 nearby dynamics nodes and using cubic spline interpo...
static bool NRSur7dq2_IsSetup(void)
Helper function which returns whether or not the global NRSur7dq2 surrogate data has been initialized...
static void PrecessingNRSur_assemble_dydt(REAL8 *dydt, REAL8 *y, REAL8 *Omega_coorb_xy, REAL8 omega, REAL8 *chiA_dot, REAL8 *chiB_dot)
static void PrecessingNRSur_LoadCoorbitalEllModes(WaveformFixedEllModeData **coorbital_mode_data, LALH5File *file, int i)
Load the WaveformFixedEllModeData from file for a single value of ell.
static void PrecessingNRSur_integrate_AB4(REAL8 *dynamics_data, REAL8 *time_steps, REAL8 *dydt0, REAL8 *dydt1, REAL8 *dydt2, REAL8 *dydt3, REAL8 normA, REAL8 normB, REAL8 q, int i_start, PrecessingNRSurData *__sur_data)
Integrates the AB4 ODE system in time forwards, and backwards if needed.
static void NRSur7dq4_Init_LALDATA(void)
This needs to be called once, before __lalsim_NRSur7dq4_data is used.
static int NRSur7dq4_effective_spins(REAL8 *chiHat, REAL8 *chi_a, const REAL8 q, const REAL8 chi1z, const REAL8 chi2z)
int XLALPrecessingNRSurDynamics(gsl_vector **t_dynamics, gsl_vector **quat0, gsl_vector **quat1, gsl_vector **quat2, gsl_vector **quat3, gsl_vector **orbphase, gsl_vector **chiAx, gsl_vector **chiAy, gsl_vector **chiAz, gsl_vector **chiBx, gsl_vector **chiBy, gsl_vector **chiBz, REAL8 q, REAL8 chiA0x, REAL8 chiA0y, REAL8 chiA0z, REAL8 chiB0x, REAL8 chiB0y, REAL8 chiB0z, REAL8 omegaRef_dimless, REAL8 init_quat0, REAL8 init_quat1, REAL8 init_quat2, REAL8 init_quat3, REAL8 init_orbphase, LALDict *LALparams, Approximant approximant)
This function evaluates the NRSur7dq2 or NRSur7dq4 surrogate model and returns the precessing frame d...
static int PrecessingNRSur_initialize_at_dynamics_node(REAL8 *dynamics_data, REAL8 t_ref, REAL8 q, REAL8 *chiA0, REAL8 *chiB0, REAL8 init_orbphase, REAL8 *init_quat, REAL8 normA, REAL8 normB, PrecessingNRSurData *__sur_data)
Initialize the dynamics ODE at a dynamics node.
static void PrecessingNRSur_LoadFitData(FitData **fit_data, LALH5File *sub, const char *name)
Loads a single fit for NRSur7dq2 or NRSur7dq4.
static void PrecessingNRSur_normalize_results(REAL8 normA, REAL8 normB, gsl_vector **quat, gsl_vector **chiA, gsl_vector **chiB)
static int PrecessingNRSur_IntegrateDynamics(REAL8 *dynamics_data, REAL8 q, REAL8 *chiA0, REAL8 *chiB0, REAL8 omega_ref, REAL8 init_orbphase, REAL8 *init_quat, LALDict *LALparams, UINT4 PrecessingNRSurVersion)
This is the main NRSur dynamics surrogate integration routine.
static void NRSur7dq4_eval_vector_fit(REAL8 *res, VectorFitData *data, REAL8 *x)
REAL8 NRSur7dq4_eval_fit(FitData *data, REAL8 *x)
static void NRSur7dq2_eval_vector_fit(REAL8 *res, VectorFitData *data, REAL8 *x)
static bool PrecessingNRSur_switch_labels_if_needed(REAL8 *m1, REAL8 *m2, REAL8 *s1x, REAL8 *s1y, REAL8 *s1z, REAL8 *s2x, REAL8 *s2y, REAL8 *s2z)
If m1<m2, swaps the labels for the two BHs such that Bh1 is always heavier.
static void PrecessingNRSur_LoadDynamicsNode(DynamicsNodeFitData **ds_node_data, LALH5File *sub, int i, UINT4 PrecessingNRSurVersion)
Loads the data for a single dynamics node into a DynamicsNodeFitData struct.
REAL8 ipow(REAL8 base, int exponent)
Helper function for integer powers.
static void PrecessingNRSur_eval_vector_fit(REAL8 *res, VectorFitData *data, REAL8 *x, PrecessingNRSurData *__sur_data)
static gsl_vector * spline_array_interp(gsl_vector *xout, gsl_vector *x, gsl_vector *y)
Do cubic spline interpolation using a gsl_interp_cspline.
static REAL8 PrecessingNRSur_get_t_ref(REAL8 omega_ref, REAL8 q, REAL8 *chiA0, REAL8 *chiB0, REAL8 *init_quat, REAL8 init_orbphase, PrecessingNRSurData *__sur_data)
Computes a reference time from a reference orbital angular frequency.
REAL8 PrecessingNRSur_eval_fit(FitData *data, REAL8 *x, PrecessingNRSurData *__sur_data)
static void PrecessingNRSur_initialize_RK4_with_half_nodes(REAL8 *dynamics_data, REAL8 *time_steps, REAL8 *dydt0, REAL8 *dydt1, REAL8 *dydt2, REAL8 *dydt3, REAL8 normA, REAL8 normB, REAL8 q, PrecessingNRSurData *__sur_data)
Initializes the AB4 ODE system from the surrogate start time by taking 3 RK4 steps,...
static void PrecessingNRSur_eval_data_piece(gsl_vector *result, REAL8 q, gsl_vector **chiA, gsl_vector **chiB, WaveformDataPiece *data, PrecessingNRSurData *__sur_data)
Evaluates a single NRSur coorbital waveoform data piece.
static void PrecessingNRSur_ds_fit_x(REAL8 *x, REAL8 q, REAL8 *y)
REAL8 PrecessingNRSur_StartFrequency(REAL8 m1, REAL8 m2, REAL8 s1x, REAL8 s1y, REAL8 s1z, REAL8 s2x, REAL8 s2y, REAL8 s2z, PrecessingNRSurData *__sur_data)
Computes the starting frequency of the NRSur7dq2 or NRSur7dq4 waveform approximant when evaluated usi...
REAL8 NRSur7dq2_eval_fit(FitData *data, REAL8 *x)
static REAL8 cubic_interp(REAL8 xout, REAL8 *x, REAL8 *y)
Cubic interpolation of 4 data points This gives a much closer result to scipy.interpolate....
static void NRSur7dq4_LoadVectorFitData(VectorFitData **vector_fit_data, LALH5File *sub, const char *name, const size_t size)
Loads a vector fit for NRSur7dq4.
SphHarmTimeSeries * XLALSphHarmTimeSeriesAddMode(SphHarmTimeSeries *appended, const COMPLEX16TimeSeries *inmode, UINT4 l, INT4 m)
Prepend a node to a linked list of SphHarmTimeSeries, or create a new head.
COMPLEX16 XLALSpinWeightedSphericalHarmonic(REAL8 theta, REAL8 phi, int s, int l, int m)
REAL8TimeSeries * XLALCreateREAL8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
void XLALDestroyCOMPLEX16TimeSeries(COMPLEX16TimeSeries *series)
COMPLEX16TimeSeries * XLALCreateCOMPLEX16TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
const LALUnit lalStrainUnit
#define XLAL_ERROR_VOID(...)
#define XLAL_ERROR_REAL8(...)
#define XLAL_PRINT_INFO(...)
#define XLAL_ERROR_NULL(...)
#define XLAL_PRINT_WARNING(...)
#define XLAL_CHECK_ABORT(assertion)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_IS_REAL8_FAIL_NAN(val)
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
int pthread_once(pthread_once_t *once_control, void(*init_routine)(void))
Data for a single dynamics node.
FitData * omega_data
A fit to the orbital angular frequency.
VectorFitData * chiB_dot_data
A 3d vector fit for the coorbital components of the time derivative of chiB taken in the coprecessing...
VectorFitData * chiA_dot_data
A 3d vector fit for the coorbital components of the time derivative of chiA taken in the coprecessing...
VectorFitData * omega_copr_data
A 2d vector fit for the x and y components of Omega^{coorb}(t) in arxiv 1705.07089.
Data used in a single scalar fit.
int n_coefs
Number of coefficients in the fit.
gsl_matrix_long * basisFunctionOrders
matrix of (n_coefs x 7) basis function orders giving the polynomial order in f(q),...
gsl_vector * coefs
coefficient vector of length n_coefs
All data needed by the full surrogate model.
gsl_vector * t_coorb
Vector of the coorbital surrogate output times.
UINT4 PrecessingNRSurVersion
One for each 2 <= ell <= LMax.
WaveformFixedEllModeData ** coorbital_mode_data
A DynamicsNodeFitData for each time in t_ds_half_times.
DynamicsNodeFitData ** ds_half_node_data
A DynamicsNodeFitData for each time in t_ds.
UINT4 setup
Indicates if this has been initialized.
DynamicsNodeFitData ** ds_node_data
Structure to carry a collection of spherical harmonic modes in COMPLEX16 time series.
Data used in a single vector fit NOTE: basisFunctionOrders, coefs, componentIndices,...
int n_coefs
Number of coefficients in the fit.
gsl_matrix_long * basisFunctionOrders
matrix of (n_coefs x 7) basis function orders giving the polynomial order in f(q),...
int vec_dim
Dimension of the vector.
gsl_vector_long * componentIndices
Each fit coefficient applies to a single component of the vector; this gives the component indices.
gsl_vector * coefs
coefficient vector of length n_coefs