Commit 4af4488d authored by Maria Haney's avatar Maria Haney
Browse files

Merge branch 'imrphenomxphm-review' into 'master'

Add IMRPhenomXP and IMRPhenomXPHM reviewed models

See merge request !1272
parents 8c2968b4 5fedeaff
Pipeline #115960 failed with stages
in 138 minutes and 24 seconds
......@@ -82,6 +82,7 @@ Bruce Allen <bruce.allen@ligo.org>
Bruce Allen <bruce.allen@ligo.org> <ballen>
Carl-Johan Haster <carl-johan.haster@ligo.org>
Carl-Johan Haster <carl-johan.haster@ligo.org> <cjhaster@star.sr.bham.ac.uk>
Cecilio Garcia <cecilio.garcia-quiros@ligo.org>
Chad Hanna <chad.hanna@ligo.org>
Chad Hanna <chad.hanna@ligo.org> <channa>
Chad Hanna <chad.hanna@ligo.org> <crh184@psu.edu>
......@@ -161,6 +162,7 @@ Frederique Marion <frederique.marion@ligo.org>
Frederique Marion <frederique.marion@ligo.org> <marionf>
Gareth Jones <gareth>
George Birthisel <gbb>
Geraint Pratten <geraint.pratten@ligo.org>
Gergely Debreczeni <Gergely.Debreczeni@cern.ch>
Gregory Mendell <gregory.mendell@ligo.org>
Gregory Mendell <gregory.mendell@ligo.org> <gmendell>
......
......@@ -757,7 +757,14 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
(--grtest-parameters dchi0,..,dxi1,..,dalpha1,..) template will assume deformations in the corresponding phase coefficients.\n\
(--ppe-parameters aPPE1,.... template will assume the presence of an arbitrary number of PPE parameters. They must be paired correctly.\n\
(--modeList lm,l-m...,lm,l-m) List of modes to be used by the model. The chosen modes ('lm') should be passed as a ',' seperated list.\n\
(--phenomXHMMband float) Threshold parameter for the Multibanding of IMRPhenomXHM. By default set to 10^-3. If set to 0 then do not use multibanding.\n\
(--phenomXHMMband float) Threshold parameter for the Multibanding of the non-precessing hlm modes in IMRPhenomXHM and IMRPhenomXPHM. If set to 0 then do not use multibanding.\n\
Options and default values can be found in https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom_x__c.html\n\
(--phenomXPHMMband float) Threshold parameter for the Multibanding of the Euler angles in IMRPhenomXPHM. If set to 0 then do not use multibanding.\n\
Options and default values can be found in https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom_x__c.html\n\
(--phenomXPFinalSpinMod int) Change version for the final spin model used in IMRPhenomXP/IMRPhenomXPHM.\n\
Options and default values can be found in https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom_x__c.html\n\
(--phenomXPrecVersion int) Change version of the Euler angles for the twisting-up of IMRPhenomXP/IMRPhenomXPHM.\n\
Options and default values can be found in https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom_x__c.html\n\
\n\
----------------------------------------------\n\
--- Starting Parameters ----------------------\n\
......@@ -1526,6 +1533,29 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
fprintf(stdout,"Template will use phenomXHMMband %s.\n",ppt->value);
}
/* Pass threshold for PhenomXPHM Multibanding of the Euler angles. */
if((ppt=LALInferenceGetProcParamVal(commandLine,"--phenomXPHMMband"))) {
REAL8 threshold = atof(ppt->value);
XLALSimInspiralWaveformParamsInsertPhenomXPHMThresholdMband(model->LALpars, threshold);
fprintf(stdout,"Template will use phenomXPHMMband %s.\n",ppt->value);
}
/* Version for the final spin model to use in IMRPhenomXP/IMRPhenomXPHM. */
if((ppt=LALInferenceGetProcParamVal(commandLine,"--phenomXPFinalSpinMod"))) {
REAL8 afversion = atoi(ppt->value);
XLALSimInspiralWaveformParamsInsertPhenomXPFinalSpinMod(model->LALpars, afversion);
fprintf(stdout,"Template will use PhenomXPFinalSpinMod %s.\n",ppt->value);
}
/* Version of the Euler angles in the twisting-up of IMRPhenomXP/IMRPhenomXPHM. */
if((ppt=LALInferenceGetProcParamVal(commandLine,"--phenomXPrecVersion"))) {
REAL8 eulerangleversion = atoi(ppt->value);
XLALSimInspiralWaveformParamsInsertPhenomXPrecVersion(model->LALpars, eulerangleversion);
fprintf(stdout,"Template will use PhenomXPrecVersion %s.\n",ppt->value);
}
fprintf(stdout,"\n\n---\t\t ---\n");
LALInferenceInitSpinVariables(state, model);
......
......@@ -317,8 +317,22 @@ amporder=0
# Control which modes to use in the waveform (default: all modes available)
#modeList=22,2-1,33,44
# Control the threshold for the multibanding of IMRPhenomXHM (by default is 10^-3). If set to 0 multibanding is switched off.
#phenomXHMMband=0
# Control the threshold for the multibanding of the non-precessing hlms in IMRPhenomXHM and IMRPhenomXPHM. If set to 0 multibanding is switched off.
# Options and default values can be found in https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom_x__c.html.
#phenomXHMMband=0.001
# Control the threshold for the multibanding of the Euler angles in IMRPhenomXPHM. If set to 0 multibanding is switched off.
# Options and default values can be found in https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom_x__c.html.
#phenomXPHMMband=0.001
# Control the version of the Euler angles for the twisting-up in IMRPhenomXP/IMRPhenomXPHM.
# Options and default values can be found in https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom_x__c.html.
#phenomXPrecVersion=223
# Change version for the final spin model used in IMRPhenomXP/IMRPhenomXPHM.
# Options and default values can be found in https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom_x__c.html.
#phenomXPFinalSpinMod=3
# maxmcmc set the maximum chain length for the nested sampling sub-chains. Default 5000
# Auto determination is on, but the length cannot be longer than that.
......
......@@ -67,7 +67,7 @@ extern "C" {
typedef enum tagIMRPhenomP_version_type {
IMRPhenomPv1_V, /**< version 1: based on IMRPhenomC */
IMRPhenomPv2_V, /**< version 2: based on IMRPhenomD */
IMRPhenomPv2NRTidal_V, /**< version Pv2_NRTidal: based on IMRPhenomPv2; NRTides added before precession; can be used with both NRTidal versions defined below */
IMRPhenomPv2NRTidal_V, /**< version Pv2_NRTidal: based on IMRPhenomPv2; NRTides added before precession; can be used with both NRTidal versions defined below */
IMRPhenomPv3_V /**< version 3: based on IMRPhenomD and the precession angles from Katerina Chatziioannou PhysRevD.95.104004 (arxiv:1703.03967) */
} IMRPhenomP_version_type;
......@@ -427,7 +427,80 @@ int XLALSimIMRPhenomXASFrequencySequence(
LALDict *lalParams
);
/* in module LALSimIMRPhenomXHM.c */
int XLALSimIMRPhenomXPMSAAngles(
REAL8Sequence *phiz_of_f, /**< [out] */
REAL8Sequence *zeta_of_f, /**< [out] */
REAL8Sequence *costhetaL_of_f, /**< [out] */
const REAL8Sequence *freqs, /**< [out] Frequency series [Hz] */
const REAL8 m1_SI, /**< mass of companion 1 (kg) */
const REAL8 m2_SI, /**< mass of companion 2 (kg) */
const REAL8 chi1x, /**< x-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
const REAL8 chi1y, /**< y-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
const REAL8 chi1z, /**< z-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
const REAL8 chi2x, /**< x-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
const REAL8 chi2y, /**< y-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
const REAL8 chi2z, /**< z-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
const REAL8 fRef_In, /**< Reference frequency (Hz) */
LALDict *lalParams /**< LAL Dictionary struct */
);
int XLALSimIMRPhenomXPPNAngles(
REAL8Sequence *alpha_of_f, /**< [out] */
REAL8Sequence *gamma_of_f, /**< [out] */
REAL8Sequence *cosbeta_of_f, /**< [out] */
const REAL8Sequence *freqs, /**< [out] Frequency series [Hz] */
const REAL8 m1_SI, /**< mass of companion 1 (kg) */
const REAL8 m2_SI, /**< mass of companion 2 (kg) */
const REAL8 chi1x, /**< x-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
const REAL8 chi1y, /**< y-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
const REAL8 chi1z, /**< z-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
const REAL8 chi2x, /**< x-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
const REAL8 chi2y, /**< y-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
const REAL8 chi2z, /**< z-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
const REAL8 fRef_In, /**< Reference frequency (Hz) */
LALDict *lalParams /**< LAL Dictionary struct */
);
int XLALSimIMRPhenomXPGenerateFD(
COMPLEX16FrequencySeries **hptilde, /**< [out] Frequency-domain waveform h+ */
COMPLEX16FrequencySeries **hctilde, /**< [out] Frequency-domain waveform hx */
REAL8 m1_SI, /**< mass of companion 1 (kg) */
REAL8 m2_SI, /**< mass of companion 2 (kg) */
REAL8 chi1x, /**< x-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
REAL8 chi1y, /**< y-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
REAL8 chi1z, /**< z-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
REAL8 chi2x, /**< x-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
REAL8 chi2y, /**< y-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
REAL8 chi2z, /**< z-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
const REAL8 distance, /**< Distance of source (m) */
const REAL8 inclination, /**< inclination of source (rad) */
const REAL8 phiRef, /**< Orbital phase (rad) at reference frequency */
REAL8 f_min, /**< Starting GW frequency (Hz) */
REAL8 f_max, /**< Ending GW frequency (Hz); Defaults to Mf = 0.3 if no f_max is specified. */
const REAL8 deltaF, /**< Sampling frequency (Hz). To use non-uniform frequency grid, set deltaF <= 0. */
REAL8 fRef_In, /**< Reference frequency (Hz) */
LALDict *lalParams /**< LAL Dictionary struct */
);
int XLALSimIMRPhenomXPFrequencySequence(
COMPLEX16FrequencySeries **hptilde, /**< [out] Frequency-domain waveform h+ */
COMPLEX16FrequencySeries **hctilde, /**< [out] Frequency-domain waveform hx */
const REAL8Sequence *freqs, /**< [out] Frequency series [Hz] */
REAL8 m1_SI, /**< mass of companion 1 (kg) */
REAL8 m2_SI, /**< mass of companion 2 (kg) */
REAL8 chi1x, /**< x-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
REAL8 chi1y, /**< y-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
REAL8 chi1z, /**< z-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
REAL8 chi2x, /**< x-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
REAL8 chi2y, /**< y-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
REAL8 chi2z, /**< z-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
const REAL8 distance, /**< Distance of source (m) */
const REAL8 inclination, /**< inclination of source (rad) */
const REAL8 phiRef, /**< Orbital phase (rad) at reference frequency */
REAL8 fRef_In, /**< Reference frequency (Hz) */
LALDict *lalParams /**< LAL Dictionary struct */
);
int XLALSimIMRPhenomXHMGenerateFDOneMode(
COMPLEX16FrequencySeries **htildelm, /**< [out] FD waveform */
REAL8 m1_SI, /**< Mass of companion 1 (kg) */
......@@ -565,6 +638,132 @@ int XLALSimIMRPhenomXHMPhase(
LALDict *lalParams /**< Extra params */
);
int XLALSimIMRPhenomXPHM(
COMPLEX16FrequencySeries **hptilde, /**< [out] Frequency-domain waveform h+ */
COMPLEX16FrequencySeries **hctilde, /**< [out] Frequency-domain waveform hx */
REAL8 m1_SI, /**< mass of companion 1 (kg) */
REAL8 m2_SI, /**< mass of companion 2 (kg) */
REAL8 chi1x, /**< x-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
REAL8 chi1y, /**< y-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
REAL8 chi1z, /**< z-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
REAL8 chi2x, /**< x-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
REAL8 chi2y, /**< y-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
REAL8 chi2z, /**< z-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
REAL8 distance, /**< Distance of source (m) */
REAL8 inclination, /**< inclination of source (rad) */
REAL8 phiRef, /**< Orbital phase (rad) at reference frequency */
REAL8 f_min, /**< Starting GW frequency (Hz) */
REAL8 f_max, /**< Ending GW frequency (Hz); Defaults to Mf = 0.3 if no f_max is specified. */
REAL8 deltaF, /**< Sampling frequency (Hz). To use non-uniform frequency grid, set deltaF <= 0. */
REAL8 fRef_In, /**< Reference frequency (Hz) */
LALDict *lalParams /**< LAL Dictionary struct */
);
int XLALSimIMRPhenomXPHMFromModes(
COMPLEX16FrequencySeries **hptilde, /**< [out] Frequency-domain waveform h+ */
COMPLEX16FrequencySeries **hctilde, /**< [out] Frequency-domain waveform hx */
REAL8 m1_SI, /**< mass of companion 1 (kg) */
REAL8 m2_SI, /**< mass of companion 2 (kg) */
REAL8 chi1x, /**< x-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
REAL8 chi1y, /**< y-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
REAL8 chi1z, /**< z-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
REAL8 chi2x, /**< x-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
REAL8 chi2y, /**< y-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
REAL8 chi2z, /**< z-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
const REAL8 distance, /**< Distance of source (m) */
const REAL8 inclination, /**< inclination of source (rad) */
const REAL8 phiRef, /**< Orbital phase (rad) at reference frequency */
REAL8 f_min, /**< Starting GW frequency (Hz) */
REAL8 f_max, /**< Ending GW frequency (Hz); Defaults to Mf = 0.3 if no f_max is specified. */
const REAL8 deltaF, /**< Sampling frequency (Hz). To use non-uniform frequency grid, set deltaF <= 0. */
REAL8 fRef_In, /**< Reference frequency (Hz) */
LALDict *lalParams /**< LAL Dictionary struct */
);
int XLALSimIMRPhenomXHMFrequencySequence(
COMPLEX16FrequencySeries **hptilde, /**< [out] Frequency-domain waveform h+ */
COMPLEX16FrequencySeries **hctilde, /**< [out] Frequency-domain waveform hx */
REAL8Sequence *freqs, /**< Input Frequency series [Hz] */
REAL8 m1_SI, /**< mass of companion 1 (kg) */
REAL8 m2_SI, /**< mass of companion 2 (kg) */
REAL8 chi1z, /**< z-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
REAL8 chi2z, /**< z-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
REAL8 distance, /**< Distance of source (m) */
REAL8 inclination, /**< inclination of source (rad) */
REAL8 phiRef, /**< Orbital phase (rad) at reference frequency */
REAL8 fRef, /**< Reference frequency (Hz) */
LALDict *lalParams /**< LAL Dictionary struct */
);
int XLALSimIMRPhenomXPHMFrequencySequence(
COMPLEX16FrequencySeries **hptilde, /**< [out] Frequency-domain waveform h+ */
COMPLEX16FrequencySeries **hctilde, /**< [out] Frequency-domain waveform hx */
REAL8Sequence *freqs, /**< input frequency series [Hz] */
REAL8 m1_SI, /**< mass of companion 1 (kg) */
REAL8 m2_SI, /**< mass of companion 2 (kg) */
REAL8 chi1x, /**< x-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
REAL8 chi1y, /**< y-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
REAL8 chi1z, /**< z-component of the dimensionless spin of object 1 w.r.t. Lhat = (0,0,1) */
REAL8 chi2x, /**< x-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
REAL8 chi2y, /**< y-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
REAL8 chi2z, /**< z-component of the dimensionless spin of object 2 w.r.t. Lhat = (0,0,1) */
REAL8 distance, /**< Distance of source (m) */
REAL8 inclination, /**< inclination of source (rad) */
REAL8 phiRef, /**< Orbital phase (rad) at reference frequency */
REAL8 fRef_In, /**< Reference frequency (Hz) */
LALDict *lalParams /**< LAL Dictionary struct */
);
int XLALSimIMRPhenomXPHMOneMode(
COMPLEX16FrequencySeries **hlmpos, /**< [out] Frequency-domain waveform hlm inertial frame positive frequencies */
COMPLEX16FrequencySeries **hlmneg, /**< [out] Frequency-domain waveform hlm inertial frame negative frequencies */
const UINT4 l, /**< First index of the (l,m) precessing mode */
const INT4 m, /**< Second index of the (l,m) precessing mode */
const REAL8 m1_SI, /**< mass of companion 1 (kg) */
const REAL8 m2_SI, /**< mass of companion 2 (kg) */
const REAL8 chi1x, /**< x-component of the dimensionless spin of object 1 */
const REAL8 chi1y, /**< y-component of the dimensionless spin of object 1 */
const REAL8 chi1z, /**< z-component of the dimensionless spin of object 1 */
const REAL8 chi2x, /**< x-component of the dimensionless spin of object 2 */
const REAL8 chi2y, /**< y-component of the dimensionless spin of object 2 */
const REAL8 chi2z, /**< z-component of the dimensionless spin of object 2 */
const REAL8 distance, /**< distance of source (m) */
const REAL8 phiRef, /**< reference orbital phase (rad) */
const REAL8 deltaF, /**< Sampling frequency (Hz) */
const REAL8 f_min, /**< Starting GW frequency (Hz) */
const REAL8 f_max, /**< End frequency; 0 defaults to ringdown cutoff freq */
const REAL8 fRef_In, /**< Reference frequency */
LALDict *lalParams /**<LAL Dictionary */
);
/*
XLAL routine to calculate the IMRPhenomX mode parameters from the source frame.
Note that this is just an IMRPhenomX compatible implementation of:
XLALSimIMRPhenomPCalculateModelParametersFromSourceFrame
*/
int XLALSimIMRPhenomXPCalculateModelParametersFromSourceFrame(
REAL8 *chi1L, /**< [out] Dimensionless aligned spin on companion 1 */
REAL8 *chi2L, /**< [out] Dimensionless aligned spin on companion 2 */
REAL8 *chi_p, /**< [out] Effective spin in the orbital plane */
REAL8 *thetaJN, /**< [out] Angle between J0 and line of sight (z-direction) */
REAL8 *alpha0, /**< [out] Initial value of alpha angle (azimuthal precession angle) */
REAL8 *phi_aligned, /**< [out] Initial phase to feed the underlying aligned-spin model */
REAL8 *zeta_polarization, /**< [out] Angle to rotate the polarizations */
const REAL8 m1_SI, /**< Mass of companion 1 (kg) */
const REAL8 m2_SI, /**< Mass of companion 2 (kg) */
const REAL8 f_ref, /**< Reference GW frequency (Hz) */
const REAL8 phiRef, /**< Reference phase */
const REAL8 incl, /**< Inclination : angle between LN and the line of sight */
const REAL8 s1x, /**< Initial value of s1x: dimensionless spin of BH 1 */
const REAL8 s1y, /**< Initial value of s1y: dimensionless spin of BH 1 */
const REAL8 s1z, /**< Initial value of s1z: dimensionless spin of BH 1 */
const REAL8 s2x, /**< Initial value of s2x: dimensionless spin of BH 2 */
const REAL8 s2y, /**< Initial value of s2y: dimensionless spin of BH 2 */
const REAL8 s2z, /**< Initial value of s2z: dimensionless spin of BH 2 */
LALDict *lalParams /**< LAL Dictionary */
);
/* in module LALSimInspiralNRWaveforms.c */
......
This diff is collapsed.
......@@ -35,12 +35,12 @@ extern "C" {
/* Dimensionless frequency (Mf) at which define the end of the waveform */
#define f_CUT 0.3
#define MAX_ALLOWED_MASS_RATIO 5000
#define MAX_ALLOWED_ETA 0.0002
#include "LALSimIMRPhenomX_internals.h"
#include "LALSimIMRPhenomXUtilities.h"
#include "LALSimIMRPhenomX_precession.h"
/* Decleration for internal function to generate aligned-spin, 22 only IMRPhenomXAS waveform */
int IMRPhenomXASGenerateFD(
COMPLEX16FrequencySeries **htilde22,
const REAL8Sequence *freqs,
......@@ -48,6 +48,7 @@ int IMRPhenomXASGenerateFD(
LALDict *lalParams
);
int IMRPhenomXCheckForUniformFrequencies(REAL8Sequence *frequencies,REAL8 df);
#ifdef __cplusplus
......
This diff is collapsed.
......@@ -38,10 +38,6 @@ extern "C" {
#include <lal/LALSimInspiral.h>
/* CONSTANTS */
#define MAX_ALLOWED_MASS_RATIO 5000
#define MAX_ALLOWED_ETA 0.0002
// /* Returns the Fourier domain strain of just one negative mode: h_l-m. This quantity is zero for negative frequencies.
// However the positive mode h_lm is zero for positive frequencies and for the negative frequencies is equal to (-1)^l h*_l-m(-f).
// This is a wrapper function that use XLALSimIMRPhenomXASGenerateFD for the 22 mode and XLALSimIMRPhenomXHMOneMode for the higher modes. */
......@@ -58,6 +54,15 @@ int IMRPhenomXHMGenerateFDOneMode(
LALDict *lalParams /**< extra params **/
);
/* Compute the frequency array and initialize htildelm to the corresponding length. */
int SetupWFArrays(
REAL8Sequence **freqs, /**< [out] frequency grid [Hz] */
COMPLEX16FrequencySeries **htildelm, /**< [out] Frequency domain hlm GW strain */
REAL8Sequence *freqs_In, /**< fmin, fmax [Hz] */
IMRPhenomXWaveformStruct *pWF, /**< Waveform structure with parameters */
LIGOTimeGPS ligotimegps_zero /**< = {0,0} */
);
#ifdef __cplusplus
}
#endif
......
......@@ -21,10 +21,6 @@
// LALSimIMRPhenomXHM_Intermediate.c
//
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_poly.h>
#include "LALSimIMRPhenomXHM_intermediate.h"
......
......@@ -22,6 +22,8 @@
// Created by Marta on 06/02/2019.
//
#include <gsl/gsl_linalg.h>
#include "LALSimIMRPhenomXHM_structs.h"
#include "LALSimIMRPhenomXHM_qnm.h"
......
......@@ -42,8 +42,6 @@ extern "C" {
#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <lal/LALStdlib.h>
#include <lal/LALConstants.h>
......@@ -56,12 +54,6 @@ extern "C" {
#include <lal/LALSimIMR.h>
/* GSL Header Files */
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#define NMAX_INSPIRAL_COEFFICIENTS 13
// You should not declare static functions here, since this file is included in other files apart form the source one.
......
......@@ -23,6 +23,8 @@
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_spline.h>
#include "LALSimIMRPhenomXHM_multiband.h"
......@@ -500,7 +502,10 @@ int XLALSimIMRPhenomXHMMultiBandOneMode(
XLAL_CHECK(XLAL_SUCCESS == status, XLAL_EFUNC, "Error: IMRPhenomXSetWaveformVariables failed.\n");
int offset = IMRPhenomXHMMultiBandOneMode(htildelm, pWF, ell, emm, lalParams);
status = IMRPhenomXHMMultiBandOneMode(htildelm, pWF, ell, emm, lalParams);
XLAL_CHECK(status == XLAL_SUCCESS, XLAL_EFUNC, "IMRPhenomXHMMultiBandOneMode failed to generate IMRPhenomXHM waveform.");
INT4 offset = (size_t) (pWF->fMin / deltaF);
if(emmIn>0){
/* (-1)^l */
......@@ -511,7 +516,7 @@ int XLALSimIMRPhenomXHMMultiBandOneMode(
#if DEBUG == 1
printf("\nTransforming to positive m by doing (-1)^l*Conjugate, frequencies must be negatives.\n");
#endif
for(UINT4 idx=0; idx<(*htildelm)->data->length; idx++){
for(UINT4 idx=offset; idx<(*htildelm)->data->length; idx++){
(*htildelm)->data->data[idx] = minus1l*conj((*htildelm)->data->data[idx]);
}
}
......@@ -539,6 +544,9 @@ int XLALSimIMRPhenomXHMMultiBandOneMode(
return offset;
}
/** @}
@} **/
int IMRPhenomXHMMultiBandOneMode(
COMPLEX16FrequencySeries **htildelm, /**< [out] FD waveform **/
......@@ -563,18 +571,24 @@ int IMRPhenomXHMMultiBandOneMode(
XLAL_CHECK ( (iStop <= npts) && (iStart <= iStop), XLAL_EDOM,
"minimum freq index %zu and maximum freq index %zu do not fulfill 0<=ind_min<=ind_max<=htilde->data>length=%zu.", iStart, iStop, npts);
#if DEBUG == 1
printf("\n***********************\n");
printf("pWF->fMin, deltaF, iStart = %.16e %.16e %zu", pWF->fMin, deltaF, iStart);
printf("\n***********************\n");
#endif
size_t offset = iStart;
/* If it is odd mode with equal black holes then return array of zeros 0 */
if(pWF->m1_SI == pWF->m2_SI && pWF->chi1L == pWF->chi2L && emm%2!=0){ // Mode zero
if(pWF->q == 1 && pWF->chi1L == pWF->chi2L && emm%2!=0){ // Mode zero
*htildelm = XLALCreateCOMPLEX16FrequencySeries("htildelm: FD waveform", &(ligotimegps_zero), 0.0, deltaF, &lalStrainUnit, iStop);
for(unsigned int idx = 0; idx < iStop; idx++){
((*htildelm)->data->data)[idx] = 0.;
}
return offset;
return XLAL_SUCCESS;
}
/** Mode non-zero **/
......@@ -994,7 +1008,8 @@ int IMRPhenomXHMMultiBandOneMode(
#endif
/* Intialize FrequencySeries for htildelm */
size_t n = count + offset;
while(count+offset > iStop){ count--; }
size_t n = iStop;
XLAL_CHECK(XLALGPSAdd(&ligotimegps_zero, -1. / deltaF), XLAL_EFUNC, "Failed to shift the coalescence time to t=0. Tried to apply a shift of -1/df with df = %g.", deltaF);
*htildelm = XLALCreateCOMPLEX16FrequencySeries("htildelm: FD waveform", &(ligotimegps_zero), 0.0, deltaF, &lalStrainUnit, n);
......@@ -1016,6 +1031,13 @@ int IMRPhenomXHMMultiBandOneMode(
((*htildelm)->data->data)[idx + offset] = minus1l * fineAmp[idx] * expphi[idx];
}
/* Sometimes rounding error can make that count+offset < iStop, meaning that we have computed less points than those we are gonna output,
here we make sure that all the elements of htildelm are initialized and we put the extra point(s) to 0. */
for(UINT4 idx = count-1; idx < iStop-offset; idx++){
/* Reconstruct waveform: h(f) = A(f) * Exp[I phi(f)] */
((*htildelm)->data->data)[idx + offset] = 0.;
}
/* Free allocated memory */
LALFree(pAmp);
LALFree(pAmp22);
......@@ -1036,7 +1058,7 @@ int IMRPhenomXHMMultiBandOneMode(
for(UINT4 idx = 0; idx < count; idx++)
{
data = expphi[idx];
fprintf(file, "%.16f %.16e %.16e\n", idx*deltaF, creal(data), cimag(data));
fprintf(file, "%.16f %.16e %.16e\n", (idx+offset)*deltaF, creal(data), cimag(data));
}
fclose(file);
#endif
......@@ -1067,10 +1089,33 @@ int IMRPhenomXHMMultiBandOneMode(
LALFree(ILphaselm);
LALFree(IntLawpoints);
return offset;
return XLAL_SUCCESS;
}
/** @addtogroup LALSimIMRPhenomX_c
* @{
* @name Routines for IMRPhenomXHM Multibanding
* @{
* @author Cecilio García Quirós, Sascha Husa
*
* @brief C code for applying Multibanding to IMRPhenomXHM_Multimode
*
* This is a technique to make the evaluation of waveform models faster by evaluating the model
* in a coarser non-uniform grid and interpolate this to the final fine uniform grid.
* We apply this technique to the fourier domain model IMRPhenomXHM as describe in this paper: https://arxiv.org/abs/2001.10897
*
* Multibanding flags:
* ThresholdMband: Determines the strength of the Multibanding algorithm.
* The lower this value is, the slower is the evaluation but more accurate is the final waveform compared to the one without multibanding.
* - 0.001: DEFAULT value
* - 0: switch off the multibanding
*
* AmpInterpol: Determines the gsl interpolation order for the amplitude.
* - 1: linear interpolation (DEFAULT)
* - 3: cubic interpolation
*
*/
/** Returns htildelm the waveform of one mode that present mode-mixing.
The multibanding is applied to the spherical part (inspiral and intermediate) and to the spheroidal part (ringdown).
......@@ -1083,7 +1128,7 @@ int XLALSimIMRPhenomXHMMultiBandOneModeMixing(
REAL8 chi1L, /**< Dimensionless aligned spin of companion 1 */
REAL8 chi2L, /**< Dimensionless aligned spin of companion 2 */
UINT4 ell, /**< l index of the mode */
INT4 emmIn, /**< m index of the mode */
INT4 emmIn, /**< m index of the mode */
REAL8 distance, /**< Luminosity distance (m) */
REAL8 f_min, /**< Starting GW frequency (Hz) */
REAL8 f_max, /**< End frequency; 0 defaults to Mf = \ref f_CUT */
......@@ -1159,7 +1204,11 @@ int XLALSimIMRPhenomXHMMultiBandOneModeMixing(
XLAL_CHECK(XLAL_SUCCESS == status, XLAL_EFUNC, "Error: IMRPhenomXSetWaveformVariables failed.\n");
int offset = IMRPhenomXHMMultiBandOneModeMixing(htildelm, htilde22, pWF, ell, emm, lalParams);
//int offset = IMRPhenomXHMMultiBandOneModeMixing(htildelm, htilde22, pWF, ell, emm, lalParams);
status = IMRPhenomXHMMultiBandOneModeMixing(htildelm, htilde22, pWF, ell, emm, lalParams);
XLAL_CHECK(status == XLAL_SUCCESS, XLAL_EFUNC, "IMRPhenomXHMMultiBandOneModeMixing failed to generate IMRPhenomXHM waveform.");
INT4 offset = (size_t) (pWF->fMin / deltaF);
if(emmIn>0){
/* (-1)^l */
......@@ -1170,7 +1219,7 @@ int XLALSimIMRPhenomXHMMultiBandOneModeMixing(
#if DEBUG == 1
printf("\nTransforming to positive m by doing (-1)^l*Conjugate, frequencies must be negatives.\n");
#endif
for(UINT4 idx=0; idx<(*htildelm)->data->length; idx++){
for(UINT4 idx=offset; idx<(*htildelm)->data->length; idx++){
(*htildelm)->data->data[idx] = minus1l*conj((*htildelm)->data->data[idx]);
}
}
......@@ -1199,6 +1248,11 @@ int XLALSimIMRPhenomXHMMultiBandOneModeMixing(
}
/** @}
* @}
*/
int IMRPhenomXHMMultiBandOneModeMixing(
COMPLEX16FrequencySeries **htildelm, /**< [out] FD waveform */
COMPLEX16FrequencySeries *htilde22, /**< Recycle the 22 mode if previously computed **/
......@@ -1781,7 +1835,8 @@ int IMRPhenomXHMMultiBandOneModeMixing(
}
/* Intialize FrequencySeries for htildelm */