migrate SEOBNRv4p branch

Dockerfile

-FROM ligo/lalsuite-runtime:jessie
-COPY /opt/lalsuite /opt/lalsuite
-ENV LD_LIBRARY_PATH="/opt/lalsuite/lib" \
-    OCTAVE_PATH="/opt/lalsuite/lib/x86_64-linux-gnu/octave/3.8.2/site/oct/x86_64-pc-linux-gnu" \
-    PATH="/opt/lalsuite/bin:${PATH}" \
-    PKG_CONFIG_PATH="/opt/lalsuite/lib/pkgconfig" \
-    PYTHONPATH="/opt/lalsuite/lib/python2.7/site-packages" \

README.md

-# LALSuite
-
-This is the main LALSuite development repository, past development is
-captured in the repository below:
-
-https://git.ligo.org/lscsoft/lalsuite-archive
-
-This new repository utilizes [git-lfs](https://wiki.ligo.org/DASWG/GitLFS#Install_the_git_LFS_client) for the managament of large files and as such `git-lfs` needs to be installed and configured to correctly clone this repository. After installing `git-lfs` it can be configured using:
-
-```
-$ git lfs install
-```
-
-This only needs to be done once for each machine you access the repository. It can then be cloned using:
-
-```
-$ git clone git@git.ligo.org:lscsoft/lalsuite.git
-```

lalframe/src/fmt.c

@@ -113,7 +113,7 @@ int main(int argc, char *argv[])
             fp = stdin;
         else
             fp = fopen(fname, "r");
-        if (!fp)
+        if (!output)
             FAILURE("could not open file %s\n", fname);
 
         readdata(&ndim, &dimlen, &data, &channames, &chanunits, fp);

lalinference/python/bayestar_samples_ppplot.py

@@ -22,6 +22,7 @@ Create a P-P plot to compare a posterior sample chain with a sky map.
 # Command line interface.
 from argparse import FileType
 from lalinference.bayestar import command
+parser = command.ArgumentParser()
 parser = command.ArgumentParser(parents=[command.figure_parser])
 parser.add_argument(
     'skymap', metavar='SKYMAP.fits[.gz]', type=FileType('rb'),

lalinference/python/lalinference/bayestar/command.py

@@ -54,11 +54,6 @@ import distutils.version
 mpl_version = distutils.version.LooseVersion(matplotlib.__version__)
 
 
-def get_version():
-    from .. import InferenceVCSInfo as vcs_info
-    return vcs_info.name + ' ' + vcs_info.version
-
-
 @contextlib.contextmanager
 def TemporaryDirectory(suffix='', prefix='tmp', dir=None, delete=True):
     try:
@@ -174,16 +169,7 @@ class MatplotlibFigureType(argparse.FileType):
 
     def __save(self):
         from matplotlib import pyplot as plt
-        _, ext = os.path.splitext(self.string)
-        ext = ext.lower()
-        program, _ = os.path.splitext(os.path.basename(sys.argv[0]))
-        cmdline = ' '.join([program] + sys.argv[1:])
-        metadata = {'Title': cmdline}
-        if ext == '.png':
-            metadata['Software'] = get_version()
-        elif ext in {'.pdf', '.ps', '.eps'}:
-            metadata['Creator'] = get_version()
-        return plt.savefig(self.string, metadata=metadata)
+        return plt.savefig(self.string)
 
     def __call__(self, string):
         from matplotlib import pyplot as plt
@@ -298,7 +284,8 @@ del group
 # Defer loading SWIG bindings until version string is needed.
 class VersionAction(argparse._VersionAction):
     def __call__(self, parser, namespace, values, option_string=None):
-        self.version = get_version()
+        from .. import InferenceVCSInfo
+        self.version = 'LALInference ' + InferenceVCSInfo.version
         super(VersionAction, self).__call__(
             parser, namespace, values, option_string)
 

lalinference/python/lalinference/rapid_pe/lalsimutils.py

@@ -109,9 +109,6 @@ class ChooseWaveformParams:
         for k, p in ChooseWaveformParams._LAL_DICT_PARAMS.iteritems():
             typfunc = ChooseWaveformParams._LAL_DICT_PTYPE[k]
             typfunc(extra_params, k, getattr(self, p))
-        # Properly add tidal parammeters
-        lalsim.SimInspiralWaveformParamsInsertTidalLambda1(extra_params, self.lambda1)
-        lalsim.SimInspiralWaveformParamsInsertTidalLambda2(extra_params, self.lambda2)
         return extra_params
 
     def copy(self):

lalinference/python/rapidpe_compute_intrinsic_grid.py

@@ -42,9 +42,9 @@ from glue.ligolw.utils import process
 import lalsimulation
 from lalinference.rapid_pe import amrlib, lalsimutils, common_cl
 
-def get_cr_from_grid(cells, weight, cr_thr=0.9, min_n=None, max_n=None):
+def get_cr_from_grid(cells, weight, cr_thr=0.9):
     """
-    Given a set of cells and the weight of that cell, calculate a N% CR including cells which contribute to that probability mass. If n is set, cr_thr is ignored and instead this many points are taken.
+    Given a set of cells and the weight of that cell, calculate a N% CR including cells which contribute to that probability mass.
     """
     if cr_thr == 0.0:
         return numpy.empty((0,))
@@ -57,13 +57,8 @@ def get_cr_from_grid(cells, weight, cr_thr=0.9, min_n=None, max_n=None):
     cell_sort[:,0] = cell_sort[:,0].cumsum()
     cell_sort[:,0] /= cell_sort[-1,0]
 
+    # find the CR probability
     idx = cell_sort[:,0].searchsorted(1-cr_thr)
-    n_select = cell_sort.shape[0] - idx
-    if min_n is not None:
-        n_select = max(n_select, min_n)
-    if max_n is not None:
-        n_select = min(n_select, max_n)
-    idx = cell_sort.shape[0] - n_select
 
     return cell_sort[idx:,1:]
 
@@ -90,7 +85,7 @@ def find_olap_index(tree, intr_prms, exact=True, **kwargs):
     pt = numpy.array([kwargs[k] for k in intr_prms])
 
     # FIXME: Replace with standard function
-    dist, m_idx = tree.query(numpy.atleast_2d(pt), k=1)
+    dist, m_idx = tree.query(pt, k=1)
     dist, m_idx = dist[0][0], int(m_idx[0][0])
 
     # FIXME: There's still some tolerance from floating point conversions
@@ -109,11 +104,6 @@ def write_to_xml(cells, intr_prms, pin_prms={}, fvals=None, fname=None, verbose=
     process.append_process_params(xmldoc, procrow, process.process_params_from_dict(opts.__dict__))
 
     rows = ["simulation_id", "process_id", "numrel_data"]
-    # Override eff_lambda to with psi0, its shoehorn column
-    if "eff_lambda" in intr_prms:
-        intr_prms[intr_prms.index("eff_lambda")] = "psi0"
-    if "deff_lambda" in intr_prms:
-        intr_prms[intr_prms.index("deff_lambda")] = "psi3"
     rows += list(intr_prms)
     rows += list(pin_prms)
     if fvals is not None:
@@ -152,7 +142,7 @@ def get_evidence_grid(points, res_pts, intr_prms, exact=False):
     grid_idx = []
     # Reorder the grid points to match their weight indices
     for res in res_pts:
-        dist, idx = grid_tree.query(numpy.atleast_2d(res), k=1)
+        dist, idx = grid_tree.query(res, k=1)
         # Stupid floating point inexactitude...
         #print res, selected[idx[0][0]]
         #assert numpy.allclose(res, selected[idx[0][0]])
@@ -192,16 +182,12 @@ grid_section.add_argument("--setup", help="Set up the initial grid based on temp
 grid_section.add_argument("-t", "--tmplt-bank", help="XML file with template bank.")
 grid_section.add_argument("-O", "--use-overlap", help="Use overlap information to define 'closeness'.")
 grid_section.add_argument("-T", "--overlap-threshold", type=float, help="Threshold on overlap value.")
-grid_section.add_argument("-s", "--points-per-side", type=int, default=10, help="Number of points per side, default is 10.")
-grid_section.add_argument("-I", "--initial-region", action="append", help="Override the initial region with a custom specification. Specify multiple times like, -I mass1=1.0,2.0 -I mass2=1.0,1.5")
 grid_section.add_argument("-D", "--deactivate", action="store_true", help="Deactivate cells initially which have no template within them.")
 grid_section.add_argument("-P", "--prerefine", help="Refine this initial grid based on overlap values.")
 
 refine_section = argp.add_argument_group("refine options", "Options for refining a pre-existing grid.")
 refine_section.add_argument("--refine", help="Refine a prexisting grid. Pass this option the grid points from previous levels (or the --setup) option.")
 refine_section.add_argument("-r", "--result-file", help="XML file containing newest result to refine.")
-refine_section.add_argument("-M", "--max-n-points", help="Refine *at most* this many points, can override confidence region thresholds.")
-refine_section.add_argument("-m", "--min-n-points", help="Refine *at least* this many points, can override confidence region thresholds.")
 
 opts = argp.parse_args()
 
@@ -340,21 +326,15 @@ if opts.refine or opts.prerefine:
     init_region, region_labels = amrlib.load_init_region(opts.refine or opts.prerefine, get_labels=True)
 else:
     ####### BEGIN INITIAL GRID CODE #########
-    if opts.initial_region is None:
     init_region, idx = determine_region(pt, pts, ovrlp, opts.overlap_threshold, expand_prms)
     region_labels = intr_prms
     # FIXME: To be reimplemented in a different way
     #if opts.expand_param is not None:
         #expand_param(init_region, opts.expand_param)
-    else:
-        # Override initial region -- use with care
-        _, init_region = common_cl.parse_param(opts.initial_region)
-        region_labels = init_region.keys()
-        init_region = amrlib.Cell(numpy.vstack(init_region[k] for k in region_labels))
 
     # TODO: Alternatively, check density of points in the region to determine
     # the points to a side
-    grid, spacing = amrlib.create_regular_grid_from_cell(init_region, side_pts=opts.points_per_side / 2, return_cells=True)
+    grid, spacing = amrlib.create_regular_grid_from_cell(init_region, side_pts=5, return_cells=True)
 
     # "Deactivate" cells not close to template points
     # FIXME: This gets more and more dangerous in higher dimensions
@@ -363,7 +343,7 @@ else:
     if opts.deactivate:
         get_idx = set()
         for pt in pts[idx:]:
-            get_idx.add(tree.query(numpy.atleast_2d(pt), k=1, return_distance=False)[0][0])
+            get_idx.add(tree.query(pt, k=1, return_distance=False)[0][0])
         selected = grid[numpy.array(list(get_idx))]
     else:
         selected = grid
@@ -395,7 +375,7 @@ if opts.result_file is not None:
 
     if opts.refine:
         # FIXME: We use overlap threshold as a proxy for confidence level
-        selected = get_cr_from_grid(selected, results, cr_thr=opts.overlap_threshold, min_n=opts.min_n_points, max_n=opts.max_n_points)
+        selected = get_cr_from_grid(selected, results, cr_thr=opts.overlap_threshold)
         print "Selected %d cells from %3.2f%% confidence region" % (len(selected), opts.overlap_threshold*100)
 
 if opts.prerefine:

lalsimulation/src/LALSimBlackHoleRingdownPrec.c

@@ -9,6 +9,7 @@
 //#include <lal/TimeSeries.h>
 //#include <lal/Units.h>
 
+#include <LALSimBlackHoleRingdown.h>
 #include <LALSimBlackHoleRingdownPrec.h>
 
 /* note: use double-precision variables, but demand single-precision accuracy */
@@ -42,6 +43,8 @@ INT4 XLALSimIMREOBFinalMassSpinPrec(
   Approximant     approximant /**<< The waveform approximant being used */
 )
 {
+  REAL8  LISCO,csi,aeff, a_tot_prec,fitpart,alpha,beta,gamma,epsilon_beta,epsilon_gamma;
+  UNUSED REAL8 epsilon_alpha;
   static const REAL8 root9ovr8minus1 = -0.057190958417936644;
   static const REAL8 root12          = 3.4641016151377544;
 
@@ -130,6 +133,7 @@ INT4 XLALSimIMREOBFinalMassSpinPrec(
       //  Precessing spins (Barausse & Rezzolla 2009, Eqs (6-10))
       chi1   = sqrt( spin1[0]*spin1[0] + spin1[1]*spin1[1] + spin1[2]*spin1[2] );
       chi2   = sqrt( spin2[0]*spin2[0] + spin2[1]*spin2[1] + spin2[2]*spin2[2] );
+      
       if ( chi1 < 1.0e-15 )
       { theta1 = 0.; }
       else
@@ -140,7 +144,40 @@ INT4 XLALSimIMREOBFinalMassSpinPrec(
       { theta2 = acos( spin2[2] / chi2 ); }
       phi1   = atan2( spin1[1], spin1[0] );
       phi2   = atan2( spin2[1], spin2[0] );
-      if ( mass1 > mass2 )
+
+
+      alpha = 0.0;
+      beta = 0.0;
+      gamma = 0.0;
+
+      if (chi1>1e-6){
+	beta=acos(spin1[2]/chi1);
+      }
+
+      if (chi2>1e-6){
+        gamma=acos(spin2[2]/chi2);
+      }
+
+     
+          REAL8 tmp = spin1[0]*spin2[0]+spin1[1]*spin2[1]+spin1[2]*spin2[2];
+          tmp /= chi1;
+          tmp /= chi2;
+
+      if (chi1>1e-4 && chi2>1e-4){
+          if ( tmp > 1. && tmp <= 1. + 1e-4) tmp = 1.;
+          if ( tmp < -1. && tmp >= -1. - 1e-4) tmp = -1.;
+        alpha=acos(tmp);
+      }
+
+
+      epsilon_alpha=0.0;
+      epsilon_beta=0.024;
+      epsilon_gamma=0.024;
+      
+      beta = beta+epsilon_beta*sin(beta);
+      gamma = gamma+epsilon_gamma*sin(gamma);
+      
+      if ( mass1 >= mass2 )
       {
         q = mass2 / mass1;
       }
@@ -164,12 +201,39 @@ INT4 XLALSimIMREOBFinalMassSpinPrec(
       eISCO = sqrt( 1. - 2./(3.*rISCO) );
       *finalMass = 1. - ( (1.-eISCO)*eta
                  + 16.*eta*eta*( 0.00258 - 0.0773/(1./((1.+1/q/q)/(1.+1/q)/(1.+1/q))*atl-1.6939) - 0.25*(1.-eISCO)) );
-      *finalSpin = tmpVar + tmpVar*eta*( s9*eta*tmpVar*tmpVar + s8*eta*eta*tmpVar + s7*eta*tmpVar
-                 + s6*tmpVar*tmpVar + s4v2*tmpVar + s5v2*eta + t0v2)
-                 + eta*( 2.*sqrt(3.) + t2v2*eta + t3v2 *eta*eta );
+      q=1./q;
+      csi = 0.474046;
+      a_tot_prec=(chi1*cos(beta)+chi2*cos(gamma)*q*q)/((1.0+q)*(1.0+q));
+      aeff=a_tot_prec+csi*eta*(chi1*cos(beta)+chi2*cos(gamma));
+      rISCO = XLALSimRadiusKerrISCO( aeff );
+      eISCO = XLALSimEnergyKerrISCO( rISCO );
+      LISCO = XLALSimAngMomKerrISCO( rISCO );
+      /* The k00, k01, ... coefficients are defined in LALSimBlackHoleRingdown.h */
+      REAL8 eta4 = eta * eta3;
 
+      if (fabs(aeff) > 0.) {
+	REAL8 aeff2 = aeff * aeff;
+  REAL8 aeff3 = aeff * aeff2;
+  REAL8 aeff4 = aeff * aeff3;
+  fitpart = k00 * eta + k01 * aeff * eta + k02 * aeff2 * eta  + k03 * aeff3 * eta + k04 * aeff4 * eta +
+    k10 * eta2 + k11 * aeff * eta2 + k12 * aeff2 * eta2 + k13 * aeff3 * eta2 +  k14 * aeff4 * eta2 +
+    k20 * eta3 + k21 * aeff * eta3 + k22 * aeff2 * eta3 +  k23 * aeff3 * eta3 + k24 * aeff4 * eta3 +
+    k30 * eta4 +  k31 * aeff * eta4 + k32 * aeff2 * eta4 + k33 * aeff3 * eta4 + k34 * aeff4 * eta4;
 
 
+      } else {
+        fitpart = k00 * eta + k10 * eta2 + k20 * eta3 + k30 * eta4;
+
+      }
+      REAL8 ell_norm=fabs(LISCO-2*a_tot_prec*(eISCO-1.0)+fitpart);
+      REAL8 prefactor=1./((1+q)*(1+q));
+      *finalSpin=prefactor*sqrt(pow(chi1,2)+pow(chi2,2)*pow(q,4)+2*chi1*chi2*pow(q,2)*cos(alpha)+
+				 2*(chi1*cos(beta)+chi2*pow(q,2)*cos(gamma))*ell_norm*q+pow(ell_norm,2)*pow(q,2));
+
+			if (*finalSpin>1){
+				*finalSpin=1-1e-6;
+			}
+
       /* Stas: Below is original Rez.-Bar. implementation,
       Above I use the aligned case fit because spin1 and spin2
       currently have only z-components equal to projection of the spins on

lalsimulation/src/LALSimIMR.h

@@ -148,7 +148,8 @@ int XLALSimIMRSpinEOBWaveformAll(
                                  const REAL8     INspin2x,
                                  const REAL8     INspin2y,
                                  const REAL8     INspin2z,
-                                 const UINT4     PrecEOBversion
+                                 const UINT4     PrecEOBversion,
+                                 const UINT4     flagEulerextension
                                  );
 
 /* in module LALSimIMREOBNRv2HMROM.c */

lalsimulation/src/LALSimIMRCalculateSpinPrecEOBHCoeffs.c

 #ifndef _LALSimIMRCalculateSpinPrecEOBHCoeffs_C
 #define _LALSimIMRCalculateSpinPrecEOBHCoeffs_C
+#include "LALSimIMRSpinEOBHamiltonian.h"
 
 /**
  * \author Craig Robinson, Yi Pan, Stas Babak, Prayush Kumar, Andrea Taracchini
@@ -101,6 +102,19 @@ static int XLALSimIMRCalculateSpinPrecEOBHCoeffs(
   coeffs->KK = KK = ifthenelse(1.5-SpinAlignedEOBversion,
                                c0 + c1*eta + c2*eta*eta,
                                c20 + c21*eta + c22*(eta*eta) + c23*(eta*eta)*eta);
+    REAL8 chi = a / (1. - 2. * eta);
+    REAL8 eta2 = eta * eta, eta3 = eta2 * eta;
+    REAL8 chi2 = chi * chi, chi3 = chi2 * chi;
+    if (SpinAlignedEOBversion == 4)
+    {
+        coeffs->KK = KK =
+        coeff00K + coeff01K * chi + coeff02K * chi2 + coeff03K * chi3 +
+        coeff10K * eta + coeff11K * eta * chi + coeff12K * eta * chi2 +
+        coeff13K * eta * chi3 + coeff20K * eta2 + coeff21K * eta2 * chi +
+        coeff22K * eta2 * chi2 + coeff23K * eta2 * chi3 + coeff30K * eta3 +
+        coeff31K * eta3 * chi + coeff32K * eta3 * chi2 + coeff33K * eta3 * chi3;
+        //      printf("KK %.16e\n", KK);
+    }
   m1PlusEtaKK = -1. + eta*KK;
   const REAL8 invm1PlusEtaKK = 1./m1PlusEtaKK;
   /* Eqs. 5.77 - 5.81 of PRD 81, 084024 (2010) */
@@ -122,13 +136,42 @@ static int XLALSimIMRCalculateSpinPrecEOBHCoeffs(
                      + ((k1p2*k1p2)-4.*(k1p2*k2)+2.*k2*k2+4.*k1*k3-4.*k4)*0.5*invm1PlusEtaKK+(256./5.)*ln2)
                     );
   coeffs->k5l= k5l= ifthenelsezero(SpinAlignedEOBversion-1.5, (m1PlusEtaKK*m1PlusEtaKK) * (64./5.));
+  if ( SpinAlignedEOBversion == 4 ) {
+      coeffs->k5 = k5 =  m1PlusEtaKK*m1PlusEtaKK
+                                       * (-4237./60.+128./5.*LAL_GAMMA+2275.*LAL_PI*LAL_PI/512.
+                                          - third*(a*a)*(k1p3-3.*(k1*k2)+3.*k3)
+                                          - ((k1p3*k1p2)-5.*(k1p3*k2)+5.*k1*k2*k2+5.*k1p2*k3-5.*k2*k3-5.*k1*k4)*fifth*invm1PlusEtaKK*invm1PlusEtaKK
+                                          + ((k1p2*k1p2)-4.*(k1p2*k2)+2.*k2*k2+4.*k1*k3-4.*k4)*0.5*invm1PlusEtaKK+(256./5.)*ln2  + (41. * LAL_PI * LAL_PI / 32. -
+                                                                                                                                    221. / 6.) * eta );
+      coeffs->k5l= k5l= (m1PlusEtaKK*m1PlusEtaKK) * (64./5.);
+   }
+
 
   /* Now calibrated parameters for spin models */
   coeffs->d1 = ifthenelsezero(1.5-SpinAlignedEOBversion, -69.5);
   coeffs->d1v2 = ifthenelsezero(SpinAlignedEOBversion-1.5, -74.71 - 156.*eta + 627.5*eta*eta);
   coeffs->dheffSS = ifthenelsezero(1.5-SpinAlignedEOBversion, 2.75);
   coeffs->dheffSSv2 = ifthenelsezero(SpinAlignedEOBversion-1.5, 8.127 - 154.2*eta + 830.8*eta*eta);
-
+  if (SpinAlignedEOBversion==4) {
+      coeffs->d1 = 0.;
+      coeffs->dheffSS = 0.;
+      // dSO
+      coeffs->d1v2 =
+      coeff00dSO + coeff01dSO * chi + coeff02dSO * chi2 + coeff03dSO * chi3 +
+      coeff10dSO * eta + coeff11dSO * eta * chi + coeff12dSO * eta * chi2 +
+      coeff13dSO * eta * chi3 + coeff20dSO * eta2 + coeff21dSO * eta2 * chi +
+      coeff22dSO * eta2 * chi2 + coeff23dSO * eta2 * chi3 + coeff30dSO * eta3 +
+      coeff31dSO * eta3 * chi + coeff32dSO * eta3 * chi2 + coeff33dSO * eta3 * chi3;
+      
+      // dSS
+      coeffs->dheffSSv2 =
+      coeff00dSS + coeff01dSS * chi + coeff02dSS * chi2 + coeff03dSS * chi3 +
+      coeff10dSS * eta + coeff11dSS * eta * chi + coeff12dSS * eta * chi2 +
+      coeff13dSS * eta * chi3 + coeff20dSS * eta2 + coeff21dSS * eta2 * chi +
+      coeff22dSS * eta2 * chi2 + coeff23dSS * eta2 * chi3 + coeff30dSS * eta3 +
+      coeff31dSS * eta3 * chi + coeff32dSS * eta3 * chi2 + coeff33dSS * eta3 * chi3;
+      //          printf("dSO %.16e, dSS %.16e\n", coeffs->d1v2,coeffs->dheffSSv2);
+  }
   return XLAL_SUCCESS;
 }
 

lalsimulation/src/LALSimIMRSpinEOBFactorizedWaveformCoefficientsPrec.c

@@ -117,11 +117,11 @@ XLALSimIMREOBCalcSpinPrecFacWaveformCoefficients(
 //        double chi = (chiS + chiA * dM / (1. - 2. * eta));
 //        coeffs->delta22v6 = eta*(1./4.*(1. - 1080.*chi - sqrt((1. - 1080.*chi)*(1. - 1080*chi) + 8.*(13.5 +270.*chi +13.5*chi*chi))));
         }
-    if (SpinAlignedEOBversion == 3 /*&& chiS + chiA * dM / (1. - 2. * eta) > 0.*/) {
-//        		coeffs->delta22v6 = -540. * eta * (chiS + chiA * dM / (1. - 2. * eta));
-        double chi = (chiS + chiA * dM / (1. - 2. * eta));
-        coeffs->delta22v6 = eta*(1./4.*(1. - 1080.*chi - sqrt((1. - 1080.*chi)*(1. - 1080*chi) + 8.*(13.5 +270.*chi +13.5*chi*chi))));
-    }
+//    if (SpinAlignedEOBversion == 3 /*&& chiS + chiA * dM / (1. - 2. * eta) > 0.*/) {
+////        		coeffs->delta22v6 = -540. * eta * (chiS + chiA * dM / (1. - 2. * eta));
+//        double chi = (chiS + chiA * dM / (1. - 2. * eta));
+//        coeffs->delta22v6 = eta*(1./4.*(1. - 1080.*chi - sqrt((1. - 1080.*chi)*(1. - 1080*chi) + 8.*(13.5 +270.*chi +13.5*chi*chi))));
+//    }
 	coeffs->rho22v2 = -43. / 42. + (55. * eta) / 84.;
 	coeffs->rho22v3 = (-2. * (chiS + chiA * dM - chiS * eta)) / 3.;
 	switch (SpinAlignedEOBversion) {
@@ -143,11 +143,25 @@ XLALSimIMREOBCalcSpinPrecFacWaveformCoefficients(
 		break;
 	}
 	coeffs->rho22v5 = (-34. * a) / 21.;
+    if (SpinAlignedEOBversion == 3)
+    {
+        coeffs->rho22v5 =
+        (-34. / 21. + 49. * eta / 18. + 209. * eta2 / 126.) * chiS +
+        (-34. / 21. - 19. * eta / 42.) * dM * chiA;
+    }
 	coeffs->rho22v6 = 1556919113. / 122245200. + (89. * a2) / 252. - (48993925. * eta) / 9779616.
 		- (6292061. * eta2) / 3259872. + (10620745. * eta3) / 39118464.
 		+ (41. * eta * LAL_PI * LAL_PI) / 192.;
 	coeffs->rho22v6l = -428. / 105.;
 	coeffs->rho22v7 = (18733. * a) / 15876. + a * a2 / 3.;
+    if (SpinAlignedEOBversion == 3)
+    {
+        coeffs->rho22v7 =
+        a3 / 3. + chiA * dM * (18733. / 15876. + (50140. * eta) / 3969. +
+                               (97865. * eta2) / 63504.) +
+        chiS * (18733. / 15876. + (74749. * eta) / 5292. -
+                (245717. * eta2) / 63504. + (50803. * eta3) / 63504.);
+    }
 	switch (SpinAlignedEOBversion) {
 	case 1:
 		coeffs->rho22v8 = eta * (-5.6 - 117.6 * eta) - 387216563023. / 160190110080. + (18353. * a2) / 21168. - a2 * a2 / 8.;
@@ -291,8 +305,7 @@ XLALSimIMREOBCalcSpinPrecFacWaveformCoefficients(
 	coeffs->delta32vh9 = -9112. / 405. + (208. * LAL_PI * LAL_PI) / 63.;
 
 	coeffs->rho32v = (4. * chiS * eta) / (-3. * m1Plus3eta);
-    /** TODO The term proportional to eta^2 a^2 in coeffs->rho32v2 is wrong, but it was used in the calibration of SEOBNRv2 */
-	coeffs->rho32v2 = (-4. * a2 * eta2) / (9. * m1Plus3eta2) + (328. - 1115. * eta
+	coeffs->rho32v2 = (328. - 1115. * eta
 				       + 320. * eta2) / (270. * m1Plus3eta);
 	//coeffs->rho32v3 = (2. * (45. * a * m1Plus3eta3
 	//     - a * eta * (328. - 2099. * eta + 5. * (733. + 20. * a2) * eta2

lalsimulation/src/LALSimIMRSpinEOBHamiltonianPrec.c

@@ -1415,6 +1415,7 @@ static int XLALSpinPrecHcapRvecDerivative(
        tplspin = 0.0;
        break;
      case 2:
+     case 4:
        tplspin = (1.-2.*eta) * chiS + (mass1 - mass2)/(mass1 + mass2) * chiA;
        break;
      default:

lalsimulation/src/LALSimIMRSpinEOBHcapNumericalDerivativePrec.c

@@ -468,6 +468,7 @@ static REAL8 XLALSpinPrecHcapNumDerivWRTParam(
 		tplspin = 0.0;
 		break;
 	case 2:
+    case 4:
         /* See below Eq. 4 of PRD 89, 061502(R) (2014)*/
 		tplspin = (1. - 2. * eta) * chiS + (mass1 - mass2) / (mass1 + mass2) * chiA;
         break;

lalsimulation/src/LALSimNeutronStarEOSTabular.c

@@ -363,51 +363,11 @@ LALSimNeutronStarEOS *XLALSimNeutronStarEOSFromFile(const char *fname)
  * @brief Creates an equation of state structure from tabulated equation
  * of state data of a known name.