Commit 30493df6 authored by John Douglas Veitch's avatar John Douglas Veitch

Merge branch 'fix_eos' into 'lalinference_o2'

Add option for fixed EOS

See merge request !302
parents cd767089 86f9f592
......@@ -61,6 +61,7 @@
#include <lal/LALString.h>
#include <lal/LALSimInspiral.h>
#include <lal/LALSimInspiralWaveformCache.h>
#include <lal/LALSimNeutronStar.h>
#include <lal/LALHashTbl.h>
#include <lal/SFTutils.h>
......@@ -467,7 +468,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)
......@@ -737,6 +738,7 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
(--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\
......@@ -903,6 +905,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");
......@@ -1350,26 +1353,15 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
}
// For EOS, 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 more than one of --tidalT and --tidal and --4PolyEOS and --4SpectralDecomp.\n");
XLAL_ERROR_NULL(XLAL_EINVAL);
} else if(LALInferenceGetProcParamVal(commandLine,"--tidalT")&&LALInferenceGetProcParamVal(commandLine,"--4PolyEOS")){
XLALPrintError("Error: cannot use more than one of --tidalT and --tidal and --4PolyEOS and --4SpectralDecomp.\n");
XLAL_ERROR_NULL(XLAL_EINVAL);
} else if(LALInferenceGetProcParamVal(commandLine,"--tidalT")&&LALInferenceGetProcParamVal(commandLine,"--4SpectralDecomp")){
XLALPrintError("Error: cannot use more than one of --tidalT and --tidal and --4PolyEOS and --4SpectralDecomp.\n");
XLAL_ERROR_NULL(XLAL_EINVAL);
} else if(LALInferenceGetProcParamVal(commandLine,"--tidal")&&LALInferenceGetProcParamVal(commandLine,"--4PolyEOS")){
XLALPrintError("Error: cannot use more than one of --tidalT and --tidal and --4PolyEOS and --4SpectralDecomp.\n");
XLAL_ERROR_NULL(XLAL_EINVAL);
} else if(LALInferenceGetProcParamVal(commandLine,"--tidal")&&LALInferenceGetProcParamVal(commandLine,"--4SpectralDecomp")){
XLALPrintError("Error: cannot use more than one of --tidalT and --tidal and --4PolyEOS and --4SpectralDecomp.\n");
XLAL_ERROR_NULL(XLAL_EINVAL);
} else if(LALInferenceGetProcParamVal(commandLine,"--4PolyEOS")&&LALInferenceGetProcParamVal(commandLine,"--4SpectralDecomp")){
XLALPrintError("Error: cannot use more than one of --tidalT and --tidal and --4PolyEOS and --4SpectralDecomp.\n");
XLAL_ERROR_NULL(XLAL_EINVAL);
} 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);
......@@ -1389,6 +1381,19 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
LALInferenceRegisterUniformVariableREAL8(state, model->params, "SDgamma2", zero, SDgamma2Min, SDgamma2Max, LALINFERENCE_PARAM_LINEAR);
LALInferenceRegisterUniformVariableREAL8(state, model->params, "SDgamma3", zero, SDgamma3Min, SDgamma3Max, 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>
......@@ -804,6 +805,53 @@ void LALInferenceTemplateXLALSimInspiralChooseWaveform(LALInferenceModel *model)
sym_mass_ratio_eta = m1*m2/((m1+m2)*(m1+m2));
LALInferenceLambdaTsEta2Lambdas(lambdaT,dLambdaT,sym_mass_ratio_eta,&lambda1,&lambda2);
}
if(model->eos_fam)
{
LALSimNeutronStarFamily *eos_fam = model->eos_fam;
REAL8 r1=0, r2=0, k2_1=0, k2_2=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);
}
/* 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.;
......
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