Commit c45d99d7 authored by Leslie Wade's avatar Leslie Wade

Add gamma check to eos->tidal function

parent 61eb8f43
......@@ -2363,31 +2363,40 @@ XLALDestroySimNeutronStarEOS(eos);
/* Find lambda1,2(m1,2|eos) for spectral EOS model */
void LALInferenceSDGammasMasses2Lambdas(REAL8 gamma[], REAL8 mass1, REAL8 mass2, REAL8 *lambda1, REAL8 *lambda2, int size){
// Convert to SI
double mass1_kg=mass1*LAL_MSUN_SI;
double mass2_kg=mass2*LAL_MSUN_SI;
// Make eos
LALSimNeutronStarEOS *eos = NULL;
LALSimNeutronStarFamily *fam = NULL;
eos = XLALSimNeutronStarEOSSpectralDecomposition(gamma,size);
fam = XLALCreateSimNeutronStarFamily(eos);
// Calculate lambda1(m1|eos)
double r = XLALSimNeutronStarRadius(mass1_kg, fam);
double k = XLALSimNeutronStarLoveNumberK2(mass1_kg, fam);
double c = mass1 * LAL_MRSUN_SI / r;
*lambda1= (2.0/3.0) * k / pow(c , 5.0);
// If unreasonable gammas, do not find lambdas
if(LALInferenceSDGammaCheck(gamma, 4) == XLAL_FAILURE){
*lambda1= 0.;
*lambda2= 0.;
}
// Else calculate lambdas
else{
// Convert to SI
double mass1_kg=mass1*LAL_MSUN_SI;
double mass2_kg=mass2*LAL_MSUN_SI;
// Make eos
LALSimNeutronStarEOS *eos = NULL;
LALSimNeutronStarFamily *fam = NULL;
eos = XLALSimNeutronStarEOSSpectralDecomposition(gamma,size);
fam = XLALCreateSimNeutronStarFamily(eos);
// Calculate lambda1(m1|eos)
double r = XLALSimNeutronStarRadius(mass1_kg, fam);
double k = XLALSimNeutronStarLoveNumberK2(mass1_kg, fam);
double c = mass1 * LAL_MRSUN_SI / r;
*lambda1= (2.0/3.0) * k / pow(c , 5.0);
// Calculate lambda2(m1|eos)
r = XLALSimNeutronStarRadius(mass2_kg, fam);
k = XLALSimNeutronStarLoveNumberK2(mass2_kg, fam);
c = mass2 * LAL_MRSUN_SI / r;
*lambda2= (2.0/3.0) * k / pow(c , 5.0);
// Calculate lambda2(m1|eos)
r = XLALSimNeutronStarRadius(mass2_kg, fam);
k = XLALSimNeutronStarLoveNumberK2(mass2_kg, fam);
c = mass2 * LAL_MRSUN_SI / r;
*lambda2= (2.0/3.0) * k / pow(c , 5.0);
// Clean up
XLALDestroySimNeutronStarFamily(fam);
XLALDestroySimNeutronStarEOS(eos);
}
// Clean up
XLALDestroySimNeutronStarFamily(fam);
XLALDestroySimNeutronStarEOS(eos);
}
......@@ -2461,6 +2470,15 @@ for (int i = 0; i < 4; ++i) {
/* determine if maximum mass has been found */
if (mdat <= mdat_prev){
fprintf(stdout,"EOS has too few points. Sample rejected.\n");
if( LALInferenceCheckVariable(params,"SDgamma0") && LALInferenceCheckVariable(params,"SDgamma1") && LALInferenceCheckVariable(params,"SDgamma2") && LALInferenceCheckVariable(params,"SDgamma3"))
{
// Retrieve EOS params from params linked list
double SDgamma0=*(double *)LALInferenceGetVariable(params,"SDgamma0");
double SDgamma1=*(double *)LALInferenceGetVariable(params,"SDgamma1");
double SDgamma2=*(double *)LALInferenceGetVariable(params,"SDgamma2");
double SDgamma3=*(double *)LALInferenceGetVariable(params,"SDgamma3");
fprintf(stdout,"%f %f %f %f\n",SDgamma0,SDgamma1,SDgamma2,SDgamma3);
}
// Clean up
XLALDestroySimNeutronStarFamily(fam);
XLALDestroySimNeutronStarEOS(eos);
......
......@@ -406,8 +406,7 @@ void LALInferenceROQWrapperForXLALSimInspiralChooseFDWaveformSequence(LALInferen
SDgamma2 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma2");
SDgamma3 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma3");
REAL8 gamma[] = {SDgamma0,SDgamma1,SDgamma2,SDgamma3};
if(LALInferenceSDGammaCheck(gamma, 4) == XLAL_SUCCESS)
LALInferenceSDGammasMasses2Lambdas(gamma,m1,m2,&lambda1,&lambda2,4);
LALInferenceSDGammasMasses2Lambdas(gamma,m1,m2,&lambda1,&lambda2,4);
XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
}
......@@ -938,8 +937,7 @@ void LALInferenceTemplateXLALSimInspiralChooseWaveform(LALInferenceModel *model)
SDgamma2 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma2");
SDgamma3 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma3");
REAL8 gamma[] = {SDgamma0,SDgamma1,SDgamma2,SDgamma3};
if(LALInferenceSDGammaCheck(gamma, 4) == XLAL_SUCCESS)
LALInferenceSDGammasMasses2Lambdas(gamma,m1,m2,&lambda1,&lambda2,4);
LALInferenceSDGammasMasses2Lambdas(gamma,m1,m2,&lambda1,&lambda2,4);
XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
}
......@@ -1476,8 +1474,7 @@ void LALInferenceTemplateXLALSimInspiralChooseWaveformPhaseInterpolated(LALInfer
SDgamma2 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma2");
SDgamma3 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma3");
REAL8 gamma[] = {SDgamma0,SDgamma1,SDgamma2,SDgamma3};
if(LALInferenceSDGammaCheck(gamma, 4) == XLAL_SUCCESS)
LALInferenceSDGammasMasses2Lambdas(gamma,m1,m2,&lambda1,&lambda2,4);
LALInferenceSDGammasMasses2Lambdas(gamma,m1,m2,&lambda1,&lambda2,4);
XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
}
......
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