Skip to content
Snippets Groups Projects
Commit 983f152a authored by Izzard, Robert Dr (Maths & Physics)'s avatar Izzard, Robert Dr (Maths & Physics)
Browse files

fix indentation, removed non-bisecting root find (that wasn't used)

parent 985f127c
No related branches found
No related tags found
No related merge requests found
......@@ -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",
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment