Commit 367b097f authored by John Douglas Veitch's avatar John Douglas Veitch
Browse files

Merge branch 'fix_eos_master' into 'master'

Fix eos master

See merge request !304
parents 91c2bad0 620ba42d
Pipeline #22931 passed with stages
in 27 minutes and 41 seconds
......@@ -62,6 +62,7 @@
#include <lal/StringInput.h>
#include <lal/LALSimInspiral.h>
#include <lal/LALSimInspiralWaveformCache.h>
#include <lal/LALSimNeutronStar.h>
#include <lal/LALHashTbl.h>
#include <lal/SFTutils.h>
......@@ -466,7 +467,8 @@ typedef struct tagLALInferenceModel
REAL8Window *window; /** A window */
REAL8 padding; /** The padding of the above window */
struct tagLALInferenceROQModel *roq; /** ROQ data */
int roq_flag;
int roq_flag; /** Is ROQ enabled */
LALSimNeutronStarFamily *eos_fam; /** Neutron Star equation of state family */
} LALInferenceModel;
......
......@@ -38,6 +38,7 @@
#include <lal/LALInferenceReadData.h>
#include <lal/LALInferenceInit.h>
#include <lal/LALInferenceCalibrationErrors.h>
#include <lal/LALSimNeutronStar.h>
static int checkParamInList(const char *list, const char *param);
static int checkParamInList(const char *list, const char *param)
......@@ -742,6 +743,8 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
(--tidal) Enables tidal corrections, only with LALSimulation.\n\
(--tidalT) Enables reparmeterized tidal corrections, only with LALSimulation.\n\
(--4PolyEOS) Enables 4-piece polytropic EOS parmeterization, only with LALSimulation.\n\
(--4SpectralDecomp) Enables 4-coeff. spectral decomposition EOS parmeterization, only with LALSimulation.\n\
(--eos EOS) Fix the neutron star EOS. Use \"--eos help\" for allowed names\n\
(--spinOrder PNorder) Specify twice the PN order (e.g. 5 <==> 2.5PN) of spin effects to use, only for LALSimulation (default: -1 <==> Use all spin effects).\n\
(--tidalOrder PNorder) Specify twice the PN order (e.g. 10 <==> 5PN) of tidal effects to use, only for LALSimulation (default: -1 <==> Use all tidal effects).\n\
(--numreldata FileName) Location of NR data file for NR waveforms (with NR_hdf5 approx).\n\
......@@ -896,6 +899,7 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
LALInferenceModel *model = XLALMalloc(sizeof(LALInferenceModel));
model->params = XLALCalloc(1, sizeof(LALInferenceVariables));
memset(model->params, 0, sizeof(LALInferenceVariables));
model->eos_fam = NULL;
UINT4 signal_flag=1;
ppt = LALInferenceGetProcParamVal(commandLine, "--noiseonly");
......@@ -1343,18 +1347,15 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
}
// Must pick to either use tidal, tidalT, or 4-piece polytrope parameters; otherwise throw error message
if(LALInferenceGetProcParamVal(commandLine,"--tidalT")&&LALInferenceGetProcParamVal(commandLine,"--tidal")){
XLALPrintError("Error: cannot use both --tidalT and --tidal.\n");
XLAL_ERROR_NULL(XLAL_EINVAL);
} else if(LALInferenceGetProcParamVal(commandLine,"--tidalT")&&LALInferenceGetProcParamVal(commandLine,"--4PolyEOS")){
XLALPrintError("Error: cannot use both --tidalT and --4PolyEOS.\n");
XLAL_ERROR_NULL(XLAL_EINVAL);
} else if(LALInferenceGetProcParamVal(commandLine,"--tidal")&&LALInferenceGetProcParamVal(commandLine,"--4PolyEOS")){
XLALPrintError("Error: cannot use both --tidal and --4PolyEOS.\n");
XLAL_ERROR_NULL(XLAL_EINVAL);
// Pull in tidalT parameters (lambdaT,dLambdaT)
} else if(LALInferenceGetProcParamVal(commandLine,"--tidalT")){
if((!!LALInferenceGetProcParamVal(commandLine,"--tidalT") + !!LALInferenceGetProcParamVal(commandLine,"--tidal")
+ !!LALInferenceGetProcParamVal(commandLine,"--4PolyEOS") + !!LALInferenceGetProcParamVal(commandLine,"--4SpectralDecomp")
+ !!LALInferenceGetProcParamVal(commandLine,"--eos")) > 1 )
{
XLALPrintError("Error: cannot use more than one of --tidalT, --tidal, --4PolyEOS, --4SpectralDecomp and --eos.\n");
XLAL_ERROR_NULL(XLAL_EINVAL);
}
if(LALInferenceGetProcParamVal(commandLine,"--tidalT")){
LALInferenceRegisterUniformVariableREAL8(state, model->params, "lambdaT", zero, lambdaTMin, lambdaTMax, LALINFERENCE_PARAM_LINEAR);
LALInferenceRegisterUniformVariableREAL8(state, model->params, "dLambdaT", zero, dLambdaTMin, dLambdaTMax, LALINFERENCE_PARAM_LINEAR);
// Pull in tidal parameters (lambda1,lambda2)
......@@ -1368,6 +1369,19 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
LALInferenceRegisterUniformVariableREAL8(state, model->params, "gamma2", zero, gamma2Min, gamma2Max, LALINFERENCE_PARAM_LINEAR);
LALInferenceRegisterUniformVariableREAL8(state, model->params, "gamma3", zero, gamma3Min, gamma3Max, LALINFERENCE_PARAM_LINEAR);
}
else if((ppt=LALInferenceGetProcParamVal(commandLine,"--eos")))
{
LALSimNeutronStarEOS *eos=NULL;
errnum=XLAL_SUCCESS;
XLAL_TRY(eos=XLALSimNeutronStarEOSByName(ppt->value), errnum);
if(errnum!=XLAL_SUCCESS)
XLAL_ERROR_NULL(errnum,"%s: %s",__func__,XLALErrorString(errnum));
XLAL_TRY(model->eos_fam = XLALCreateSimNeutronStarFamily(eos),errnum);
if(errnum!=XLAL_SUCCESS)
XLAL_ERROR_NULL(errnum,"%s: %s",__func__,XLALErrorString(errnum));
if(!model->eos_fam) XLAL_ERROR_NULL(XLAL_EINVAL, "Unable to initialise EOS family");
}
LALSimInspiralSpinOrder spinO = LAL_SIM_INSPIRAL_SPIN_ORDER_ALL;
ppt=LALInferenceGetProcParamVal(commandLine, "--spinOrder");
......
......@@ -40,6 +40,7 @@
#include <lal/LALSimInspiral.h>
#include <lal/LALInferenceTemplate.h>
#include <lal/LALInferenceMultibanding.h>
#include <lal/LALSimNeutronStar.h>
/* LIB imports*/
#include <lal/LALInferenceBurstRoutines.h>
......@@ -74,6 +75,21 @@ const char list_extra_parameters[34][16] = {"dchi0","dchi1","dchi2","dchi3","dch
const UINT4 N_extra_params = 34;
/* Return the quadrupole moment of a neutron star given its lambda
* We use the relations defined here. https://arxiv.org/pdf/1302.4499.pdf.
* Note that the convention we use is that:
* \f$\mathrm{dquadmon} = \bar{Q} - 1.\f$
* Where \f$\bar{Q}\f$ (dimensionless) is the reduced quadrupole moment.
*/
static REAL8 dquadmon_from_lambda(REAL8 lambdav);
static REAL8 dquadmon_from_lambda(REAL8 lambdav)
{
double ll = log(lambdav);
double ai = .194, bi = .0936, ci = 0.0474, di = -4.21e-3, ei = 1.23e-4;
double ln_quad_moment = ai + bi*ll + ci*ll*ll + di*pow(ll,3.0) + ei*pow(ll,4.0);
return(exp(ln_quad_moment) - 1.0);
}
static int InterpolateWaveform(REAL8Vector *freqs, COMPLEX16FrequencySeries *src, COMPLEX16FrequencySeries *dest);
static int InterpolateWaveform(REAL8Vector *freqs, COMPLEX16FrequencySeries *src, COMPLEX16FrequencySeries *dest)
......@@ -793,6 +809,60 @@ void LALInferenceTemplateXLALSimInspiralChooseWaveform(LALInferenceModel *model)
XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
}
if(model->eos_fam)
{
LALSimNeutronStarFamily *eos_fam = model->eos_fam;
REAL8 r1=0, r2=0, k2_1=0, k2_2=0, lambda1=0, lambda2=0;
REAL8 mass_max = XLALSimNeutronStarMaximumMass(eos_fam) / LAL_MSUN_SI;
REAL8 mass_min = XLALSimNeutronStarFamMinimumMass(eos_fam) / LAL_MSUN_SI;
if(m1<mass_max && m1>mass_min)
{
/* Compute l1, l2 from mass and EOS */
r1 = XLALSimNeutronStarRadius(m1*LAL_MSUN_SI, eos_fam);
k2_1 = XLALSimNeutronStarLoveNumberK2(m1*LAL_MSUN_SI, eos_fam);
lambda1 = (2./3.)*k2_1 * pow(r1/(m1*LAL_MRSUN_SI), 5.0);
}
if(m2<mass_max && m2>mass_min)
{
r2 = XLALSimNeutronStarRadius(m2*LAL_MSUN_SI, eos_fam);
k2_2 = XLALSimNeutronStarLoveNumberK2(m2*LAL_MSUN_SI, eos_fam);
lambda2 = (2./3.)*k2_2 * pow(r2/(m2*LAL_MRSUN_SI), 5.0);
}
/* Set waveform params */
XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
REAL8 dQuadMon1 = dquadmon_from_lambda(lambda1);
REAL8 dQuadMon2 = dquadmon_from_lambda(lambda2);
XLALSimInspiralWaveformParamsInsertdQuadMon1(model->LALpars, dQuadMon1);
XLALSimInspiralWaveformParamsInsertdQuadMon2(model->LALpars, dQuadMon2);
/* Calculate maximum frequency */
/* If both lambdas are non-zero compute EOS-dependent f_max */
if((lambda1 > 0) && (lambda2 > 0) && (approximant == TaylorF2))
{
/* Start with ISCO */
f_max = 1. / (pow(6,1.5) * LAL_PI * (m1*LAL_MTSUN_SI + m2*LAL_MTSUN_SI));
REAL8 log_lambda1 = log(lambda1);
REAL8 log_lambda2 = log(lambda2);
REAL8 compactness1 = 0.371 - 0.0391 * log_lambda1 + 0.001056 * log_lambda1 * log_lambda1;
REAL8 compactness2 = 0.371 - 0.0391 * log_lambda2 + 0.001056 * log_lambda2 * log_lambda2;
REAL8 rad1 = m1*LAL_MTSUN_SI / compactness1;
REAL8 rad2 = m2*LAL_MTSUN_SI / compactness2;
/* Use the smaller of ISCO and the EOS-dependent termination */
REAL8 fmax_eos = 1. / LAL_PI * pow((m1*LAL_MTSUN_SI + m2*LAL_MTSUN_SI) / pow((rad1 + rad2),3.0),0.5);
if (fmax_eos < f_max)
{
f_max = fmax_eos;
}
}
/* Add derived quantities for output */
LALInferenceAddVariable(model->params, "radius1", &r1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
LALInferenceAddVariable(model->params, "radius2", &r2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
LALInferenceAddVariable(model->params, "lambda1", &lambda1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
LALInferenceAddVariable(model->params, "lambda2", &lambda2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
}
/* ==== 4-PIECE POLYTROPE EOS PARAMETERS ==== */
REAL8 logp1 = 0.;
......@@ -810,6 +880,7 @@ void LALInferenceTemplateXLALSimInspiralChooseWaveform(LALInferenceModel *model)
LALInferenceLogp1GammasMasses2Lambdas(logp1,gamma1,gamma2,gamma3,m1,m2,&lambda1,&lambda2);
XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
}
/* Only use GR templates */
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment