Maintenance will be performed on git.ligo.org, chat.ligo.org, containers.ligo.org, and docs.ligo.org on the morning of Tuesday 11th August 2020, starting at approximately 9am PDT. It is expected to take around 20 minutes and there will be a short period of downtime (less than five minutes) towards the end of the maintenance window. Please direct any comments, questions, or concerns to computing-help@ligo.org.

Commit f7078061 authored by Jolien Creighton's avatar Jolien Creighton

Merge branch 'spectral_sample' into 'master'

Spectral EOS sampling in LALInferenceMCMC

See merge request !722
parents d2c40344 787929fb
Pipeline #68555 failed with stages
in 83 minutes and 10 seconds
......@@ -1163,7 +1163,7 @@ if __name__=='__main__':
else:
fixedBurnins = None
from lalinference.bayespputils import massParams,spinParams,cosmoParam,strongFieldParams,distParams,incParams,polParams,skyParams,phaseParams,timeParams,endTimeParams,statsParams,calibParams,snrParams,tidalParams,eosParams
from lalinference.bayespputils import massParams,spinParams,cosmoParam,strongFieldParams,distParams,incParams,polParams,skyParams,phaseParams,timeParams,endTimeParams,statsParams,calibParams,snrParams,tidalParams,fourPiecePolyParams,spectralParams
oneDMenus={'Masses':None,'SourceFrame':None,'Timing':None,'Extrinsic':None,'Spins':None,'StrongField':None,'Others':None}
......@@ -1230,12 +1230,19 @@ if __name__=='__main__':
for tp in tidalParams:
if not (mp == tp):
twoDGreedyMenu.append([mp, tp])
for tp in fourPiecePolyParams:
if not (mp == tp):
twoDGreedyMenu.append([mp, tp])
for tp in spectralParams:
if not (mp == tp):
twoDGreedyMenu.append([mp, tp])
for sp1,sp2 in combinations(snrParams,2):
twoDGreedyMenu.append([sp1,sp2])
twoDGreedyMenu.append(['lambda1','lambda2'])
twoDGreedyMenu.append(['lam_tilde','dlam_tilde'])
twoDGreedyMenu.append(['lambdat','dlambdat'])
twoDGreedyMenu.append(['logp1','gamma1','gamma2','gamma3'])
twoDGreedyMenu.append(['SDgamma0','SDgamma1','SDgamma2','SDgamma3'])
for psip in polParams:
for phip in phaseParams:
twoDGreedyMenu.append([psip,phip])
......
......@@ -143,9 +143,10 @@ tigerParams=['dchi%i'%(i) for i in range(8)] + ['dchi%il'%(i) for i in [5,6] ] +
bransDickeParams=['omegaBD','ScalarCharge1','ScalarCharge2']
massiveGravitonParams=['lambdaG']
tidalParams=['lambda1','lambda2','lam_tilde','dlam_tilde','lambdat','dlambdat','lambdas','bluni']
eosParams=['logp1','gamma1','gamma2','gamma3']
fourPiecePolyParams=['logp1','gamma1','gamma2','gamma3']
spectralParams=['sdgamma0','sdgamma1','sdgamma2','sdgamma3']
energyParams=['e_rad', 'l_peak','e_rad_maxldist']
strongFieldParams=ppEParams+tigerParams+bransDickeParams+massiveGravitonParams+tidalParams+eosParams+energyParams
strongFieldParams=ppEParams+tigerParams+bransDickeParams+massiveGravitonParams+tidalParams+fourPiecePolyParams+spectralParams+energyParams
#Extrinsic
distParams=['distance','distMPC','dist','distance_maxl']
......@@ -174,7 +175,9 @@ for param in tigerParams + bransDickeParams + massiveGravitonParams:
greedyBinSizes[param]=0.01
for param in tidalParams:
greedyBinSizes[param]=2.5
for param in eosParams:
for param in fourPiecePolyParams:
greedyBinSizes[param]=2.5
for param in spectralParams:
greedyBinSizes[param]=2.5
#Confidence levels
for loglname in statsParams:
......@@ -357,6 +360,10 @@ def get_prior(name):
'gamma1':None,
'gamma2':None,
'gamma3':None,
'sdgamma0': None,
'sdgamma1': None,
'sdgamma2': None,
'sdgamma3': None,
'calamp_h1' : 'uniform',
'calamp_l1' : 'uniform',
'calpha_h1' : 'uniform',
......@@ -464,6 +471,10 @@ def plot_label(param):
'gamma1':r'$\Gamma_1$',
'gamma2':r'$\Gamma_2$',
'gamma3':r'$\Gamma_3$',
'sdgamma0' : r'$\gamma_{0}$',
'sdgamma1' : r'$\gamma_{1}$',
'sdgamma2' : r'$\gamma_{2}$',
'sdgamma3' : r'$\gamma_{3}$',
'calamp_h1' : r'$\delta A_{H1}$',
'calamp_l1' : r'$\delta A_{L1}$',
'calpha_h1' : r'$\delta \phi_{H1}$',
......
......@@ -2330,7 +2330,7 @@ void LALInferenceLambdaTsEta2Lambdas(REAL8 lambdaT, REAL8 dLambdaT, REAL8 eta, R
return;
}
/* Find lambda1,2(m1,2|eos) */
/* Find lambda1,2(m1,2|eos) for 4-piece polytrope EOS model */
void LALInferenceLogp1GammasMasses2Lambdas(REAL8 logp1,REAL8 gamma1,REAL8 gamma2,REAL8 gamma3, REAL8 mass1, REAL8 mass2, REAL8 *lambda1, REAL8 *lambda2){
// Convert to SI
REAL8 logp1_si=logp1-1.0;
......@@ -2343,13 +2343,13 @@ LALSimNeutronStarFamily *fam = NULL;
eos = XLALSimNeutronStarEOS4ParameterPiecewisePolytrope(logp1_si, gamma1, gamma2, gamma3);
fam = XLALCreateSimNeutronStarFamily(eos);
// Calculating lambda1(m1|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);
// Calculating lambda2(m2|eos)
// Calculate lambda2(m2|eos)
r = XLALSimNeutronStarRadius(mass2_kg, fam);
k = XLALSimNeutronStarLoveNumberK2(mass2_kg, fam);
c = mass2 * LAL_MRSUN_SI / r;
......@@ -2360,6 +2360,46 @@ XLALDestroySimNeutronStarFamily(fam);
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){
// 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);
// Clean up
XLALDestroySimNeutronStarFamily(fam);
XLALDestroySimNeutronStarEOS(eos);
}
}
/* Checks if EOS allows for acausal speed of sound and unphysical maximum masses */
int LALInferenceEOSPhysicalCheck(LALInferenceVariables *params, ProcessParamsTable *commandLine){
int ret;
......@@ -2381,7 +2421,24 @@ if(LALInferenceCheckVariable(params, "logp1") && LALInferenceCheckVariable(param
// Make 4-piece polytrope eos
eos = XLALSimNeutronStarEOS4ParameterPiecewisePolytrope(logp1_si,gamma1,gamma2,gamma3);
fam = XLALCreateSimNeutronStarFamily(eos);
// Else if using 4-coeff spectral eos params...
}
else 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");
double gamma[]={SDgamma0,SDgamma1,SDgamma2,SDgamma3};
if(LALInferenceSDGammaCheck(gamma, 4) == XLAL_FAILURE)
return XLAL_FAILURE;
// Make 4-piece polytrope eos
eos = XLALSimNeutronStarEOSSpectralDecomposition(gamma,4);
}
// Else fail, since you need an eos
else {
......@@ -2389,6 +2446,51 @@ else {
return XLAL_FAILURE;
}
/* FIXME: This is a little clunky,
Check to make sure family will contain
enough pts for interpolation */
double pdat;
double mdat;
double mdat_prev;
double rdat;
double kdat;
/* Initialize previous value for mdat comparison, set to something that will always
make (mdat <= mdat_prev) == true. */
mdat_prev = 0.0;
// Ensure mass turnover does not happen too soon
const double logpmin = 75.5;
double logpmax = log(XLALSimNeutronStarEOSMaxPressure(eos));
double dlogp = (logpmax - logpmin) / 100.;
// Need at least 8 points
for (int i = 0; i < 4; ++i) {
pdat = exp(logpmin + i * dlogp);
XLALSimNeutronStarTOVODEIntegrate(&rdat, &mdat, &kdat, pdat, eos);
/* 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,"spectral: %f %f %f %f\n",SDgamma0,SDgamma1,SDgamma2,SDgamma3);
}
// Clean up
XLALDestroySimNeutronStarFamily(fam);
XLALDestroySimNeutronStarEOS(eos);
return XLAL_FAILURE;
}
mdat_prev = mdat;
}
// Make family
fam = XLALCreateSimNeutronStarFamily(eos);
// Determine which mass parameterization is used
double mass1 = 0.;
double mass2 = 0.;
......@@ -2414,6 +2516,9 @@ else if(LALInferenceCheckVariable(params, "chirpmass") && LALInferenceCheckVaria
else {
// Else fail
fprintf(stdout,"ERROR: NO MASS PARAMETERS FOUND\n");
// Clean up
XLALDestroySimNeutronStarFamily(fam);
XLALDestroySimNeutronStarEOS(eos);
return XLAL_FAILURE;
}
......@@ -2455,6 +2560,65 @@ return ret;
}
// Determines if current spectral parameters are within Adiabatic 'prior' range
int LALInferenceSDGammaCheck(double gamma[], int size) {
int i;
int ret = XLAL_SUCCESS;
double p0 = 4.43784199e-13;
double xmax = 12.3081;
double pmax = p0 * exp(xmax);
size_t ndat = 500;
double *pdat;
double *adat;
double *xdat;
pdat = XLALCalloc(ndat, sizeof(*pdat));
adat = XLALCalloc(ndat, sizeof(*adat));
xdat = XLALCalloc(ndat, sizeof(*xdat));
// Generating higher density portion of EOS with spectral decomposition
double logpmax = log(pmax);
double logp0 = log(p0);
double dlogp = (logpmax-logp0)/ndat;
// Calculating pressure and adiabatic index table
for(i = 0;i < (int) ndat;i++)
{
pdat[i] = exp(logp0 + dlogp*i);
xdat[i] = log(pdat[i]/p0);
adat[i] = AdiabaticIndex(gamma, xdat[i], size);
}
// Make sure all values are within Adiabatic Prior range
for(i = 0;i<(int) ndat;i++) {
// FIXME: Should make range adjustable from the commandline
if(adat[i] < 0.6 || adat[i] > 4.5) {
ret = XLAL_FAILURE;
break;
}
}
LALFree(pdat);
LALFree(xdat);
LALFree(adat);
LALCheckMemoryLeaks();
return ret;
}
/* Specral decomposition of eos's adiabatic index */
double AdiabaticIndex(double gamma[],double x, int size)
{
double Gamma, logGamma = 0;
int i;
for(i=0;i<size;i++)
{
logGamma += gamma[i]*pow(x,(double) i);
}
Gamma = exp(logGamma);
return Gamma;
}
static void deleteCell(LALInferenceKDTree *cell) {
if (cell == NULL) {
return; /* Our work here is done. */
......
......@@ -888,9 +888,18 @@ void LALInferenceLambdaTsEta2Lambdas(REAL8 lambdaT, REAL8 dLambdaT, REAL8 eta, R
/** Calculate lambda1,2(m1,2|eos(logp1,gamma1,gamma2,gamma3)) */
void LALInferenceLogp1GammasMasses2Lambdas(REAL8 logp1, REAL8 gamma1, REAL8 gamma2, REAL8 gamma3, REAL8 mass1, REAL8 mass2, REAL8 *lambda1, REAL8 *lambda2);
/** Check for causality violation and mass conflict given masses and eos **/
/** Convert from spectral parameters to lambda1, lambda2 */
void LALInferenceSDGammasMasses2Lambdas(REAL8 gamma[], REAL8 mass1, REAL8 mass2, REAL8 *lambda1, REAL8 *lambda2, int size);
/** Check for causality violation and mass conflict given masses and eos */
int LALInferenceEOSPhysicalCheck(LALInferenceVariables *params, ProcessParamsTable *commandLine);
/** Specral decomposition of eos's adiabatic index */
double AdiabaticIndex(double gamma[],double x, int size);
/** Determine if the Adiabatic index is within 'prior' range */
int LALInferenceSDGammaCheck(double gamma[], int size);
/**
* The kD trees in LALInference are composed of cells. Each cell
* represents a rectangular region in parameter space, defined by
......
......@@ -792,6 +792,11 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
gamma3 gamma3.\n\
(requires --BinaryLove)\n\
lambdaS Symmetric tidal deformability.\n\
(requires --4SpectralDecomp):\n\
SDgamma0 SDgamma0.\n\
SDgamma1 SDgamma1.\n\
SDgamma2 SDgamma2.\n\
SDgamma3 SDgamma3.\n\
* \n\
----------------------------------------------\n\
--- Prior Ranges -----------------------------\n\
......@@ -882,6 +887,14 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
REAL8 gamma3Max=4.5;
REAL8 lambdaSMin=0.0;
REAL8 lambdaSMax=5000.0;
REAL8 SDgamma0Min=0.2;
REAL8 SDgamma0Max=2.0;
REAL8 SDgamma1Min=-1.6;
REAL8 SDgamma1Max=1.7;
REAL8 SDgamma2Min=-0.6;
REAL8 SDgamma2Max=0.6;
REAL8 SDgamma3Min=-0.02;
REAL8 SDgamma3Max=0.02;
gsl_rng *GSLrandom=state->GSLrandom;
REAL8 endtime=0.0, timeParam=0.0;
REAL8 timeMin=endtime-dt,timeMax=endtime+dt;
......@@ -1387,6 +1400,12 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
LALInferenceRegisterUniformVariableREAL8(state, model->params, "gamma1", zero, gamma1Min, gamma1Max, LALINFERENCE_PARAM_LINEAR);
LALInferenceRegisterUniformVariableREAL8(state, model->params, "gamma2", zero, gamma2Min, gamma2Max, LALINFERENCE_PARAM_LINEAR);
LALInferenceRegisterUniformVariableREAL8(state, model->params, "gamma3", zero, gamma3Min, gamma3Max, LALINFERENCE_PARAM_LINEAR);
// Pull in spectral decomposition parameters (SDgamma0,SDgamma1,SDgamma2,SDgamma3)
} else if(LALInferenceGetProcParamVal(commandLine,"--4SpectralDecomp")){
LALInferenceRegisterUniformVariableREAL8(state, model->params, "SDgamma0", zero, SDgamma0Min, SDgamma0Max, LALINFERENCE_PARAM_LINEAR);
LALInferenceRegisterUniformVariableREAL8(state, model->params, "SDgamma1", zero, SDgamma1Min, SDgamma1Max, LALINFERENCE_PARAM_LINEAR);
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;
......
......@@ -634,6 +634,13 @@ REAL8 LALInferenceInspiralPrior(LALInferenceRunState *runState, LALInferenceVari
return -INFINITY;
}
}
else if((LALInferenceCheckVariable(params,"SDgamma0")&&LALInferenceCheckVariable(params,"SDgamma1")&&LALInferenceCheckVariable(params,"SDgamma2")&&LALInferenceCheckVariable(params,"SDgamma3")))
{
/*If EOS params and masses are aphysical, return -INFINITY to ensure point is rejected*/
if(LALInferenceEOSPhysicalCheck(params,runState->commandLine)==XLAL_FAILURE){
return -INFINITY;
}
}
}/* end prior for signal model parameters */
......
......@@ -96,7 +96,8 @@ const char *const splineCalibrationProposalName = "SplineCalibration";
const char *const distanceLikelihoodProposalName = "DistanceLikelihood";
static const char *intrinsicNames[] = {"chirpmass", "q", "eta", "mass1", "mass2", "a_spin1", "a_spin2",
"tilt_spin1", "tilt_spin2", "phi12", "phi_jl", "frequency", "quality", "duration","polar_angle", "phase", "polar_eccentricity","dchi0","dchi1","dchi2","dchi3","dchi4","dchi5","dchi5l","dchi6","dchi6l","dchi7","aPPE","alphaPPE","bPPE","betaPPE","betaStep","fStep","dxi1","dxi2","dxi3","dxi4","dxi5","dxi6","dalpha1","dalpha2","dalpha3","dalpha4","dalpha5","dbeta1","dbeta2","dbeta3","dsigma1","dsigma2","dsigma3","dsigma4",NULL};
"tilt_spin1", "tilt_spin2", "phi12", "phi_jl", "frequency", "quality", "duration","polar_angle", "phase", "polar_eccentricity","dchi0","dchi1","dchi2","dchi3","dchi4","dchi5","dchi5l","dchi6","dchi6l","dchi7","aPPE","alphaPPE","bPPE","betaPPE","betaStep","fStep","dxi1","dxi2","dxi3","dxi4","dxi5","dxi6","dalpha1","dalpha2","dalpha3","dalpha4","dalpha5","dbeta1","dbeta2","dbeta3","dsigma1","dsigma2","dsigma3","dsigma4","lambda1","lambda2","lambdaT","dlambdaT","logp1", "gamma1", "gamma2", "gamma3", "SDgamma0","SDgamma1","SDgamma2","SDgamma3",NULL};
static const char *extrinsicNames[] = {"rightascension", "declination", "cosalpha", "azimuth", "polarisation", "distance",
"logdistance", "time", "costheta_jn", "t0", "theta","hrss", "loghrss", NULL};
......@@ -1274,7 +1275,8 @@ REAL8 LALInferenceDrawApproxPrior(LALInferenceThreadState *thread,
const char *flat_params[] = {"q", "eta", "t0", "azimuth", "cosalpha", "time", "phase", "polarisation",
"rightascension", "costheta_jn", "phi_jl",
"phi12", "a_spin1", "a_spin2", "logp1", "gamma1", "gamma2", "gamma3", NULL};
"phi12", "a_spin1", "a_spin2", "logp1", "gamma1", "gamma2", "gamma3",
"SDgamma0","SDgamma1","SDgamma2","SDgamma3", NULL};
LALInferenceVariables *args = thread->proposalArgs;
......
......@@ -536,6 +536,10 @@ static void LALInferencePrintDataWithInjection(LALInferenceIFOData *IFOdata, Pro
(--inj-gamma1) value of gamma1 to be injected\n\
(--inj-gamma2) value of gamma2 to be injected\n\
(--inj-gamma3) value of gamma3 to be injected\n\
(--inj-SDgamma0) value of SDgamma0 to be injected (0)\n\
(--inj-SDgamma1) value of SDgamma1 to be injected (0)\n\
(--inj-SDgamma2) value of SDgamma2 to be injected (0)\n\
(--inj-SDgamma3) value of SDgamma3 to be injected (0)\n\
(--inj-spinOrder PNorder) Specify twice the injection PN order (e.g. 5 <==> 2.5PN)\n\
of spin effects effects to use, only for LALSimulation\n\
(default: -1 <==> Use all spin effects).\n\
......@@ -1625,6 +1629,27 @@ void LALInferenceInjectInspiralSignal(LALInferenceIFOData *IFOdata, ProcessParam
*/
}
// Inject 4-coef. spectral eos
REAL8 SDgamma0=0.0;
REAL8 SDgamma1=0.0;
REAL8 SDgamma2=0.0;
REAL8 SDgamma3=0.0;
if(LALInferenceGetProcParamVal(commandLine,"--inj-SDgamma0") && LALInferenceGetProcParamVal(commandLine,"--inj-SDgamma1") && LALInferenceGetProcParamVal(commandLine,"--inj-SDgamma2") && LALInferenceGetProcParamVal(commandLine,"--inj-SDgamma3")){
SDgamma0= atof(LALInferenceGetProcParamVal(commandLine,"--inj-SDgamma0")->value);
SDgamma1= atof(LALInferenceGetProcParamVal(commandLine,"--inj-SDgamma1")->value);
SDgamma2= atof(LALInferenceGetProcParamVal(commandLine,"--inj-SDgamma2")->value);
SDgamma3= atof(LALInferenceGetProcParamVal(commandLine,"--inj-SDgamma3")->value);
REAL8 gamma[]={SDgamma0,SDgamma1,SDgamma2,SDgamma3};
LALInferenceSDGammasMasses2Lambdas(gamma,m1,m2,&lambda1,&lambda2,4);
fprintf(stdout,"Injection SDgamma0 set to %lf\n",SDgamma0);
fprintf(stdout,"Injection SDgamma1 set to %lf\n",SDgamma1);
fprintf(stdout,"Injection SDgamma2 set to %lf\n",SDgamma2);
fprintf(stdout,"Injection SDgamma3 set to %lf\n",SDgamma3);
fprintf(stdout,"lambda1 set to %f\n",lambda1);
fprintf(stdout,"lambda2 set to %f\n",lambda2);
}
REAL8 fref = 100.;
if(LALInferenceGetProcParamVal(commandLine,"--inj-fref")) {
fref = atoi(LALInferenceGetProcParamVal(commandLine,"--inj-fref")->value);
......@@ -1894,6 +1919,27 @@ void InjectFD(LALInferenceIFOData *IFOdata, SimInspiralTable *inj_table, Process
*/
}
// Inject 4-coef. spectral eos
REAL8 SDgamma0=0.0;
REAL8 SDgamma1=0.0;
REAL8 SDgamma2=0.0;
REAL8 SDgamma3=0.0;
if(LALInferenceGetProcParamVal(commandLine,"--inj-SDgamma0") && LALInferenceGetProcParamVal(commandLine,"--inj-SDgamma1") && LALInferenceGetProcParamVal(commandLine,"--inj-SDgamma2") && LALInferenceGetProcParamVal(commandLine,"--inj-SDgamma3")){
SDgamma0= atof(LALInferenceGetProcParamVal(commandLine,"--inj-SDgamma0")->value);
SDgamma1= atof(LALInferenceGetProcParamVal(commandLine,"--inj-SDgamma1")->value);
SDgamma2= atof(LALInferenceGetProcParamVal(commandLine,"--inj-SDgamma2")->value);
SDgamma3= atof(LALInferenceGetProcParamVal(commandLine,"--inj-SDgamma3")->value);
REAL8 gamma[]={SDgamma0,SDgamma1,SDgamma2,SDgamma3};
LALInferenceSDGammasMasses2Lambdas(gamma,inj_table->mass1,inj_table->mass2,&lambda1,&lambda2,4);
fprintf(stdout,"Injection SDgamma0 set to %lf\n",SDgamma0);
fprintf(stdout,"Injection SDgamma1 set to %lf\n",SDgamma1);
fprintf(stdout,"Injection SDgamma2 set to %lf\n",SDgamma2);
fprintf(stdout,"Injection SDgamma3 set to %lf\n",SDgamma3);
fprintf(stdout,"lambda1 set to %f\n",lambda1);
fprintf(stdout,"lambda2 set to %f\n",lambda2);
}
/* FIXME: One also need to code the same manipulations done to f_max in LALInferenceTemplateXLALSimInspiralChooseWaveform line 856 circa*/
/* Set up LAL dictionary */
......
......@@ -392,6 +392,25 @@ void LALInferenceROQWrapperForXLALSimInspiralChooseFDWaveformSequence(LALInferen
XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
}
/* ==== SPECTRAL DECOMPOSITION PARAMETERS ==== */
REAL8 SDgamma0 = 0.;
REAL8 SDgamma1 = 0.;
REAL8 SDgamma2 = 0.;
REAL8 SDgamma3 = 0.;
/* Checks for 4 spectral parameters */
if(!LALInferenceCheckVariable(model->params, "logp1")&&LALInferenceCheckVariable(model->params, "SDgamma0")&&LALInferenceCheckVariable(model->params, "SDgamma1")&&LALInferenceCheckVariable(model->params, "SDgamma2")&&LALInferenceCheckVariable(model->params,"SDgamma3")){
REAL8 lambda1 = 0.;
REAL8 lambda2 = 0.;
SDgamma0 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma0");
SDgamma1 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma1");
SDgamma2 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma2");
SDgamma3 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma3");
REAL8 gamma[] = {SDgamma0,SDgamma1,SDgamma2,SDgamma3};
LALInferenceSDGammasMasses2Lambdas(gamma,m1,m2,&lambda1,&lambda2,4);
XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
}
/* ==== BINARY_LOVE PARAMETERS ==== */
if(LALInferenceCheckVariable(model->params, "lambdaS")){
REAL8 lambda1=0.;
......@@ -902,7 +921,25 @@ void LALInferenceTemplateXLALSimInspiralChooseWaveform(LALInferenceModel *model)
LALInferenceLogp1GammasMasses2Lambdas(logp1,gamma1,gamma2,gamma3,m1,m2,&lambda1,&lambda2);
XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
}
/* ==== SPECTRAL DECOMPOSITION PARAMETERS ==== */
REAL8 SDgamma0 = 0.;
REAL8 SDgamma1 = 0.;
REAL8 SDgamma2 = 0.;
REAL8 SDgamma3 = 0.;
/* Checks for 4 spectral parameters */
if(!LALInferenceCheckVariable(model->params, "logp1")&&LALInferenceCheckVariable(model->params, "SDgamma0")&&LALInferenceCheckVariable(model->params, "SDgamma1")&&LALInferenceCheckVariable(model->params, "SDgamma2")&&LALInferenceCheckVariable(model->params,"SDgamma3")){
REAL8 lambda1 = 0.;
REAL8 lambda2 = 0.;
SDgamma0 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma0");
SDgamma1 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma1");
SDgamma2 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma2");
SDgamma3 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma3");
REAL8 gamma[] = {SDgamma0,SDgamma1,SDgamma2,SDgamma3};
LALInferenceSDGammasMasses2Lambdas(gamma,m1,m2,&lambda1,&lambda2,4);
XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
}
/* ==== BINARY_LOVE PARAMETERS ==== */
......@@ -1424,6 +1461,25 @@ void LALInferenceTemplateXLALSimInspiralChooseWaveformPhaseInterpolated(LALInfer
XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
}
/* ==== SPECTRAL DECOMPOSITION PARAMETERS ==== */
REAL8 SDgamma0 = 0.;
REAL8 SDgamma1 = 0.;
REAL8 SDgamma2 = 0.;
REAL8 SDgamma3 = 0.;
if(!LALInferenceCheckVariable(model->params, "logp1")&&LALInferenceCheckVariable(model->params, "SDgamma0")&&LALInferenceCheckVariable(model->params, "SDgamma1")&&LALInferenceCheckVariable(model->params, "SDgamma2")&&LALInferenceCheckVariable(model->params,"SDgamma3")){
REAL8 lambda1 = 0.;
REAL8 lambda2 = 0.;
SDgamma0 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma0");
SDgamma1 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma1");
SDgamma2 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma2");
SDgamma3 = *(REAL8*) LALInferenceGetVariable(model->params,"SDgamma3");
REAL8 gamma[] = {SDgamma0,SDgamma1,SDgamma2,SDgamma3};
LALInferenceSDGammasMasses2Lambdas(gamma,m1,m2,&lambda1,&lambda2,4);
XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
}
/* Only use GR templates */
/* Fill in the extra parameters for testing GR, if necessary */
for (UINT4 k=0; k<N_extra_params; k++)
......
......@@ -352,6 +352,12 @@ adapt-temps=
#Common EoS tidal model for BNS systems (requires BinaryLove=)
#lambdaS Symmetric tidal deformability, (lambda1+lambda2)/2 [0,5000]
#Equation of State parameters (requires 4SpectralDecomp=):
#SDgamma0 SDgamma0 [0.2,2.0]
#SDgamma1 SDgamma1 [-1.6,1.7]
#SDgamma2 SDgamma2 [-0.6,0.6]
#SDgamma3 SDgamma3 [-0.02,0.02]
# Settings for allowed component masses in Solar Masses, with default values
#comp-max=30.0
#comp-min=1.0
......
......@@ -285,6 +285,8 @@ LALSimNeutronStarEOS *XLALSimNeutronStarEOSSpectralDecomposition(double gamma[],
gamma[0], gamma[1], gamma[2], gamma[3]) >= (int) sizeof(eos->name))
XLAL_PRINT_WARNING("EOS name too long");
LALFree(edat);
LALFree(pdat);
return eos;
}
......
......@@ -28,6 +28,7 @@
#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_min.h>
GSL_VAR const gsl_interp_type * lal_gsl_interp_steffen;
#include <lal/LALStdlib.h>
#include <lal/LALSimNeutronStar.h>
......@@ -157,10 +158,14 @@ LALSimNeutronStarFamily * XLALCreateSimNeutronStarFamily(
&fam->kdat[i], fam->pdat[i], eos);
/* resize arrays */
if(fam->pdat[i] == fam->pdat[i-1])
if(fam->pdat[i] <= fam->pdat[i-1]){
fam->pdat[i-1] = fam->pdat[i];
ndat = i;
else
}
else{
ndat = i + 1;
}
fam->pdat = LALRealloc(fam->pdat, ndat * sizeof(*fam->pdat));
fam->mdat = LALRealloc(fam->mdat, ndat * sizeof(*fam->mdat));
fam->rdat = LALRealloc(fam->rdat, ndat * sizeof(*fam->rdat));
......@@ -175,8 +180,8 @@ LALSimNeutronStarFamily * XLALCreateSimNeutronStarFamily(
fam->k_of_m_acc = gsl_interp_accel_alloc();
fam->p_of_m_interp = gsl_interp_alloc(gsl_interp_cspline, ndat);
fam->r_of_m_interp = gsl_interp_alloc(gsl_interp_cspline, ndat);
fam->k_of_m_interp = gsl_interp_alloc(gsl_interp_cspline, ndat);
fam->r_of_m_interp = gsl_interp_alloc(lal_gsl_interp_steffen, ndat);
fam->k_of_m_interp = gsl_interp_alloc(lal_gsl_interp_steffen, ndat);
gsl_interp_init(fam->p_of_m_interp, fam->mdat, fam->pdat, ndat);
gsl_interp_init(fam->r_of_m_interp, fam->mdat, fam->rdat, ndat);
......
......@@ -319,7 +319,9 @@ liblalsimulation_la_SOURCES = \
LALSimUniversalRelations.c \
LALSimInspiralEOS.c \
LALSimNRHybSurUtilities.c \
LALSimIMRNRHybSur3dq8.c
LALSimIMRNRHybSur3dq8.c \
gsl_interpolation_integ_eval.h \
gsl_interpolation_steffen.c
nodist_liblalsimulation_la_SOURCES = \
LALSimulationBuildInfoHeader.h \
......
/* This file is a copy of GSL's
* interpolation/integ_eval_macro.h
*
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
/* function for doing the spline integral evaluation
which is common to both the cspline and akima methods
*/
static inline double
integ_eval (double ai, double bi, double ci, double di, double xi, double a,
double b)
{