Commit a80f4b83 authored by Reinhard Prix's avatar Reinhard Prix

XLALAddBinaryTimes: fixed root-finding accuracy to a sensible requirement

- require maximal CW-phase error of 1e-3
- fixes #1357
Original: df54e62c43ff8b67bb899c8121eb5d3586b1419f
parent 6f322e6a
......@@ -30,7 +30,6 @@
#include <gsl/gsl_roots.h>
/*---------- local DEFINES ----------*/
#define EA_ACC 1E-9 /* the timing accuracy of LALGetBinaryTimes in seconds */
/*----- Macros ----- */
......@@ -297,6 +296,7 @@ XLALAddBinaryTimes ( SSBtimes **tSSBOut, //!< [out] reference-time offsets in
REAL8 sinw = sin ( argp ); /* the sin and cos of the argument of periapsis */
REAL8 cosw = cos ( argp );
REAL8 n = LAL_TWOPI / Porb;
REAL8 Freq = Doppler->fkdot[0];
REAL8 refTimeREAL8 = XLALGPSGetREAL8 ( &tSSBIn->refTime );
......@@ -305,7 +305,10 @@ XLALAddBinaryTimes ( SSBtimes **tSSBOut, //!< [out] reference-time offsets in
REAL8 B = n * asini * sinw; // see Eq.(eq:defB)
REAL8 tp = XLALGPSGetREAL8 ( &(Doppler->tp) );
REAL8 Tp = tp + asini * sinw * (1 - e); // see Eq.(eq:defTPeri)
REAL8 acc = LAL_TWOPI * EA_ACC / Porb; // root-finding accuracy, EA_ACC represents the required timing precision in seconds
REAL8 maxPhaseErr = 1e-3; // maximal allowed CW phase-error due to error in E
REAL8 epsabs = maxPhaseErr / ( Freq * Porb); // absolute root-finding accuracy required on E
REAL8 epsrel = 0; // no constraint on relative accuracy
/* loop over the SFTs i */
for ( UINT4 i = 0; i < numSteps; i++ )
......@@ -340,14 +343,14 @@ XLALAddBinaryTimes ( SSBtimes **tSSBOut, //!< [out] reference-time offsets in
E_i = gsl_root_fsolver_root(s);
E_lo = gsl_root_fsolver_x_lower (s);
E_hi = gsl_root_fsolver_x_upper (s);
status = gsl_root_test_interval ( E_lo, E_hi, acc, 0 );
status = gsl_root_test_interval ( E_lo, E_hi, epsabs, epsrel );
if (status == GSL_SUCCESS) { XLALPrintInfo ("Converged:\n"); }
XLALPrintInfo ("%5d [%.7f, %.7f] %.7f %+10.7g %10.7g\n", iter, E_lo, E_hi, E_i, acc, E_hi - E_lo);
XLALPrintInfo ("%5d [%.7f, %.7f] %.7f %+10.7g %10.7g\n", iter, E_lo, E_hi, E_i, epsabs, E_hi - E_lo);
} while ( (status == GSL_CONTINUE) && (iter < max_iter) );
XLAL_CHECK ( status == GSL_SUCCESS, XLAL_EMAXITER, "Eccentric anomaly: failed to converge within %d iterations\n", max_iter );
XLAL_CHECK ( status == GSL_SUCCESS, XLAL_EMAXITER, "Eccentric anomaly: failed to converge to epsabs=%g within %d iterations\n", epsabs, max_iter );
gsl_root_fsolver_free(s);
} // gsl-root finding block
......
......@@ -183,6 +183,7 @@ main ( int argc, char *argv[] )
MultiSSBtimes *multiBinary_test = NULL;
PulsarDopplerParams doppler;
memset(&doppler, 0, sizeof(doppler));
doppler.fkdot[0] = 100;
doppler.tp = orbit.tp;
doppler.argp = orbit.argp;
doppler.asini = orbit.asini;
......
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