diff --git a/src/tides/tidal_lambda2.c b/src/tides/tidal_lambda2.c index ee5eb839a3f46839f886fe8d0a48b824a5791cfd..8e7518805e78d49f3f341d39a21065326efde606 100644 --- a/src/tides/tidal_lambda2.c +++ b/src/tides/tidal_lambda2.c @@ -116,31 +116,6 @@ double tidal_lambda2(struct stardata_t * const stardata, const int itmax = 200; const Boolean uselog = FALSE; const double convergence_alpha = 1.0; - - if(0) - { - /* - * Loop to crudely find where the root is near - */ - double bestx = -666.0; - double bestf = 1e200; - for(double x = min; - x <= max; - x += (max-min)*0.01) - { - const double f = Zahn_Eq_14b(x,RHS); - if(1)printf("Func(%g) = %g\n", - x, - f - ); - if(fabs(f) < bestf) - { - bestx = x; - bestf = fabs(f); - } - } - printf("Root is near x = %g\n",bestx); - } const double xa = generic_bisect(&err, BISECT_USE_MONOCHECKS_FATAL, BISECTOR_TIDES_XA, // make negative for logging @@ -159,7 +134,7 @@ double tidal_lambda2(struct stardata_t * const stardata, Dprint("Check xa falls into CZ with xa>xb xa = %g, xb %g\n", xa,xb); - if(xa>xb) + if(xa > xb) { const double integral1 = GSL_integrator(xa, // lower 1.0, // upper @@ -186,18 +161,18 @@ double tidal_lambda2(struct stardata_t * const stardata, } else { - /* - * There is a root but if falls below convective zone, single integral in Eq. 13 - */ - integral = GSL_integrator(xb, // lower - 1.0, // upper - integration_tolerance, - &integration_error, - GSL_INTEGRATOR_QAG, - &Zahn_integral_13); - Dprint("irrelevant root ignored: xb = %g, integral = %g\n", - xb, - integral); + /* + * There is a root but if falls below convective zone, single integral in Eq. 13 + */ + integral = GSL_integrator(xb, // lower + 1.0, // upper + integration_tolerance, + &integration_error, + GSL_INTEGRATOR_QAG, + &Zahn_integral_13); + Dprint("irrelevant root ignored: xb = %g, integral = %g\n", + xb, + integral); } } else @@ -223,8 +198,7 @@ double tidal_lambda2(struct stardata_t * const stardata, { lambda = 0.0; } -// const double Eq15 = 0.019 * pow(alpha,4.0/3.0) * sqrt(320.0 / (320.0 + Pow2(Pi / (2.0 * tf)))); -// ratio above is the wrong way + const double Eq15 = 0.019 * pow(alpha,4.0/3.0) * sqrt(320.0 / (320.0 + Pow2((2.0 * tf) / Pi))); Dprint("lambda = %g, approximate Eq15 gives %g\n",