Maintenance will be performed on git.ligo.org, chat.ligo.org, containers.ligo.org, and docs.ligo.org on Tuesday 22nd September 2020 starting at approximately 9am MST.It is expected to take around 15 minutes and there will be a short period of downtime towards the end of the maintenance window. Please address any comments, questions, or concerns to computing-help@igwn.org.

Commit fd55c663 authored by Leslie Wade's avatar Leslie Wade

Changed eos interpolation to steffen (copied from gsl), fixed memory bug

parent 88bd03cf
......@@ -2386,34 +2386,12 @@ else{
double c = mass1 * LAL_MRSUN_SI / r;
*lambda1= (2.0/3.0) * k / pow(c , 5.0);
if(r<0){
printf("Warning: Negative radius, r1 = %e\n",r);
printf("Setting lambda1 = 0.0\n");
*lambda1=0.0;
}
if(k<0){
printf("Warning: Negative love number, k1 = %e\n",k);
printf("Setting lambda1 = 0.0\n");
*lambda1=0.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);
if(r<0){
printf("Warning: Negative radius, r2 = %e\n",r);
printf("Setting lambda2 = 0.0\n");
*lambda2=0.0;
}
if(k<0){
printf("Warning: Negative love number, k2 = %e\n",k);
printf("Setting lambda2 = 0.0\n");
*lambda2=0.0;
}
// Clean up
XLALDestroySimNeutronStarFamily(fam);
XLALDestroySimNeutronStarEOS(eos);
......@@ -2484,9 +2462,9 @@ 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) / 1000.;
double dlogp = (logpmax - logpmin) / 100.;
// Need at least 8 points
for (int i = 0; i < 8; ++i) {
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 */
......@@ -2499,7 +2477,7 @@ for (int i = 0; i < 8; ++i) {
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);
fprintf(stdout,"spectral: %f %f %f %f\n",SDgamma0,SDgamma1,SDgamma2,SDgamma3);
}
// Clean up
XLALDestroySimNeutronStarFamily(fam);
......
......@@ -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>
......@@ -98,7 +99,7 @@ LALSimNeutronStarFamily * XLALCreateSimNeutronStarFamily(
LALSimNeutronStarEOS * eos)
{
LALSimNeutronStarFamily * fam;
const size_t ndatmax = 1000;
const size_t ndatmax = 100;
const double logpmin = 75.5;
double logpmax;
double dlogp;
......@@ -179,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_linear, ndat);
fam->k_of_m_interp = gsl_interp_alloc(gsl_interp_linear, 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);
......
......@@ -320,7 +320,9 @@ liblalsimulation_la_SOURCES = \
LALSimInspiralEOS.c \
LALSimNRHybSurUtilities.c \
LALSimIMRNRHybSur3dq8.c
gsl_interpolation_integ_eval.h \
gsl_interpolation_steffen.c
nodist_liblalsimulation_la_SOURCES = \
LALSimulationBuildInfoHeader.h \
LALSimulationVCSInfo.c \
......
/* 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)
{
const double r1 = a - xi;
const double r2 = b - xi;
const double r12 = r1 + r2;
const double bterm = 0.5 * bi * r12;
const double cterm = (1.0 / 3.0) * ci * (r1 * r1 + r2 * r2 + r1 * r2);
const double dterm = 0.25 * di * r12 * (r1 * r1 + r2 * r2);
return (b - a) * (ai + bterm + cterm + dterm);
}
This diff is collapsed.
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