39#include <lal/XLALError.h>
41#include <gsl/gsl_math.h>
42#include <gsl/gsl_multifit.h>
43#include <gsl/gsl_interp.h>
44#include <gsl/gsl_fit.h>
48#ifdef LAL_HDF5_ENABLED
49#include <lal/H5FileIO.h>
52UNUSED
static int read_vector(
const char dir[],
const char fname[], gsl_vector *v);
53UNUSED
static int read_matrix(
const char dir[],
const char fname[], gsl_matrix *
m);
58UNUSED
static UINT4 blend_functions(gsl_vector * freqs_out, gsl_vector * out_fun, gsl_vector * freq_in_1, gsl_vector * fun_in_1, gsl_vector * freq_in_2, gsl_vector * fun_in_2,
REAL8 freq_1,
REAL8 freq_2);
67#ifdef LAL_HDF5_ENABLED
68UNUSED
static int CheckVectorFromHDF5(
LALH5File *
file,
const char name[],
const double *v,
size_t n);
71UNUSED
static int ReadHDF5LongVectorDataset(
LALH5File *
file,
const char *
name, gsl_vector_long **
data);
72UNUSED
static int ReadHDF5LongMatrixDataset(
LALH5File *
file,
const char *
name, gsl_matrix_long **
data);
73UNUSED
static void PrintInfoStringAttribute(
LALH5File *
file,
const char attribute[]);
83 gsl_bspline_workspace *bwx,
84 gsl_bspline_workspace *bwy,
85 gsl_bspline_workspace *bwz
94 gsl_bspline_workspace *bwx,
95 gsl_bspline_workspace *bwy
98UNUSED
static gsl_vector *
Fit_cubic(
const gsl_vector *xi,
const gsl_vector *yi);
104 const double Mtot_sec,
111 const double Mtot_sec,
119static int read_vector(
const char dir[],
const char fname[], gsl_vector *v) {
120 size_t size = strlen(dir) + strlen(fname) + 2;
122 snprintf(
path, size,
"%s/%s", dir, fname);
124 FILE *f = fopen(
path,
"rb");
127 int ret = gsl_vector_fread(f, v);
137static int read_matrix(
const char dir[],
const char fname[], gsl_matrix *
m) {
138 size_t size = strlen(dir) + strlen(fname) + 2;
140 snprintf(
path, size,
"%s/%s", dir, fname);
142 FILE *f = fopen(
path,
"rb");
145 int ret = gsl_matrix_fread(f,
m);
155#ifdef LAL_HDF5_ENABLED
156static int CheckVectorFromHDF5(
LALH5File *
file,
const char name[],
const double *v,
size_t n) {
157 gsl_vector *temp = NULL;
158 ReadHDF5RealVectorDataset(
file,
name, &temp);
163 for (
size_t i=0;
i<temp->size;
i++)
164 if (gsl_vector_get(temp,
i) != v[
i])
166 gsl_vector_free(temp);
188 if (dimLength == NULL) {
192 if (dimLength->
length != 1) {
197 n = dimLength->
data[0];
201 *
data = gsl_vector_alloc(n);
207 else if ((*data)->size != n) {
240 if (dimLength == NULL) {
244 if (dimLength->
length != 2) {
249 n1 = dimLength->
data[0];
250 n2 = dimLength->
data[1];
254 *
data = gsl_matrix_alloc(n1, n2);
260 else if ((*data)->size1 != n1 || (*data)->size2 != n2) {
293 if (dimLength == NULL) {
297 if (dimLength->
length != 1) {
302 n = dimLength->
data[0];
306 *
data = gsl_vector_long_alloc(n);
312 else if ((*data)->size != n) {
345 if (dimLength == NULL) {
349 if (dimLength->
length != 2) {
354 n1 = dimLength->
data[0];
355 n2 = dimLength->
data[1];
359 *
data = gsl_matrix_long_alloc(n1, n2);
365 else if ((*data)->size1 != n1 || (*data)->size2 != n2) {
380static void PrintInfoStringAttribute(
LALH5File *
file,
const char attribute[]) {
400 gsl_bspline_workspace *bwx,
401 gsl_bspline_workspace *bwy,
402 gsl_bspline_workspace *bwz
405 gsl_vector *Bx4 = gsl_vector_alloc(4);
406 gsl_vector *By4 = gsl_vector_alloc(4);
407 gsl_vector *Bz4 = gsl_vector_alloc(4);
409 size_t isx, isy, isz;
426 for (
int i=0;
i<4;
i++)
427 for (
int j=0; j<4; j++)
428 for (
int k=0; k<4; k++) {
432 double cijk = gsl_vector_get(v, (ii*ncy + jj)*ncz + kk);
433 sum += cijk * gsl_vector_get(Bx4,
i) * gsl_vector_get(By4, j) * gsl_vector_get(Bz4, k);
436 gsl_vector_free(Bx4);
437 gsl_vector_free(By4);
438 gsl_vector_free(Bz4);
452 gsl_bspline_workspace *bwx,
453 gsl_bspline_workspace *bwy
455 gsl_matrix
c = gsl_matrix_view_vector(v, ncx, ncy).matrix;
458 gsl_vector *Bx4 = gsl_vector_alloc(4);
459 gsl_vector *By4 = gsl_vector_alloc(4);
472 for (
int i=0;
i<4;
i++)
473 for (
int j=0; j<4; j++)
474 sum += gsl_matrix_get(&
c, isx +
i, isy + j) * gsl_vector_get(Bx4,
i) * gsl_vector_get(By4, j);
476 gsl_vector_free(Bx4);
477 gsl_vector_free(By4);
483static gsl_vector *
Fit_cubic(
const gsl_vector *xi,
const gsl_vector *yi) {
484 const int n = xi->size;
487 gsl_matrix *X = gsl_matrix_alloc(n,
p);
488 gsl_vector *
y = gsl_vector_alloc(n);
489 gsl_vector *
c = gsl_vector_alloc(
p);
490 gsl_matrix *cov = gsl_matrix_alloc(
p,
p);
493 for (
int i=0;
i < n;
i++) {
494 double xval = gsl_vector_get(xi,
i);
495 double xval2 = xval*xval;
496 double xval3 = xval2*xval;
497 gsl_matrix_set(X,
i, 0, 1.0);
498 gsl_matrix_set(X,
i, 1, xval);
499 gsl_matrix_set(X,
i, 2, xval2);
500 gsl_matrix_set(X,
i, 3, xval3);
502 double yval = gsl_vector_get(yi,
i);
503 gsl_vector_set (
y,
i, yval);
506 gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(n,
p);
507 gsl_multifit_linear(X,
y,
c, cov, &chisq, work);
510 gsl_multifit_linear_free(work);
512 gsl_matrix_free(cov);
521 return !gsl_fcmp(
x,
y, epsilon);
534 if (fabs(*
x - X) < epsilon)
541 const double Mtot_sec,
549 double m1_SI = Mtot_SI *
q / (1. +
q);
550 double m2_SI = Mtot_SI * 1. / (1. +
q);
556 const double Mtot_sec,
563 double q = (1. + sqrt(1. - 4. * eta) - 2. * eta) / (2. * eta);
579 for(
i=0 ;
i<count ;
i++)
580 len += strlen(va_arg(ap,
char*));
589 for(
i=0 ;
i<count ;
i++)
591 char *
s = va_arg(ap,
char*);
592 strcpy(merged+null_pos,
s);
593 null_pos += strlen(
s);
610 return_value = 1 / (1 + exp_value);
617 for(
unsigned int i = 0;
i < freqs->size;
i++){
618 if(freqs->data[
i] <= freq_1){
619 out_fun->data[
i] = 0.;
621 else if((freqs->data[
i]>freq_1) && (freqs->data[
i]<freq_2)){
622 out_fun->data[
i] =
sigmoid(- (freq_2 - freq_1)/(freqs->data[
i]-freq_1) - (freq_2 - freq_1)/(freqs->data[
i]-freq_2));
625 out_fun->data[
i] = 1.;
634 for(
unsigned int i = 0;
i < freqs_in_1->size;
i++){
635 if(freqs_in_1->data[
i]<freq_1){
639 for(
unsigned int i = freqs_in_2->size -1;
i > 0;
i--){
640 if(freqs_in_2->data[
i]>=freq_1){
650 REAL8 diff = fabs(freq-freq_array->data[index]);
651 for(
unsigned int i = 1;
i < freq_array->size;
i++){
652 if(fabs(freq-freq_array->data[
i])< diff){
653 diff = fabs(freq-freq_array->data[
i]);
661UINT4 blend_functions(gsl_vector * freqs_out, gsl_vector * out_fun, gsl_vector * freq_in_1, gsl_vector * fun_in_1, gsl_vector * freq_in_2, gsl_vector * fun_in_2,
REAL8 freq_1,
REAL8 freq_2){
663 gsl_vector *bld_2 = gsl_vector_calloc(out_fun->size);
666 UNUSED
UINT4 ret =
blend(freqs_out, bld_2,freq_1,freq_2);
668 gsl_interp_accel *acc = gsl_interp_accel_alloc ();
669 gsl_spline *spline = gsl_spline_alloc (gsl_interp_cspline, freq_in_1->size);
670 gsl_spline_init (spline, freq_in_1->data, fun_in_1->data, freq_in_1->size);
679 while(freqs_out->data[
i]<=freq_2){
681 inter_func = gsl_spline_eval (spline, freqs_out->data[
i], acc);
683 out_fun->data[
i] = inter_func*(1.-bld_2->data[
i]);
686 gsl_spline_free (spline);
687 gsl_interp_accel_free(acc);
692 spline = gsl_spline_alloc (gsl_interp_cspline, freq_in_2->size);
693 acc = gsl_interp_accel_alloc ();
694 gsl_spline_init (spline, freq_in_2->data, fun_in_2->data, freq_in_2->size);
695 i = freqs_out->size-1;
696 while(freqs_out->data[
i]>freq_2){
697 inter_func = gsl_spline_eval (spline, freqs_out->data[
i], acc);
698 out_fun->data[
i] = inter_func;
704 while((freqs_out->data[
i]>=freq_1)&&(freqs_out->data[
i]<=freq_2)){
705 inter_func = gsl_spline_eval (spline, freqs_out->data[
i], acc);
706 out_fun->data[
i] += inter_func*bld_2->data[
i];
712 gsl_vector_free(bld_2);
713 gsl_spline_free (spline);
714 gsl_interp_accel_free(acc);
726 if((f_align_start<f_array_1->
data[0]) || (f_align_start<f_array_2->
data[0]) || (f_align_end>f_array_1->data[f_array_1->size-1]) || (f_align_end>f_array_2->data[f_array_2->size-1])) {
733 gsl_interp_accel *acc_1 = gsl_interp_accel_alloc();
734 gsl_spline *spline_1 = gsl_spline_alloc(gsl_interp_cspline, phase_1->size);
735 gsl_spline_init(spline_1, f_array_1->data, phase_1->data, phase_1->size);
736 gsl_interp_accel *acc_2 = gsl_interp_accel_alloc();
737 gsl_spline *spline_2 = gsl_spline_alloc(gsl_interp_cspline, phase_2->size);
738 gsl_spline_init(spline_2, f_array_2->data, phase_2->data, phase_2->size);
739 gsl_vector* delta_phi = gsl_vector_alloc(npt_fit);
740 gsl_vector* freq_delta_phi = gsl_vector_alloc(npt_fit);
743 for(
unsigned int i=0;
i<freq_delta_phi->size;
i++){
744 f = f_align_start + ((double)
i)/(freq_delta_phi->size - 1) * (f_align_end-f_align_start);
745 freq_delta_phi->data[
i] = f;
746 delta_phi->data[
i] = gsl_spline_eval(spline_2, f, acc_2) - gsl_spline_eval(spline_1, f, acc_1);
748 REAL8 c0,
c1, cov00, cov01, cov11, chisq;
749 gsl_fit_linear(freq_delta_phi->data, 1, delta_phi->data, 1, freq_delta_phi->size, &
c0, &
c1, &cov00, &cov01, &cov11, &chisq);
754 for(
unsigned int i = 0;
i < f_array_1->size;
i++){
755 phase_1->data[
i] = phase_1->data[
i] + (
c1*f_array_1->data[
i] +
c0);
758 gsl_vector_free(delta_phi);
759 gsl_vector_free(freq_delta_phi);
760 gsl_spline_free(spline_1);
761 gsl_interp_accel_free(acc_1);
762 gsl_spline_free(spline_2);
763 gsl_interp_accel_free(acc_2);
778 gsl_vector_fprintf(
fp, v,
"%g");
784 for (
size_t i=0;
i <
s->size;
i++) {
792 const char filename[],
796 gsl_vector *v = gsl_vector_alloc(mode->
data->
length);
797 double (*part_fun_ptr)(
double complex);
799 part_fun_ptr = creal;
801 part_fun_ptr = cimag;
803 gsl_vector_set(v,
i, part_fun_ptr(mode->
data->
data[
i]));
816 if((f_align_start<f_array_1->
data[0]) || (f_align_start<f_array_2->
data[0]) || (f_align_end>f_array_1->data[f_array_1->size-1]) || (f_align_end>f_array_2->data[f_array_2->size-1])) {
823 gsl_interp_accel *acc_1 = gsl_interp_accel_alloc();
824 gsl_spline *spline_1 = gsl_spline_alloc(gsl_interp_cspline, phase_1->size);
825 gsl_spline_init(spline_1, f_array_1->data, phase_1->data, phase_1->size);
826 gsl_interp_accel *acc_2 = gsl_interp_accel_alloc();
827 gsl_spline *spline_2 = gsl_spline_alloc(gsl_interp_cspline, phase_2->size);
828 gsl_spline_init(spline_2, f_array_2->data, phase_2->data, phase_2->size);
829 gsl_vector* delta_phi = gsl_vector_alloc(npt_fit);
830 gsl_vector* freq_delta_phi = gsl_vector_alloc(npt_fit);
833 for(
unsigned int i=0;
i<freq_delta_phi->size;
i++){
834 f = f_align_start + ((double)
i)/(freq_delta_phi->size - 1) * (f_align_end-f_align_start);
835 delta_phi->data[
i] = gsl_spline_eval(spline_2, f, acc_2) - gsl_spline_eval(spline_1, f, acc_1) - (2*
LAL_PI*Deltat22*f + modeM/2.*Deltaphi22);
839 REAL8 av_phase_diff = 0.;
840 for(
unsigned int i = 0;
i < freq_delta_phi->size;
i++){
841 av_phase_diff += delta_phi->data[
i];
843 av_phase_diff = av_phase_diff / (delta_phi->size);
849 for(
unsigned int i = 0;
i < f_array_1->size;
i++){
850 phase_1->data[
i] = phase_1->data[
i] + (2*
LAL_PI*Deltat22*f_array_1->data[
i] + modeM/2.*Deltaphi22 + pishift);
854 gsl_vector_free(delta_phi);
855 gsl_vector_free(freq_delta_phi);
856 gsl_spline_free(spline_1);
857 gsl_interp_accel_free(acc_1);
858 gsl_spline_free(spline_2);
859 gsl_interp_accel_free(acc_2);
876 modefreqVec.
data = &modeFreq;
877 REAL8 spin1[3] = {0., 0., chi1z};
878 REAL8 spin2[3] = {0., 0., chi2z};
881 return Ms * creal(modefreqVec.
data[0]);
895 modefreqVec.
data = &modeFreq;
896 REAL8 spin1[3] = {0., 0., chi1z};
897 REAL8 spin2[3] = {0., 0., chi2z};
900 return Ms * creal(modefreqVec.
data[0]);
906 gsl_vector* phaseout,
909 int N = phasein->size;
910 double*
p = phasein->data;
911 double* pmod = (
double*) malloc(
sizeof(
double) *
N);
912 int* jumps = (
int*) malloc(
sizeof(
int) *
N);
913 int* cumul = (
int*) malloc(
sizeof(
int) *
N);
916 for(
int i=0;
i<
N;
i++) {
923 for(
int i=1;
i<
N;
i++) {
924 d = pmod[
i] - pmod[
i-1];
933 for(
int i=1;
i<
N;
i++) {
939 double* pout = phaseout->data;
940 for(
int i=0;
i<
N;
i++) {
struct tagLALH5File LALH5File
struct tagLALH5Dataset LALH5Dataset
INT4 XLALSimIMREOBGenerateQNMFreqV5(COMPLEX16Vector *modefreqs, const REAL8 mass1, const REAL8 mass2, const REAL8 spin1[3], const REAL8 spin2[3], UINT4 l, INT4 m, UINT4 nmodes, Approximant approximant)
These functions generate the quasinormal mode frequencies for a black hole ringdown,...
static UNUSED void test_save_gsl_vector(const char filename[], gsl_vector *v)
static UNUSED UINT8 compute_i_at_f(gsl_vector *freq_array, REAL8 freq)
static UNUSED gsl_vector * Fit_cubic(const gsl_vector *xi, const gsl_vector *yi)
static UNUSED void nudge(REAL8 *x, REAL8 X, REAL8 epsilon)
static UNUSED UINT4 align_wfs_window_from_22(gsl_vector *f_array_1, gsl_vector *f_array_2, gsl_vector *phase_1, gsl_vector *phase_2, REAL8 f_align_start, REAL8 f_align_end, REAL8 Deltat22, REAL8 Deltaphi22, INT4 modeM)
static UNUSED REAL8 Get_omegaQNM_SEOBNRv5(REAL8 q, REAL8 chi1z, REAL8 chi2z, UINT4 l, UINT4 m)
static UNUSED double SEOBNRROM_Ringdown_Mf_From_Mtot_q(const double Mtot_sec, const double q, const double chi1, const double chi2, Approximant apx)
static UNUSED void test_save_gsl_spline(const char filename[], gsl_spline *s)
static UNUSED UINT4 compute_i_max_LF_i_min_HF(INT8 *i_max_LF, INT8 *i_min_LF, gsl_vector *freqs_in_1, gsl_vector *freqs_in_2, REAL8 freq_1)
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 REAL8 Get_omegaQNM_SEOBNRv4(REAL8 q, REAL8 chi1z, REAL8 chi2z, UINT4 l, UINT4 m)
static UNUSED UINT4 unwrap_phase(gsl_vector *phaseout, gsl_vector *phasein)
static UNUSED int read_vector(const char dir[], const char fname[], gsl_vector *v)
static UNUSED void test_save_cmplx_freq_series(COMPLEX16FrequencySeries *mode, const char filename[], bool save_real_part)
static UNUSED UINT4 align_wfs_window(gsl_vector *f_array_1, gsl_vector *f_array_2, gsl_vector *phase_1, gsl_vector *phase_2, REAL8 *Deltat, REAL8 *Deltaphi, REAL8 f_align_start, REAL8 f_align_end)
static UNUSED UINT4 blend(gsl_vector *freqs, gsl_vector *out_fun, REAL8 freq_1, REAL8 freq_2)
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 char * concatenate_strings(int count,...)
static UNUSED REAL8 sigmoid(REAL8 x)
static UNUSED UINT4 blend_functions(gsl_vector *freqs_out, gsl_vector *out_fun, gsl_vector *freq_in_1, gsl_vector *fun_in_1, gsl_vector *freq_in_2, gsl_vector *fun_in_2, REAL8 freq_1, REAL8 freq_2)
static UNUSED int read_matrix(const char dir[], const char fname[], gsl_matrix *m)
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 bool approximately_equal(REAL8 x, REAL8 y, REAL8 epsilon)
double XLALSimInspiralGetFinalFreq(REAL8 m1, REAL8 m2, const REAL8 S1x, const REAL8 S1y, const REAL8 S1z, const REAL8 S2x, const REAL8 S2y, const REAL8 S2z, Approximant approximant)
Function that gives the default ending frequencies of the given approximant.
static vector d(const double L_norm, const double J_norm, const vector roots)
Internal function that returns the coefficients "d_0", "d_2" and "d_4" from 1703.03967 corresponding ...
static REAL8 sum(REAL8 *arr1, REAL8 *arr2, int k)
void XLALH5DatasetFree(LALH5Dataset UNUSED *dset)
int XLALH5AttributeQueryStringValue(char UNUSED *value, size_t UNUSED size, const LALH5Generic UNUSED object, const char UNUSED *key)
LALTYPECODE XLALH5DatasetQueryType(LALH5Dataset UNUSED *dset)
LALH5Dataset * XLALH5DatasetRead(LALH5File UNUSED *file, const char UNUSED *name)
int XLALH5DatasetQueryData(void UNUSED *data, LALH5Dataset UNUSED *dset)
int XLALH5FileQueryStringAttributeValue(char UNUSED *value, size_t UNUSED size, LALH5File UNUSED *file, const char UNUSED *key)
UINT4Vector * XLALH5DatasetQueryDims(LALH5Dataset UNUSED *dset)
void * XLALMalloc(size_t n)
void * XLALCalloc(size_t m, size_t n)
INT4 XLALSimIMREOBGenerateQNMFreqV2(COMPLEX16Vector *modefreqs, const REAL8 mass1, const REAL8 mass2, const REAL8 spin1[3], const REAL8 spin2[3], UINT4 l, INT4 m, UINT4 nmodes, Approximant approximant)
These functions generate the quasinormal mode frequencies for a black hole ringdown.
Approximant
Enum that specifies the PN approximant to be used in computing the waveform.
@ SEOBNRv4
Spin nonprecessing EOBNR model v4.
@ SEOBNRv5_ROM
Time domain, precessing phenomenological IMR waveform model with subdominant modes ([arXiv: 20XY....
void XLALDestroyUINT4Vector(UINT4Vector *vector)
#define XLAL_PRINT_INFO(...)
int int int XLALPrintInfo(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
static int gsl_bspline_basis(const double x, gsl_vector *Bk, size_t *istart, gsl_bspline_workspace *w)