22 #define UNUSED __attribute__ ((unused))
37 #include <gsl/gsl_errno.h>
38 #include <gsl/gsl_bspline.h>
39 #include <gsl/gsl_blas.h>
40 #include <gsl/gsl_min.h>
41 #include <gsl/gsl_spline.h>
42 #include <gsl/gsl_poly.h>
43 #include <lal/Units.h>
44 #include <lal/SeqFactories.h>
45 #include <lal/LALConstants.h>
46 #include <lal/XLALError.h>
47 #include <lal/FrequencySeries.h>
49 #include <lal/StringInput.h>
50 #include <lal/Sequence.h>
51 #include <lal/LALStdio.h>
52 #include <lal/FileIO.h>
55 #ifdef LAL_HDF5_ENABLED
56 #include <lal/H5FileIO.h>
57 static const char ROMDataHDF5[] =
"SEOBNRv4ROM_v2.0.hdf5";
58 static const INT4 ROMDataHDF5_VERSION_MAJOR = 2;
59 static const INT4 ROMDataHDF5_VERSION_MINOR = 0;
60 static const INT4 ROMDataHDF5_VERSION_MICRO = 0;
63 #include <lal/LALSimInspiral.h>
64 #include <lal/LALSimIMR.h>
68 #include <lal/LALConfig.h>
69 #ifdef LAL_PTHREAD_LOCK
75 #ifdef LAL_PTHREAD_LOCK
76 static pthread_once_t SEOBNRv4ROM_is_initialized = PTHREAD_ONCE_INIT;
81 typedef struct tagSEOBNRROMdataDS_coeff
110 SEOBNRROMdataDS_submodel*
sub1;
111 SEOBNRROMdataDS_submodel*
sub2;
112 SEOBNRROMdataDS_submodel*
sub3;
118 typedef int (*
load_dataPtr)(
const char*, gsl_vector *, gsl_vector *, gsl_matrix *, gsl_matrix *, gsl_vector *);
120 typedef struct tagSplineData
122 gsl_bspline_workspace *bwx;
123 gsl_bspline_workspace *bwy;
124 gsl_bspline_workspace *bwz;
140 gsl_vector *cvec_amp,
141 gsl_vector *cvec_phi,
149 const double *etavec,
150 const double *chi1vec,
151 const double *chi2vec,
158 UNUSED SEOBNRROMdataDS_submodel **submodel,
159 UNUSED
const char dir[],
160 UNUSED
const char grp_name[]
195 static size_t NextPow2(
const size_t n);
202 const double *etavec,
203 const double *chi1vec,
204 const double *chi2vec
208 gsl_spline **spline_phi,
209 gsl_interp_accel **acc_phi,
226 gsl_bspline_workspace *bwx,
227 gsl_bspline_workspace *bwy
232 SEOBNRROMdataDS_submodel *submodel_lo,
233 SEOBNRROMdataDS_submodel *submodel_hi,
234 gsl_vector* amp_f_lo,
235 gsl_vector* amp_f_hi,
240 gsl_interp_accel **acc_amp,
241 gsl_spline **spline_amp
246 SEOBNRROMdataDS_submodel *submodel_lo,
247 SEOBNRROMdataDS_submodel *submodel_hi,
248 gsl_vector* phi_f_lo,
249 gsl_vector* phi_f_hi,
252 gsl_interp_accel **acc_phi_out,
253 gsl_spline **spline_phi_out
291 const double *etavec,
292 const double *chi1vec,
293 const double *chi2vec
296 if(!splinedata) exit(1);
302 const size_t nbreak_x = ncx-2;
303 const size_t nbreak_y = ncy-2;
304 const size_t nbreak_z = ncz-2;
307 gsl_bspline_workspace *bwx = gsl_bspline_alloc(4, nbreak_x);
308 gsl_bspline_workspace *bwy = gsl_bspline_alloc(4, nbreak_y);
309 gsl_bspline_workspace *bwz = gsl_bspline_alloc(4, nbreak_z);
312 gsl_vector *breakpts_x = gsl_vector_alloc(nbreak_x);
313 gsl_vector *breakpts_y = gsl_vector_alloc(nbreak_y);
314 gsl_vector *breakpts_z = gsl_vector_alloc(nbreak_z);
316 gsl_vector_set(breakpts_x,
i, etavec[
i]);
317 for (
UINT4 j=0; j<nbreak_y; j++)
318 gsl_vector_set(breakpts_y, j, chi1vec[j]);
319 for (
UINT4 k=0; k<nbreak_z; k++)
320 gsl_vector_set(breakpts_z, k, chi2vec[k]);
322 gsl_bspline_knots(breakpts_x, bwx);
323 gsl_bspline_knots(breakpts_y, bwy);
324 gsl_bspline_knots(breakpts_z, bwz);
326 gsl_vector_free(breakpts_x);
327 gsl_vector_free(breakpts_y);
328 gsl_vector_free(breakpts_z);
330 (*splinedata)->bwx=bwx;
331 (*splinedata)->bwy=bwy;
332 (*splinedata)->bwz=bwz;
337 if(!splinedata)
return;
338 if(splinedata->
bwx) gsl_bspline_free(splinedata->
bwx);
339 if(splinedata->
bwy) gsl_bspline_free(splinedata->
bwy);
340 if(splinedata->
bwz) gsl_bspline_free(splinedata->
bwz);
350 gsl_vector *cvec_amp,
351 gsl_vector *cvec_phi,
358 const double *etavec,
359 const double *chi1vec,
360 const double *chi2vec,
376 gsl_bspline_workspace *bwx=splinedata->
bwx;
377 gsl_bspline_workspace *bwy=splinedata->
bwy;
378 gsl_bspline_workspace *bwz=splinedata->
bwz;
382 for (
int k=0; k<
nk_amp; k++) {
383 gsl_vector v = gsl_vector_subvector(cvec_amp, k*
N,
N).vector;
385 gsl_vector_set(c_amp, k, csum);
389 for (
int k=0; k<
nk_phi; k++) {
390 gsl_vector v = gsl_vector_subvector(cvec_phi, k*
N,
N).vector;
392 gsl_vector_set(c_phi, k, csum);
402 SEOBNRROMdataDS_submodel **submodel,
403 UNUSED
const char dir[],
404 UNUSED
const char grp_name[]
408 if(!submodel) exit(1);
411 *submodel =
XLALCalloc(1,
sizeof(SEOBNRROMdataDS_submodel));
415 #ifdef LAL_HDF5_ENABLED
416 size_t size = strlen(dir) + strlen(ROMDataHDF5) + 2;
418 snprintf(
path, size,
"%s/%s", dir, ROMDataHDF5);
424 ReadHDF5RealVectorDataset(sub,
"Amp_ciall", & (*submodel)->cvec_amp);
425 ReadHDF5RealVectorDataset(sub,
"Phase_ciall", & (*submodel)->cvec_phi);
428 ReadHDF5RealMatrixDataset(sub,
"Bamp", & (*submodel)->Bamp);
429 ReadHDF5RealMatrixDataset(sub,
"Bphase", & (*submodel)->Bphi);
432 ReadHDF5RealVectorDataset(sub,
"Mf_grid_Amp", & (*submodel)->gA);
433 ReadHDF5RealVectorDataset(sub,
"Mf_grid_Phi", & (*submodel)->gPhi);
436 ReadHDF5RealVectorDataset(sub,
"etavec", & (*submodel)->etavec);
437 ReadHDF5RealVectorDataset(sub,
"chi1vec", & (*submodel)->chi1vec);
438 ReadHDF5RealVectorDataset(sub,
"chi2vec", & (*submodel)->chi2vec);
441 (*submodel)->nk_amp = (*submodel)->gA->size;
442 (*submodel)->nk_phi = (*submodel)->gPhi->size;
443 (*submodel)->ncx = (*submodel)->etavec->size + 2;
444 (*submodel)->ncy = (*submodel)->chi1vec->size + 2;
445 (*submodel)->ncz = (*submodel)->chi2vec->size + 2;
448 (*submodel)->eta_bounds[0] = gsl_vector_get((*submodel)->etavec, 0);
449 (*submodel)->eta_bounds[1] = gsl_vector_get((*submodel)->etavec, (*submodel)->etavec->size - 1);
450 (*submodel)->chi1_bounds[0] = gsl_vector_get((*submodel)->chi1vec, 0);
451 (*submodel)->chi1_bounds[1] = gsl_vector_get((*submodel)->chi1vec, (*submodel)->chi1vec->size - 1);
452 (*submodel)->chi2_bounds[0] = gsl_vector_get((*submodel)->chi2vec, 0);
453 (*submodel)->chi2_bounds[1] = gsl_vector_get((*submodel)->chi2vec, (*submodel)->chi2vec->size - 1);
467 if(submodel->cvec_amp) gsl_vector_free(submodel->cvec_amp);
468 if(submodel->cvec_phi) gsl_vector_free(submodel->cvec_phi);
469 if(submodel->Bamp) gsl_matrix_free(submodel->Bamp);
470 if(submodel->Bphi) gsl_matrix_free(submodel->Bphi);
471 if(submodel->gA) gsl_vector_free(submodel->gA);
472 if(submodel->gPhi) gsl_vector_free(submodel->gPhi);
473 if(submodel->etavec) gsl_vector_free(submodel->etavec);
474 if(submodel->chi1vec) gsl_vector_free(submodel->chi1vec);
475 if(submodel->chi2vec) gsl_vector_free(submodel->chi2vec);
480 UNUSED SEOBNRROMdataDS *romdata,
481 UNUSED
const char dir[])
487 XLALPrintError(
"WARNING: You tried to setup the SEOBNRv4ROM model that was already initialised. Ignoring\n");
491 #ifdef LAL_HDF5_ENABLED
493 size_t size = strlen(dir) + strlen(ROMDataHDF5) + 2;
495 snprintf(
path, size,
"%s/%s", dir, ROMDataHDF5);
499 PrintInfoStringAttribute(
file,
"Email");
500 PrintInfoStringAttribute(
file,
"Description");
501 ret = ROM_check_version_number(
file, ROMDataHDF5_VERSION_MAJOR,
502 ROMDataHDF5_VERSION_MINOR,
503 ROMDataHDF5_VERSION_MICRO);
532 (romdata)->sub1 = NULL;
535 (romdata)->sub2 = NULL;
538 (romdata)->sub3 = NULL;
544 if(!romdatacoeff) exit(1);
551 (*romdatacoeff)->c_amp = gsl_vector_alloc(
nk_amp);
552 (*romdatacoeff)->c_phi = gsl_vector_alloc(
nk_phi);
557 if(romdatacoeff->
c_amp) gsl_vector_free(romdatacoeff->
c_amp);
558 if(romdatacoeff->
c_phi) gsl_vector_free(romdatacoeff->
c_phi);
565 return 1 << (size_t) ceil(log2(n));
570 SEOBNRROMdataDS_submodel *submodel_lo,
571 SEOBNRROMdataDS_submodel *submodel_hi,
572 gsl_vector* amp_f_lo,
573 gsl_vector* amp_f_hi,
579 gsl_interp_accel **acc_amp,
580 gsl_spline **spline_amp
585 for (jA_lo=0; jA_lo < submodel_lo->nk_amp; jA_lo++) {
586 if (gsl_vector_get(submodel_lo->gA, jA_lo) > Mfm) {
594 for (jA_hi=0; jA_hi < submodel_hi->nk_amp; jA_hi++)
595 if (gsl_vector_get(submodel_hi->gA, jA_hi) > Mfm)
598 int nA = 1 + jA_lo + (submodel_hi->nk_amp - jA_hi);
600 gsl_vector *gAU = gsl_vector_alloc(nA);
601 gsl_vector *amp_f = gsl_vector_alloc(nA);
604 for (
int i=0;
i<=jA_lo;
i++) {
605 gsl_vector_set(gAU,
i, gsl_vector_get(submodel_lo->gA,
i));
606 double A = amp_pre_lo * gsl_vector_get(amp_f_lo,
i);
607 gsl_vector_set(amp_f,
i, A);
610 for (
int i=jA_lo+1;
i<nA;
i++) {
611 int k = jA_hi - (jA_lo+1) +
i;
612 gsl_vector_set(gAU,
i, gsl_vector_get(submodel_hi->gA, k));
613 double A = amp_pre_hi * gsl_vector_get(amp_f_hi, k);
614 gsl_vector_set(amp_f,
i, A);
618 *acc_amp = gsl_interp_accel_alloc();
619 *spline_amp = gsl_spline_alloc(gsl_interp_cspline, nA);
621 gsl_spline_init(*spline_amp, gsl_vector_const_ptr(gAU,0),
622 gsl_vector_const_ptr(amp_f,0), nA);
624 gsl_vector_free(gAU);
625 gsl_vector_free(amp_f);
626 gsl_vector_free(amp_f_lo);
627 gsl_vector_free(amp_f_hi);
633 SEOBNRROMdataDS_submodel *submodel_lo,
634 SEOBNRROMdataDS_submodel *submodel_hi,
635 gsl_vector* phi_f_lo,
636 gsl_vector* phi_f_hi,
639 gsl_interp_accel **acc_phi,
640 gsl_spline **spline_phi
645 for (jP_lo=0; jP_lo < submodel_lo->nk_phi; jP_lo++) {
646 if (gsl_vector_get(submodel_lo->gPhi, jP_lo) > Mfm) {
654 for (jP_hi=0; jP_hi < submodel_hi->nk_phi; jP_hi++)
655 if (gsl_vector_get(submodel_hi->gPhi, jP_hi) > Mfm)
658 int nP = 1 + jP_lo + (submodel_hi->nk_phi - jP_hi);
659 gsl_vector *gPU = gsl_vector_alloc(nP);
660 gsl_vector *phi_f = gsl_vector_alloc(nP);
662 for (
int i=0;
i<=jP_lo;
i++) {
663 gsl_vector_set(gPU,
i, gsl_vector_get(submodel_lo->gPhi,
i));
664 double P = gsl_vector_get(phi_f_lo,
i);
665 gsl_vector_set(phi_f,
i, P);
668 for (
int i=jP_lo+1;
i<nP;
i++) {
669 int k = jP_hi - (jP_lo+1) +
i;
670 gsl_vector_set(gPU,
i, gsl_vector_get(submodel_hi->gPhi, k));
671 double P = gsl_vector_get(phi_f_hi, k);
672 gsl_vector_set(phi_f,
i, P);
680 gsl_interp_accel *acc_phi_lo = gsl_interp_accel_alloc();
681 gsl_spline *spline_phi_lo = gsl_spline_alloc(gsl_interp_cspline, submodel_lo->nk_phi);
682 gsl_spline_init(spline_phi_lo, gsl_vector_const_ptr(submodel_lo->gPhi,0),
683 gsl_vector_const_ptr(phi_f_lo,0), submodel_lo->nk_phi);
686 gsl_vector_const_view gP_hi_data = gsl_vector_const_subvector(submodel_hi->gPhi, jP_hi - nn, 2*nn+1);
687 gsl_vector_const_view P_hi_data = gsl_vector_const_subvector(phi_f_hi, jP_hi - nn, 2*nn+1);
688 gsl_vector *P_lo_data = gsl_vector_alloc(2*nn+1);
689 for (
int i=0;
i<2*nn+1;
i++) {
690 double P = gsl_spline_eval(spline_phi_lo, gsl_vector_get(&gP_hi_data.vector,
i), acc_phi_lo);
691 gsl_vector_set(P_lo_data,
i, P);
695 gsl_vector *cP_lo =
Fit_cubic(&gP_hi_data.vector, P_lo_data);
696 gsl_vector *cP_hi =
Fit_cubic(&gP_hi_data.vector, &P_hi_data.vector);
698 double P_lo_derivs[2];
699 double P_hi_derivs[2];
700 gsl_poly_eval_derivs(cP_lo->data, 4, Mfm, P_lo_derivs, 2);
701 gsl_poly_eval_derivs(cP_hi->data, 4, Mfm, P_hi_derivs, 2);
703 double delta_omega = P_hi_derivs[1] - P_lo_derivs[1];
704 double delta_phi = P_hi_derivs[0] - P_lo_derivs[0] - delta_omega * Mfm;
706 for (
int i=jP_lo+1;
i<nP;
i++) {
707 int k = jP_hi - (jP_lo+1) +
i;
708 double f = gsl_vector_get(submodel_hi->gPhi, k);
709 gsl_vector_set(gPU,
i, f);
710 double P = gsl_vector_get(phi_f_hi, k) - delta_omega * f - delta_phi;
711 gsl_vector_set(phi_f,
i, P);
715 gsl_vector_free(P_lo_data);
716 gsl_vector_free(cP_lo);
717 gsl_vector_free(cP_hi);
718 gsl_vector_free(phi_f_lo);
719 gsl_vector_free(phi_f_hi);
722 *acc_phi = gsl_interp_accel_alloc();
723 *spline_phi = gsl_spline_alloc(gsl_interp_cspline, nP);
725 gsl_spline_init(*spline_phi, gsl_vector_const_ptr(gPU,0),
726 gsl_vector_const_ptr(phi_f,0), nP);
730 gsl_vector_free(phi_f);
731 gsl_vector_free(gPU);
732 gsl_spline_free(spline_phi_lo);
733 gsl_interp_accel_free(acc_phi_lo);
767 if(!hptilde || !hctilde)
773 "Error setting up SEOBNRv4ROM data - check your $LAL_DATA_PATH\n");
776 if(*hptilde || *hctilde) {
777 XLALPrintError(
"(*hptilde) and (*hctilde) are supposed to be NULL, but got %p and %p",
778 (*hptilde), (*hctilde));
784 if (eta > 0.25)
nudge(&eta, 0.25, 1
e-6);
785 if (eta < 0.01)
nudge(&eta, 0.01, 1
e-6);
787 if ( chi1 < -1.0 || chi2 < -1.0 || chi1 > 1.0 || chi2 > 1.0) {
788 XLALPrintError(
"XLAL Error - %s: chi1 or chi2 smaller than -1.0 or larger than 1.0!\n"
789 "SEOBNRv4ROM is only available for spins in the range -1 <= a/M <= 1.0.\n",
794 if (eta<0.01 || eta > 0.25) {
795 XLALPrintError(
"XLAL Error - %s: eta (%f) smaller than 0.01 or unphysical!\n"
796 "SEOBNRv4ROM is only available for eta in the range 0.01 <= eta <= 0.25.\n",
802 SEOBNRROMdataDS_submodel *submodel_hi;
803 SEOBNRROMdataDS_submodel *submodel_lo;
804 submodel_lo = romdata->sub1;
807 if (chi1 < romdata->sub3->chi1_bounds[0] || eta > romdata->sub3->eta_bounds[1])
808 submodel_hi = romdata->sub2;
810 submodel_hi = romdata->sub3;
815 double fLow = freqs_in->
data[0];
816 double fHigh = freqs_in->
data[freqs_in->
length - 1];
823 double Mf_ROM_min = fmax(gsl_vector_get(submodel_lo->gA, 0),
824 gsl_vector_get(submodel_lo->gPhi,0));
826 double Mf_ROM_max = fmin(gsl_vector_get(submodel_hi->gA, submodel_hi->nk_amp-1),
827 gsl_vector_get(submodel_hi->gPhi, submodel_hi->nk_phi-1));
828 double fLow_geom = fLow * Mtot_sec;
829 double fHigh_geom = fHigh * Mtot_sec;
830 double fRef_geom = fRef * Mtot_sec;
831 double deltaF_geom = deltaF * Mtot_sec;
836 if (fHigh_geom == 0 || fHigh_geom >
Mf_ROM_max)
840 if (fHigh_geom <= fLow_geom)
841 XLAL_ERROR(
XLAL_EDOM,
"End frequency %g is smaller than (or equal to) starting frequency %g!\n", fHigh_geom, fLow_geom);
843 XLALPrintWarning(
"Reference frequency Mf_ref=%g is greater than maximal frequency in ROM Mf=%g. Starting at maximal frequency in ROM.\n", fRef_geom,
Mf_ROM_max);
847 XLALPrintWarning(
"Reference frequency Mf_ref=%g is smaller than lowest frequency in ROM Mf=%g. Starting at lowest frequency in ROM.\n", fLow_geom,
Mf_ROM_min);
859 REAL8 amp_pre_lo = 1.0;
860 REAL8 amp_pre_hi = 1.0;
867 submodel_lo->cvec_amp,
868 submodel_lo->cvec_phi,
875 gsl_vector_const_ptr(submodel_lo->etavec, 0),
876 gsl_vector_const_ptr(submodel_lo->chi1vec, 0),
877 gsl_vector_const_ptr(submodel_lo->chi2vec, 0),
878 romdata_coeff_lo->
c_amp,
879 romdata_coeff_lo->
c_phi
892 submodel_hi->cvec_amp,
893 submodel_hi->cvec_phi,
900 gsl_vector_const_ptr(submodel_hi->etavec, 0),
901 gsl_vector_const_ptr(submodel_hi->chi1vec, 0),
902 gsl_vector_const_ptr(submodel_hi->chi2vec, 0),
903 romdata_coeff_hi->
c_amp,
904 romdata_coeff_hi->
c_phi
916 gsl_vector* amp_f_lo = gsl_vector_alloc(submodel_lo->nk_amp);
917 gsl_vector* phi_f_lo = gsl_vector_alloc(submodel_lo->nk_phi);
918 gsl_blas_dgemv(CblasTrans, 1.0, submodel_lo->Bamp, romdata_coeff_lo->
c_amp, 0.0, amp_f_lo);
919 gsl_blas_dgemv(CblasTrans, 1.0, submodel_lo->Bphi, romdata_coeff_lo->
c_phi, 0.0, phi_f_lo);
921 gsl_vector* amp_f_hi = gsl_vector_alloc(submodel_hi->nk_amp);
922 gsl_vector* phi_f_hi = gsl_vector_alloc(submodel_hi->nk_phi);
923 gsl_blas_dgemv(CblasTrans, 1.0, submodel_hi->Bamp, romdata_coeff_hi->
c_amp, 0.0, amp_f_hi);
924 gsl_blas_dgemv(CblasTrans, 1.0, submodel_hi->Bphi, romdata_coeff_hi->
c_phi, 0.0, phi_f_hi);
926 const double Mfm = 0.01;
929 gsl_interp_accel *acc_amp;
930 gsl_spline *spline_amp;
931 GlueAmplitude(submodel_lo, submodel_hi, amp_f_lo, amp_f_hi, amp_pre_lo, amp_pre_hi, Mfm,
932 &acc_amp, &spline_amp
936 gsl_interp_accel *acc_phi;
937 gsl_spline *spline_phi;
938 GluePhasing(submodel_lo, submodel_hi, phi_f_lo, phi_f_hi, Mfm,
939 &acc_phi, &spline_phi
949 npts =
NextPow2(fHigh_geom / deltaF_geom) + 1;
950 if (fHigh_geom < fHigh * Mtot_sec)
951 npts =
NextPow2(fHigh * Mtot_sec / deltaF_geom) + 1;
959 double fHigh_temp = fHigh_geom / Mtot_sec;
961 UINT4 iStop = (
UINT4) ceil(fHigh_temp / deltaF);
967 freqs->
data[
i-iStart] =
i*deltaF_geom;
981 freqs->
data[
i] = freqs_in->
data[
i] * Mtot_sec;
984 if (!(*hptilde) || !(*hctilde)) {
986 gsl_spline_free(spline_amp);
987 gsl_spline_free(spline_phi);
988 gsl_interp_accel_free(acc_amp);
989 gsl_interp_accel_free(acc_phi);
994 memset((*hptilde)->data->data, 0, npts *
sizeof(
COMPLEX16));
995 memset((*hctilde)->data->data, 0, npts *
sizeof(
COMPLEX16));
1000 COMPLEX16 *pdata=(*hptilde)->data->data;
1001 COMPLEX16 *cdata=(*hctilde)->data->data;
1003 REAL8 cosi = cos(inclination);
1004 REAL8 pcoef = 0.5*(1.0 + cosi*cosi);
1009 double amp0 = Mtot * Mtot_sec *
LAL_MRSUN_SI / (distance);
1012 double phase_change = gsl_spline_eval(spline_phi, fRef_geom, acc_phi) - 2*phiRef;
1018 const REAL8 factor = sqrt(1. - 4.*eta);
1019 const REAL8 m1 = 0.5*Mtot*(1.+ factor);
1020 const REAL8 m2 = 0.5*Mtot*(1.- factor);
1027 XLAL_CHECK(
XLAL_SUCCESS == ret, ret,
"Failed to generate tidal amplitude series to construct SEOBNRv4_ROM_NRTidalv2 waveform.");
1030 double f = freqs->
data[
i];
1031 double ampT = amp_tidal->data[
i];
1035 double A = gsl_spline_eval(spline_amp, f, acc_amp);
1036 double phase = gsl_spline_eval(spline_phi, f, acc_phi) - phase_change;
1037 COMPLEX16 htilde =
s*amp0*(A+ampT) * (cos(phase) + I*sin(phase));
1038 pdata[j] = pcoef * htilde;
1039 cdata[j] = -I * ccoef * htilde;
1043 double f = freqs->
data[
i];
1046 double A = gsl_spline_eval(spline_amp, f, acc_amp);
1047 double phase = gsl_spline_eval(spline_phi, f, acc_phi) - phase_change;
1048 COMPLEX16 htilde =
s*amp0*A * (cos(phase) + I*sin(phase));
1050 pdata[j] = pcoef * htilde;
1051 cdata[j] = -I * ccoef * htilde;
1069 gsl_spline_free(spline_amp);
1070 gsl_spline_free(spline_phi);
1071 gsl_interp_accel_free(acc_amp);
1072 gsl_interp_accel_free(acc_phi);
1080 REAL8 t_corr = gsl_spline_eval_deriv(spline_phi, Mf_final, acc_phi) / (2*
LAL_PI);
1084 double f = freqs->
data[
i] - fRef_geom;
1086 double phase_factor = -2*
LAL_PI * f * t_corr;
1087 COMPLEX16 t_factor = (cos(phase_factor) + I*sin(phase_factor));
1088 pdata[j] *= t_factor;
1089 cdata[j] *= t_factor;
1095 gsl_spline_free(spline_amp);
1096 gsl_spline_free(spline_phi);
1097 gsl_interp_accel_free(acc_amp);
1098 gsl_interp_accel_free(acc_phi);
1157 struct tagCOMPLEX16FrequencySeries **hptilde,
1158 struct tagCOMPLEX16FrequencySeries **hctilde,
1176 double m1temp = m1SI;
1177 double chi1temp = chi1;
1187 double Mtot = mass1+mass2;
1188 double eta = mass1 * mass2 / (Mtot*Mtot);
1194 #ifdef LAL_PTHREAD_LOCK
1202 "Error setting up SEOBNRv4ROM data - check your $LAL_DATA_PATH\n");
1207 int retcode =
SEOBNRv4ROMCore(hptilde, hctilde, phiRef, fRef, distance,
1208 inclination, Mtot_sec, eta, chi1, chi2, freqs,
1209 0, nk_max, LALparams, NRTidal_version);
1222 struct tagCOMPLEX16FrequencySeries **hptilde,
1223 struct tagCOMPLEX16FrequencySeries **hctilde,
1243 double m1temp = m1SI;
1244 double chi1temp = chi1;
1254 double Mtot = mass1+mass2;
1255 double eta = mass1 * mass2 / (Mtot*Mtot);
1262 #ifdef LAL_PTHREAD_LOCK
1272 freqs->
data[0] = fLow;
1273 freqs->
data[1] = fHigh;
1275 int retcode =
SEOBNRv4ROMCore(hptilde, hctilde, phiRef, fRef, distance,
1276 inclination, Mtot_sec, eta, chi1, chi2, freqs,
1277 deltaF, nk_max, LALparams, NRTidal_version);
1288 gsl_spline **spline_phi,
1289 gsl_interp_accel **acc_phi,
1303 double Mtot = mass1 + mass2;
1304 double eta = mass1 * mass2 / (Mtot*Mtot);
1308 if (eta > 0.25)
nudge(&eta, 0.25, 1
e-6);
1309 if (eta < 0.01)
nudge(&eta, 0.01, 1
e-6);
1311 if ( chi1 < -1.0 || chi2 < -1.0 || chi1 > 1.0 || chi2 > 1.0) {
1312 XLALPrintError(
"XLAL Error - %s: chi1 or chi2 smaller than -1.0 or larger than 1.0!\nSEOBNRv4ROM is only available for spins in the range -1 <= a/M <= 1.0.\n", __func__);
1316 if (eta < 0.01 || eta > 0.25) {
1317 XLALPrintError(
"XLAL Error - %s: eta (%f) smaller than 0.01 or unphysical!\nSEOBNRv4ROM is only available for symmetric mass-ratios in the range 0.01 <= eta <= 0.25.\n", __func__,eta);
1322 #ifdef LAL_PTHREAD_LOCK
1331 "Error setting up SEOBNRv4ROM data - check your $LAL_DATA_PATH\n");
1335 SEOBNRROMdataDS_submodel *submodel_hi;
1336 SEOBNRROMdataDS_submodel *submodel_lo;
1337 submodel_lo = romdata->sub1;
1340 if (chi1 < romdata->sub3->chi1_bounds[0] || eta > romdata->sub3->eta_bounds[1])
1341 submodel_hi = romdata->sub2;
1343 submodel_hi = romdata->sub3;
1349 submodel_lo->nk_phi);
1351 submodel_hi->nk_phi);
1354 *
Mf_ROM_min = fmax(gsl_vector_get(submodel_lo->gA, 0),
1355 gsl_vector_get(submodel_lo->gPhi,0));
1357 *
Mf_ROM_max = fmin(gsl_vector_get(submodel_hi->gA, submodel_hi->nk_amp-1),
1358 gsl_vector_get(submodel_hi->gPhi, submodel_hi->nk_phi-1));
1367 submodel_lo->cvec_amp,
1368 submodel_lo->cvec_phi,
1369 submodel_lo->nk_amp,
1370 submodel_lo->nk_phi,
1375 gsl_vector_const_ptr(submodel_lo->etavec, 0),
1376 gsl_vector_const_ptr(submodel_lo->chi1vec, 0),
1377 gsl_vector_const_ptr(submodel_lo->chi2vec, 0),
1378 romdata_coeff_lo->
c_amp,
1379 romdata_coeff_lo->
c_phi
1392 submodel_hi->cvec_amp,
1393 submodel_hi->cvec_phi,
1394 submodel_hi->nk_amp,
1395 submodel_hi->nk_phi,
1400 gsl_vector_const_ptr(submodel_hi->etavec, 0),
1401 gsl_vector_const_ptr(submodel_hi->chi1vec, 0),
1402 gsl_vector_const_ptr(submodel_hi->chi2vec, 0),
1403 romdata_coeff_hi->
c_amp,
1404 romdata_coeff_hi->
c_phi
1414 gsl_vector* phi_f_lo = gsl_vector_alloc(submodel_lo->nk_phi);
1415 gsl_blas_dgemv(CblasTrans, 1.0, submodel_lo->Bphi, romdata_coeff_lo->
c_phi,
1418 gsl_vector* phi_f_hi = gsl_vector_alloc(submodel_hi->nk_phi);
1419 gsl_blas_dgemv(CblasTrans, 1.0, submodel_hi->Bphi, romdata_coeff_hi->
c_phi,
1422 const double Mfm = 0.01;
1425 GluePhasing(submodel_lo, submodel_hi, phi_f_lo, phi_f_hi, Mfm,
1464 double m1temp = m1SI;
1465 double chi1temp = chi1;
1473 gsl_spline *spline_phi = NULL;
1474 gsl_interp_accel *acc_phi = NULL;
1475 double Mf_final = NAN, Mtot_sec;
1478 &Mtot_sec, m1SI, m2SI, chi1, chi2,
1483 if (spline_phi == NULL) {
1484 XLALPrintError(
"XLAL Error - %s: `spline_phi` is not initialized\n", __func__);
1488 if (acc_phi == NULL) {
1489 XLALPrintError(
"XLAL Error - %s: `acc_phi` is not initialized\n", __func__);
1493 if (isnan(Mf_final)) {
1494 XLALPrintError(
"XLAL Error - %s: `Mf_final` is not initialized\n", __func__);
1499 XLALPrintError(
"XLAL Error - %s: `Mf_ROM_min` is not initialized\n", __func__);
1504 XLALPrintError(
"XLAL Error - %s: `Mf_ROM_max` is not initialized\n", __func__);
1509 double t_corr = gsl_spline_eval_deriv(spline_phi, Mf_final, acc_phi) / (2*
LAL_PI);
1512 double Mf = frequency * Mtot_sec;
1513 if (Mf < Mf_ROM_min || Mf >
Mf_ROM_max || Mf > Mf_final) {
1514 gsl_spline_free(spline_phi);
1515 gsl_interp_accel_free(acc_phi);
1517 "Min / max / final Mf values are %g, %g, %g\n", frequency, Mf,
Mf_ROM_min,
Mf_ROM_max, Mf_final);
1521 double time_M = gsl_spline_eval_deriv(spline_phi, frequency * Mtot_sec, acc_phi) / (2*
LAL_PI) - t_corr;
1522 *t = time_M * Mtot_sec;
1524 gsl_spline_free(spline_phi);
1525 gsl_interp_accel_free(acc_phi);
1557 double m1temp = m1SI;
1558 double chi1temp = chi1;
1566 gsl_spline *spline_phi = NULL;
1567 gsl_interp_accel *acc_phi = NULL;
1568 double Mf_final = NAN, Mtot_sec;
1571 &Mtot_sec, m1SI, m2SI, chi1, chi2,
1576 if (spline_phi == NULL) {
1577 XLALPrintError(
"XLAL Error - %s: `spline_phi` is not initialized\n", __func__);
1581 if (acc_phi == NULL) {
1582 XLALPrintError(
"XLAL Error - %s: `acc_phi` is not initialized\n", __func__);
1586 if (isnan(Mf_final)) {
1587 XLALPrintError(
"XLAL Error - %s: `Mf_final` is not initialized\n", __func__);
1592 XLALPrintError(
"XLAL Error - %s: `Mf_ROM_min` is not initialized\n", __func__);
1597 XLALPrintError(
"XLAL Error - %s: `Mf_ROM_max` is not initialized\n", __func__);
1602 double t_corr = gsl_spline_eval_deriv(spline_phi, Mf_final, acc_phi) / (2*
LAL_PI);
1608 double log_f_pts[
N];
1609 double log_t_pts[
N];
1610 double log_f_min = log(
Mf_ROM_min * 1.000001);
1611 double log_f_rng_2 = log(Mf_final/2.0);
1612 double dlog_f = (log_f_rng_2 - log_f_min) / (
N-1);
1615 for (
int i=0;
i<
N;
i++) {
1616 log_f_pts[
i] = log_f_rng_2 -
i*dlog_f;
1618 double time_M = gsl_spline_eval_deriv(spline_phi, exp(log_f_pts[
i]), acc_phi) / (2*
LAL_PI) - t_corr;
1619 log_t_pts[
i] = log(time_M * Mtot_sec);
1623 double t_rng_2 = exp(log_t_pts[0]);
1624 double t_min = exp(log_t_pts[
N-1]);
1625 if (t < t_rng_2 || t > t_min) {
1626 gsl_spline_free(spline_phi);
1627 gsl_interp_accel_free(acc_phi);
1628 XLAL_ERROR(
XLAL_EDOM,
"The frequency of time %g is outside allowed frequency range.\n", t);
1632 gsl_interp_accel *acc = gsl_interp_accel_alloc();
1633 gsl_spline *spline = gsl_spline_alloc(gsl_interp_cspline,
N);
1634 gsl_spline_init(spline, log_t_pts, log_f_pts,
N);
1636 *frequency = exp(gsl_spline_eval(spline, log(t), acc)) / Mtot_sec;
1638 gsl_spline_free(spline);
1639 gsl_interp_accel_free(acc);
1640 gsl_spline_free(spline_phi);
1641 gsl_interp_accel_free(acc_phi);
1654 #ifdef LAL_HDF5_ENABLED
1655 #define datafile ROMDataHDF5
1659 char *dir = dirname(
path);
struct tagLALH5File LALH5File
int XLALSimNRTunedTidesFDTidalAmplitudeFrequencySeries(const REAL8Sequence *amp_tidal, const REAL8Sequence *fHz, REAL8 m1_SI, REAL8 m2_SI, REAL8 lambda1, REAL8 lambda2)
Function to call amplitude tidal series only; done for convenience to use for PhenomD_NRTidalv2 and S...
static const REAL8 Mf_ROM_max
static const REAL8 Mf_ROM_min
static void nudge(REAL8 *x, REAL8 X, REAL8 epsilon)
Auxiliary functions for SEOBNRv1/v2 reduced order modeling codes in LALSimIMRSEOBNRv2ROMDoubleSpin....
static UNUSED gsl_vector * Fit_cubic(const gsl_vector *xi, const gsl_vector *yi)
static UNUSED REAL8 Interpolate_Coefficent_Tensor(gsl_vector *v, REAL8 eta, REAL8 chi1, REAL8 chi2, int ncy, int ncz, gsl_bspline_workspace *bwx, gsl_bspline_workspace *bwy, gsl_bspline_workspace *bwz)
static UNUSED double SEOBNRROM_Ringdown_Mf_From_Mtot_Eta(const double Mtot_sec, const double eta, const double chi1, const double chi2, Approximant apx)
static UNUSED int SEOBNRROMdataDS_Init(SEOBNRROMdataDS *romdata, const char dir[])
int XLALSimIMRSEOBNRv4ROMFrequencyOfTime(REAL8 *frequency, REAL8 t, REAL8 m1SI, REAL8 m2SI, REAL8 chi1, REAL8 chi2)
Compute the starting frequency so that the given amount of 'time' elapses in the ROM waveform from th...
static UNUSED void GluePhasing(SEOBNRROMdataDS_submodel *submodel_lo, SEOBNRROMdataDS_submodel *submodel_hi, gsl_vector *phi_f_lo, gsl_vector *phi_f_hi, const double Mfm, gsl_interp_accel **acc_phi_out, gsl_spline **spline_phi_out)
static UNUSED void SEOBNRv4ROM_Init_LALDATA(void)
Setup SEOBNRv4ROM model using data files installed in $LAL_DATA_PATH.
static UNUSED int SEOBNRROMdataDS_Init_submodel(UNUSED SEOBNRROMdataDS_submodel **submodel, UNUSED const char dir[], UNUSED const char grp_name[])
int(* load_dataPtr)(const char *, gsl_vector *, gsl_vector *, gsl_matrix *, gsl_matrix *, gsl_vector *)
static UNUSED int SEOBNRv4ROM_Init(const char dir[])
Setup SEOBNRv4ROM model using data files installed in dir.
static UNUSED void SplineData_Destroy(SplineData *splinedata)
static UNUSED REAL8 Interpolate_Coefficent_Matrix(gsl_vector *v, REAL8 eta, REAL8 chi, int ncx, int ncy, gsl_bspline_workspace *bwx, gsl_bspline_workspace *bwy)
static UNUSED int SEOBNRv4ROMTimeFrequencySetup(gsl_spline **spline_phi, gsl_interp_accel **acc_phi, REAL8 *Mf_final, REAL8 *Mtot_sec, REAL8 m1SI, REAL8 m2SI, REAL8 chi1, REAL8 chi2, REAL8 *Mf_ROM_min, REAL8 *Mf_ROM_max)
static UNUSED void SplineData_Init(SplineData **splinedata, int ncx, int ncy, int ncz, const double *etavec, const double *chi1vec, const double *chi2vec)
static SEOBNRROMdataDS __lalsim_SEOBNRv4ROMDS_data
static UNUSED int SEOBNRv4ROMCore(COMPLEX16FrequencySeries **hptilde, COMPLEX16FrequencySeries **hctilde, double phiRef, double fRef, double distance, double inclination, double Mtot_sec, double eta, double chi1, double chi2, const REAL8Sequence *freqs, double deltaF, int nk_max, LALDict *LALparams, NRTidal_version_type NRTidal_version)
Core function for computing the ROM waveform.
static UNUSED void SEOBNRROMdataDS_coeff_Init(SEOBNRROMdataDS_coeff **romdatacoeff, int nk_amp, int nk_phi)
static UNUSED void SEOBNRROMdataDS_Cleanup_submodel(SEOBNRROMdataDS_submodel *submodel)
static UNUSED void SEOBNRROMdataDS_coeff_Cleanup(SEOBNRROMdataDS_coeff *romdatacoeff)
static UNUSED void SEOBNRROMdataDS_Cleanup(SEOBNRROMdataDS *romdata)
static UNUSED void GlueAmplitude(SEOBNRROMdataDS_submodel *submodel_lo, SEOBNRROMdataDS_submodel *submodel_hi, gsl_vector *amp_f_lo, gsl_vector *amp_f_hi, double amp_pre_lo, double amp_pre_hi, const double Mfm, gsl_interp_accel **acc_amp, gsl_spline **spline_amp)
static int TP_Spline_interpolation_3d(REAL8 eta, REAL8 chi1, REAL8 chi2, gsl_vector *cvec_amp, gsl_vector *cvec_phi, int nk_amp, int nk_phi, int nk_max, int ncx, int ncy, int ncz, const double *etavec, const double *chi1vec, const double *chi2vec, gsl_vector *c_amp, gsl_vector *c_phi)
static size_t NextPow2(const size_t n)
static UNUSED bool SEOBNRv4ROM_IsSetup(void)
Helper function to check if the SEOBNRv4ROM model has been initialised.
int XLALSimIMRSEOBNRv4ROMTimeOfFrequency(REAL8 *t, REAL8 frequency, REAL8 m1SI, REAL8 m2SI, REAL8 chi1, REAL8 chi2)
Compute the 'time' elapsed in the ROM waveform from a given starting frequency until the ringdown.
#define XLAL_FILE_RESOLVE_PATH(fname)
COMPLEX16FrequencySeries * XLALCreateCOMPLEX16FrequencySeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaF, const LALUnit *sampleUnits, size_t length)
void XLALH5FileClose(LALH5File UNUSED *file)
LALH5File * XLALH5GroupOpen(LALH5File UNUSED *file, const char UNUSED *name)
LALH5File * XLALH5FileOpen(const char UNUSED *path, const char UNUSED *mode)
void * XLALCalloc(size_t m, size_t n)
void * XLALMalloc(size_t n)
@ NRTidalv2_V
version NRTidalv2: https://arxiv.org/abs/1905.06011
int XLALSimIMRSEOBNRv4ROM(struct tagCOMPLEX16FrequencySeries **hptilde, struct tagCOMPLEX16FrequencySeries **hctilde, REAL8 phiRef, REAL8 deltaF, REAL8 fLow, REAL8 fHigh, REAL8 fRef, REAL8 distance, REAL8 inclination, REAL8 m1SI, REAL8 m2SI, REAL8 chi1, REAL8 chi2, INT4 nk_max, LALDict *LALparams, NRTidal_version_type NRTidal_version)
Compute waveform in LAL format for the SEOBNRv4_ROM model.
int XLALSimIMRSEOBNRv4ROMFrequencySequence(struct tagCOMPLEX16FrequencySeries **hptilde, struct tagCOMPLEX16FrequencySeries **hctilde, const REAL8Sequence *freqs, REAL8 phiRef, REAL8 fRef, REAL8 distance, REAL8 inclination, REAL8 m1SI, REAL8 m2SI, REAL8 chi1, REAL8 chi2, INT4 nk_max, LALDict *LALparams, NRTidal_version_type NRTidal_version)
Compute waveform in LAL format at specified frequencies for the SEOBNRv4_ROM model.
@ SEOBNRv4
Spin nonprecessing EOBNR model v4.
void XLALDestroyREAL8Sequence(REAL8Sequence *sequence)
REAL8Sequence * XLALCreateREAL8Sequence(size_t length)
const LALUnit lalStrainUnit
const LALUnit lalSecondUnit
LALUnit * XLALUnitMultiply(LALUnit *output, const LALUnit *unit1, const LALUnit *unit2)
#define XLAL_ERROR_VOID(...)
int int int XLALPrintInfo(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_CHECK(assertion,...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
int int XLALPrintWarning(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
int pthread_once(pthread_once_t *once_control, void(*init_routine)(void))
gsl_bspline_workspace * bwx
gsl_bspline_workspace * bwz
gsl_bspline_workspace * bwy
SEOBNRROMdataDS_submodel * sub2
SEOBNRROMdataDS_submodel * sub1
SEOBNRROMdataDS_submodel * sub3