diff --git a/src/MINT/MINT_load_MS_grid.c b/src/MINT/MINT_load_MS_grid.c
index 17102a3a904f2d786a54f2fac84a293a060e09e2..0f0dbc8422dcc8344b64a384b79f3166cb6b16ed 100644
--- a/src/MINT/MINT_load_MS_grid.c
+++ b/src/MINT/MINT_load_MS_grid.c
@@ -29,7 +29,7 @@ Boolean MINT_load_MS_grid(struct stardata_t * const stardata)
{ MINT_MS_DATA_ITEMS_LIST };
#undef X
- return MINT_Load_Table(
+ Boolean ret = MINT_Load_Table(
stardata,
MINT_TABLE_MS,
NULL,
@@ -39,6 +39,8 @@ Boolean MINT_load_MS_grid(struct stardata_t * const stardata)
MINT_MS_parameter_actions,
MINT_MS_data_actions,
vb);
+
+ return ret;
}
#endif // MINT
diff --git a/src/MINT/MINT_mix_stars.c b/src/MINT/MINT_mix_stars.c
index 116360b8ae4fa9064971e4be4e0421b0d878c032..7d1ed2dd6154f32591c1d4a5cc962df683fccdaa 100644
--- a/src/MINT/MINT_mix_stars.c
+++ b/src/MINT/MINT_mix_stars.c
@@ -45,7 +45,8 @@ void MINT_mix_stars(struct stardata_t * const stardata)
*/
merged_star->angular_momentum =
stardata->preferences->merger_angular_momentum_factor *
- breakup_angular_momentum(merged_star);
+ breakup_angular_momentum(stardata,
+ merged_star);
calculate_rotation_variables(merged_star,
merged_star->radius);
diff --git a/src/MINT/MINT_stellar_structure_CHeB.c b/src/MINT/MINT_stellar_structure_CHeB.c
index 859b7596a9fafb60b7ccfe8d3353e9ecc01430d6..e17eb710f758fdf8259b9bea34069f9c14fa77dd 100644
--- a/src/MINT/MINT_stellar_structure_CHeB.c
+++ b/src/MINT/MINT_stellar_structure_CHeB.c
@@ -282,12 +282,12 @@ Stellar_type MINT_stellar_structure_CHeB(struct stardata_t * RESTRICT const star
newstar->age,
newstar->epoch);
- newstar->k2 = get_data(MINT_CHeB_K2);
+ //newstar->moment_of_inertia_factor = get_data(MINT_CHeB_MOMENT_OF_INERTIA_FACTOR);
/*
- * The line above should be replaced with
- * newstar->k2 = get_data(MINT_CHeB_MOMENT_OF_INERTIA_FACTOR);
+ * The line above should be replaced with
+ * newstar->moment_of_inertia_factor = get_data(MINT_CHeB_MOMENT_OF_INERTIA_FACTOR);
* once that quantity is available in CHeB grids
- * The K2 above is MESA's and not I/(MR^2)
+ * The MOMENT_OF_INERTIA_FACTOR above is MESA's and not I/(MR^2)
*/
newstar->E2 =
available[MINT_CHeB_TIDAL_E2] ?
diff --git a/src/MINT/MINT_stellar_structure_MS.c b/src/MINT/MINT_stellar_structure_MS.c
index 363a033959980c330ce3c66ff8bd1fd81badd59e..a9d1e540038e000631cc7d70e91de177dd9d52a5 100644
--- a/src/MINT/MINT_stellar_structure_MS.c
+++ b/src/MINT/MINT_stellar_structure_MS.c
@@ -293,10 +293,10 @@ Stellar_type MINT_stellar_structure_MS(struct stardata_t * RESTRICT const starda
newstar->starnum,
stardata->model.dt);
- newstar->k2 = get_data(MINT_MS_MOMENT_OF_INERTIA_FACTOR);
- if(vb)printf("Moment of inertia factor k2=I/(MR^2) : star = %d, k22 = %g\n",
+ newstar->moment_of_inertia_factor = get_data(MINT_MS_MOMENT_OF_INERTIA_FACTOR);
+ if(vb)printf("Moment of inertia factor moment_of_inertia_factor=I/(MR^2) : star = %d, moment_of_inertia_factor2 = %g\n",
newstar->starnum,
- stardata->k2);
+ newstar->moment_of_inertia_factor);
newstar->E2 =
available[MINT_MS_TIDAL_E2] ?
diff --git a/src/RLOF/RLOF_dynamical_comenv.c b/src/RLOF/RLOF_dynamical_comenv.c
index b173cc6f0abbf9532b07478eb289c083a56cac7d..c4921f75a180cf656b0fa459276269599cc43efb 100644
--- a/src/RLOF/RLOF_dynamical_comenv.c
+++ b/src/RLOF/RLOF_dynamical_comenv.c
@@ -58,7 +58,7 @@ Boolean RLOF_dynamical_comenv(struct stardata_t * const RESTRICT stardata,
if(Is_zero(rc)) rc = star->radius;
- star->omega = star->angular_momentum/(K3_CORE_GYRATION*Pow2(rc)*star->mass);
+ star->omega = star->angular_momentum/(CORE_MOMENT_OF_INERTIA_FACTOR*Pow2(rc)*star->mass);
}
}
#endif//BSE
diff --git a/src/binary_c_macros.h b/src/binary_c_macros.h
index 74c3b51669ae877d68665463eda0f722a84b091f..6f7b331a42ee0b31a82f8c65e00b9218dd84efc4 100644
--- a/src/binary_c_macros.h
+++ b/src/binary_c_macros.h
@@ -138,7 +138,7 @@
#define SPECIFIC_ANGULAR_MOMENTUM_CGS (R_SUN*R_SUN/YEAR_LENGTH_IN_SECONDS)
#define K2_GYRATION 0.1
-#define K3_CORE_GYRATION 0.21
+#define CORE_MOMENT_OF_INERTIA_FACTOR 0.21
#define MR23YR 0.4311
#define MINIMUM_STELLAR_ANGMOM 1e-50
#define MAXIMUM_STELLAR_ANGMOM 1e100
diff --git a/src/binary_c_parameters.h b/src/binary_c_parameters.h
index dcf45b9e92efd776ae8ebf2d805b5bfc06621f3a..fdc59fc2100b69154d1cae002b9a8b0d52b25985 100644
--- a/src/binary_c_parameters.h
+++ b/src/binary_c_parameters.h
@@ -288,10 +288,10 @@
//#define HE_STAR_ACCRETING_H_FIX
/*
- * The DTR_FUDGES are a set of fudges to fix this function
+ * The DTR_CHECKS are a set of fudges to fix this function
* so time_remaining>0 at all times. Otherwise, things gets really slow
*/
-#define DTR_FUDGES
+#define DTR_CHECKS
/* choose one of these for WD cooling */
#define MESTEL_COOLING
@@ -402,15 +402,6 @@
*/
//#define REJUVENATE_HERZTSPRUNG_GAP_STARS
-/*
- * Fudge such that K2 does not change too fast for non-compact stars.
- * This prevents it changing by > K2FUDGE_MAX per timestep.
- * Enable K2FUDGE_WARNING to output a warning when this happens.
- */
-#define K2FUDGE
-#define K2FUDGE_MAX 0.05
-//#define K2FUDGE_WARNING
-
/*
* minimum luminosity for TPAGB stars (fudge/check!)
*/
diff --git a/src/binary_c_structures.h b/src/binary_c_structures.h
index 7b86de4e46d034946306a38700b5db3f9dabb0c7..ee4c6922caad67dc291b81432c0497d22a31b8d1 100644
--- a/src/binary_c_structures.h
+++ b/src/binary_c_structures.h
@@ -1372,6 +1372,8 @@ struct preferences_t {
#endif // MINT
Boolean load_data_grids;
+
+ Boolean apply_Darwin_Radau_correction;
};
#ifdef EVOLUTION_SPLITTING
@@ -1695,7 +1697,6 @@ struct star_t {
double adaptive2_Rexcess_high;
Boolean adaptive2_converged;
#endif//ADAPTIVE_RLOF2
- /* RLOF interpolation variables */
double rol0; /* roche lobe radius at the start of RLOF */
double aj0; /* age at the start of RLOF */
double mass_at_RLOF_start; /* used in RLOF interpolation */
@@ -1706,7 +1707,9 @@ struct star_t {
double tm,tn; /* Main sequence and nuclear timescales */
Time tkh; /** Kelvin-Helmholtz timescale **/
double max_EAGB_He_core_mass;
- double menv,renv,k2; // envelope mass, radius, gyration constant,tidal constant
+ double menv; /* convective envelope mass */
+ double renv; /* convective envelope radius */
+ double moment_of_inertia_factor; /* I/(MR^2) [== k2 in BSE] */
double E2; /* set in grid data */
double E_Zahn; /* set in grid data, used to calculate lambda2_Zahn */
double compact_core_mass; /* mass in a compact core */
diff --git a/src/binary_star_functions/calculate_spins.c b/src/binary_star_functions/calculate_spins.c
index 885c80aa6b7fbf63b4de2edd0a0299b20dc89f26..e52ce9221dae5ac54e08f974c0cbf3f6369810da 100644
--- a/src/binary_star_functions/calculate_spins.c
+++ b/src/binary_star_functions/calculate_spins.c
@@ -3,11 +3,9 @@ No_empty_translation_unit_warning;
void calculate_spins(struct stardata_t * RESTRICT const stardata)
{
- Star_number k;
- Evolving_Starloop(k)
+ Foreach_evolving_star(star)
{
- SETstar(k);
- star->omega = omegaspin(star);
+ star->omega = omegaspin(stardata,star);
Dprint("Spin star %d (type %d) from j=%30.12e r=%g menv = %g - %g = %g I = %30.12e omega = %30.12e (omega orb %30.12e)\n",
star->starnum,
star->stellar_type,
@@ -16,7 +14,7 @@ void calculate_spins(struct stardata_t * RESTRICT const stardata)
star->mass,
Outermost_core_mass(star),
star->mass - Outermost_core_mass(star),
- moment_of_inertia(star,star->effective_radius),
+ moment_of_inertia(stardata,star,star->effective_radius),
star->omega,
stardata->common.orbit.angular_frequency
);
diff --git a/src/binary_star_functions/contact_system.c b/src/binary_star_functions/contact_system.c
index 1de912cd9ab82fc22db47720fc7186ad3ca8d356..826ce777279e08416aa71370cc428ef3ec1cb71e 100644
--- a/src/binary_star_functions/contact_system.c
+++ b/src/binary_star_functions/contact_system.c
@@ -194,7 +194,7 @@ static void contact_system_call_corerd(const Star_number i,
star->mass,
star->phase_start_core_mass,
stardata->common.metallicity_parameters[i]);
- star->omega = star->angular_momentum/(K3_CORE_GYRATION*rc*rc*star->mass);
+ star->omega = star->angular_momentum/(CORE_MOMENT_OF_INERTIA_FACTOR*rc*rc*star->mass);
}
}
#endif//BSE
diff --git a/src/binary_star_functions/prevent_overspin.c b/src/binary_star_functions/prevent_overspin.c
index c1e59e192fed5eb09ed32877e4d80eb49bb931ae..2c4f32e9d67016feb5106c992c3029031abfb995 100644
--- a/src/binary_star_functions/prevent_overspin.c
+++ b/src/binary_star_functions/prevent_overspin.c
@@ -55,7 +55,8 @@ Boolean prevent_overspin(struct star_t * const star,
Is_zero(star->omega)
) return FALSE;
- double jcrit = breakup_angular_momentum(star);
+ double jcrit = breakup_angular_momentum(stardata,
+ star);
/*
* Calculate the angular momentum the star would
diff --git a/src/common_envelope/common_envelope_event_handler.c b/src/common_envelope/common_envelope_event_handler.c
index f58f3ee12fc9c05c2ab4dab45362cde566b475c7..9066f503da8d96a67fca9d9c251a785d419f94ae 100644
--- a/src/common_envelope/common_envelope_event_handler.c
+++ b/src/common_envelope/common_envelope_event_handler.c
@@ -69,7 +69,7 @@ Event_handler_function common_envelope_event_handler(void * eventp,
if(Is_zero(rc)) rc = star->radius;
- star->omega = star->angular_momentum/(K3_CORE_GYRATION*Pow2(rc)*star->mass);
+ star->omega = star->angular_momentum/(CORE_MOMENT_OF_INERTIA_FACTOR*Pow2(rc)*star->mass);
}
}
#endif//BSE
diff --git a/src/common_envelope/common_envelope_evolution.h b/src/common_envelope/common_envelope_evolution.h
index c5080c90bafdfee08b49f9b75bb608a1bc233ec7..f5dec055affc914c0bff6997fd5823820b685f3c 100644
--- a/src/common_envelope/common_envelope_evolution.h
+++ b/src/common_envelope/common_envelope_evolution.h
@@ -140,8 +140,8 @@ No_empty_translation_unit_warning;
&Outermost_core_mass(accretor), \
&rc1, \
&rc2, \
- &k21, \
- &k22, \
+ &moment_of_inertia_factor1, \
+ &moment_of_inertia_factor2, \
&tn, \
common_envelope_mass, \
&donor->stellar_type, \
@@ -268,7 +268,7 @@ No_empty_translation_unit_warning;
if(stardata->model.coalesce>0) \
{ \
Append_logstring(LOG_COMENV_PRE, \
- "Merge%u [rc1=%g%crl1=%g or R2=%g%crl2=%g] ", \
+ "Merge%u [rc1=%g%crl1=%g or R2=%g%crl2=%g] ", \
stardata->model.coalesce, \
rc1, \
(stardata->model.coalesce==1 ? '>' : '<'), \
@@ -384,7 +384,7 @@ No_empty_translation_unit_warning;
r1 = donor->radius; \
rc1 = donor->core_radius; \
menv = donor->menv; \
- k21 = donor->k2; \
+ moment_of_inertia_factor1 = donor->moment_of_inertia_factor; \
memcpy(&tscls1[0], \
donor->timescales, \
sizeof(double)*TSCLS_ARRAY_SIZE); \
@@ -399,14 +399,14 @@ No_empty_translation_unit_warning;
* Debugging output for star 1
*/
-#define PRECALL1 printf("Precall &donor->age=%g m01=%g m1=%15.10e mc1=%15.10e (m1-mc1=%g) ntp=%g :: k21=%g l1=%g menv=%g r1=%15.10e rc1=%g renv=%g tm1=%g tn=%g st=%d\n", \
+#define PRECALL1 printf("Precall &donor->age=%g m01=%g m1=%15.10e mc1=%15.10e (m1-mc1=%g) ntp=%g :: moment_of_inertia_factor1=%g l1=%g menv=%g r1=%15.10e rc1=%g renv=%g tm1=%g tn=%g st=%d\n", \
donor->age, \
donor->phase_start_mass, \
donor->mass, \
Outermost_core_mass(donor), \
donor->mass-Outermost_core_mass(donor), \
stardata->star[donor->starnum].num_thermal_pulses, \
- k21, \
+ moment_of_inertia_factor1, \
l1, \
menv, \
r1, \
@@ -418,14 +418,14 @@ No_empty_translation_unit_warning;
);
-#define POSTCALL1 printf("Postcall &donor->age=%g m01=%g m1=%15.10e mc1=%15.10e (m1-mc1=%g) ntp=%g :: k21=%g l1=%g menv=%g r1=%15.10e rc1=%g renv=%g tm1=%g tn=%g st=%d\n", \
+#define POSTCALL1 printf("Postcall &donor->age=%g m01=%g m1=%15.10e mc1=%15.10e (m1-mc1=%g) ntp=%g :: moment_of_inertia_factor1=%g l1=%g menv=%g r1=%15.10e rc1=%g renv=%g tm1=%g tn=%g st=%d\n", \
donor->age, \
donor->phase_start_mass, \
donor->mass, \
Outermost_core_mass(donor), \
donor->mass-Outermost_core_mass(donor), \
stardata->star[donor->starnum].num_thermal_pulses, \
- k21, \
+ moment_of_inertia_factor1, \
l1, \
menv, \
r1, \
@@ -466,7 +466,7 @@ No_empty_translation_unit_warning;
rc1 = donor->core_radius; \
menv = donor->menv; \
/*renv = donor->renv;*/ \
- k21 = donor->k2; \
+ moment_of_inertia_factor1 = donor->moment_of_inertia_factor; \
memcpy(&tscls1[0], \
donor->timescales, \
sizeof(double)*TSCLS_ARRAY_SIZE); \
@@ -503,7 +503,7 @@ No_empty_translation_unit_warning;
r2 = accretor->radius; \
rc2 = accretor->core_radius; \
menv = accretor->menv; \
- k22 = accretor->k2; \
+ moment_of_inertia_factor2 = accretor->moment_of_inertia_factor; \
memcpy(&tscls2[0], \
accretor->timescales, \
sizeof(double) * TSCLS_ARRAY_SIZE); \
diff --git a/src/common_envelope/common_envelope_evolution_BSE.c b/src/common_envelope/common_envelope_evolution_BSE.c
index 690b110ccd84d2183f3e1698a569ec5d417a7215..7a6ca2f5182ec613657c24148054881ed807bd0e 100644
--- a/src/common_envelope/common_envelope_evolution_BSE.c
+++ b/src/common_envelope/common_envelope_evolution_BSE.c
@@ -116,7 +116,7 @@ int common_envelope_evolution_BSE (
double EbindI,EbindF=0.0,EorbI=0.0,EorbF=0.0,Ecirc,lambda1,lambda2;
double r1=0.0,r2=0.0,rc1=0.0,rc2=0.0,rl1=0.0,rl2=0.0,final_separation=0.0;
double Roche_separation=0.0;
- double tm1,tm2,tn,fage1,fage2,menv,mc22,mc3,mf,dM1,dM2,dm,k21,k22;
+ double tm1,tm2,tn,fage1,fage2,menv,mc22,mc3,mf,dM1,dM2,dm,moment_of_inertia_factor1,moment_of_inertia_factor2;
const double jorb_cgs = (donor->mass*M_SUN)*(accretor->mass*M_SUN)/
(M_SUN*(donor->mass+accretor->mass))*
@@ -286,7 +286,7 @@ int common_envelope_evolution_BSE (
r1 = donor->radius;
rc1 = donor->core_radius;
menv = donor->menv;
- k21 = donor->k2;
+ moment_of_inertia_factor1 = donor->moment_of_inertia_factor;
memcpy(&tscls1[0],
donor->timescales,
@@ -529,8 +529,8 @@ int common_envelope_evolution_BSE (
/*
* calculate stellar spins
*/
- ospin[0] = donor->angular_momentum / (k21*Pow2(r1)*((donor->mass)-(Outermost_core_mass(donor)))+K3_CORE_GYRATION*Pow2(rc1)* Outermost_core_mass(donor));
- ospin[1] = accretor->angular_momentum / (k22*Pow2(r2)*((accretor->mass)-(Outermost_core_mass(accretor)))+K3_CORE_GYRATION*Pow2(rc2)* Outermost_core_mass(accretor));
+ ospin[0] = donor->angular_momentum / (moment_of_inertia_factor1*Pow2(r1)*((donor->mass)-(Outermost_core_mass(donor)))+CORE_MOMENT_OF_INERTIA_FACTOR*Pow2(rc1)* Outermost_core_mass(donor));
+ ospin[1] = accretor->angular_momentum / (moment_of_inertia_factor2*Pow2(r2)*((accretor->mass)-(Outermost_core_mass(accretor)))+CORE_MOMENT_OF_INERTIA_FACTOR*Pow2(rc2)* Outermost_core_mass(accretor));
double ntp,ntp_since_mcmin;
@@ -651,7 +651,7 @@ int common_envelope_evolution_BSE (
sqrt(final_separation/(AU_IN_SOLAR_RADII*(Outermost_core_mass(donor) + accretor->mass)));
orbital_angular_frequency = TWOPI/orbital_period_years;
double dely = 0.50* accretor->mass *r2*r2*(ospin[1]*ospin[1] - orbital_angular_frequency*orbital_angular_frequency)/3.91e+08;
- dely *= k22;
+ dely *= moment_of_inertia_factor2;
EbindI = Max(0.0,EbindI - dely);
EorbF = EorbI + EbindI/alpha_ce;
final_separation = Outermost_core_mass(donor) * accretor->mass/(2.0*EorbF);
@@ -918,7 +918,7 @@ int common_envelope_evolution_BSE (
orbital_period_years = final_separation/AU_IN_SOLAR_RADII*sqrt(final_separation/(AU_IN_SOLAR_RADII*(Outermost_core_mass(donor) + accretor->mass)));
orbital_angular_frequency = TWOPI/orbital_period_years;
const double dely = 0.50* accretor->mass * Pow2(r2) * (Pow2(ospin[1]) - Pow2(orbital_angular_frequency))/3.91e+08
- * K3_CORE_GYRATION;
+ * CORE_MOMENT_OF_INERTIA_FACTOR;
EbindI = Max(0.0,EbindI - dely);
EorbF = EorbI + EbindI/alpha_ce;
final_separation = Outermost_core_mass(donor) * accretor->mass/(2.0*EorbF);
@@ -1496,11 +1496,11 @@ int common_envelope_evolution_BSE (
*/
donor->angular_momentum =
orbital_angular_frequency *
- (k21*Pow2(r1) * (donor->mass-Outermost_core_mass(donor))+
- K3_CORE_GYRATION *Pow2(rc1)* Outermost_core_mass(donor));
- Dprint("Set JPSIN1 (2) = omega_orb = %g * (k21=%g * (r1=%g)^2 * (m1=%g - mc1=%g) + (rc1=%g)^2 * mc1=%g = %g\n",
+ (moment_of_inertia_factor1*Pow2(r1) * (donor->mass-Outermost_core_mass(donor))+
+ CORE_MOMENT_OF_INERTIA_FACTOR *Pow2(rc1)* Outermost_core_mass(donor));
+ Dprint("Set JPSIN1 (2) = omega_orb = %g * (moment_of_inertia_factor1=%g * (r1=%g)^2 * (m1=%g - mc1=%g) + (rc1=%g)^2 * mc1=%g = %g\n",
orbital_angular_frequency,
- k21,
+ moment_of_inertia_factor1,
r1,
donor->mass,
Outermost_core_mass(donor),
@@ -1543,7 +1543,8 @@ int common_envelope_evolution_BSE (
*
* This is wrong if no material is ejected!
*/
- donor->angular_momentum = breakup_angular_momentum(donor);
+ donor->angular_momentum = breakup_angular_momentum(stardata,
+ donor);
}
else
{
@@ -1551,7 +1552,9 @@ int common_envelope_evolution_BSE (
"Unknown comenv merger spin method %d",
stardata->preferences->comenv_merger_spin_method);
}
- donor->omega = donor->angular_momentum / moment_of_inertia(donor,donor->radius);
+ donor->omega = donor->angular_momentum / moment_of_inertia(stardata,
+ donor,
+ donor->radius);
accretor->angular_momentum = 0.0;
if(merged_stellar_type == MASSLESS_REMNANT)
@@ -1646,8 +1649,8 @@ int common_envelope_evolution_BSE (
/*
* leave the original stellar spins and let tides do the rest later
*/
- donor->angular_momentum = ospin[0]*(k21*Pow2(r1)*(donor->mass-Outermost_core_mass(donor))+K3_CORE_GYRATION*Pow2(rc1)* Outermost_core_mass(donor));
- accretor->angular_momentum = ospin[1]*(k22*Pow2(r2)*(accretor->mass-Outermost_core_mass(accretor))+K3_CORE_GYRATION*Pow2(rc2)* Outermost_core_mass(accretor));
+ donor->angular_momentum = ospin[0]*(moment_of_inertia_factor1*Pow2(r1)*(donor->mass-Outermost_core_mass(donor))+CORE_MOMENT_OF_INERTIA_FACTOR*Pow2(rc1)* Outermost_core_mass(donor));
+ accretor->angular_momentum = ospin[1]*(moment_of_inertia_factor2*Pow2(r2)*(accretor->mass-Outermost_core_mass(accretor))+CORE_MOMENT_OF_INERTIA_FACTOR*Pow2(rc2)* Outermost_core_mass(accretor));
}
else if(stardata->preferences->comenv_ejection_spin_method ==
COMENV_EJECTION_SPIN_METHOD_SYNCHRONIZE)
@@ -1657,9 +1660,9 @@ int common_envelope_evolution_BSE (
* enough angular momentum to do this.
*/
donor->angular_momentum = orbital_angular_frequency*
- (k21*Pow2(r1)*(donor->mass-Outermost_core_mass(donor))+K3_CORE_GYRATION*Pow2(rc1)* Outermost_core_mass(donor));
+ (moment_of_inertia_factor1*Pow2(r1)*(donor->mass-Outermost_core_mass(donor))+CORE_MOMENT_OF_INERTIA_FACTOR*Pow2(rc1)* Outermost_core_mass(donor));
accretor->angular_momentum = orbital_angular_frequency*
- (k22*Pow2(r2)*(accretor->mass-Outermost_core_mass(accretor))+K3_CORE_GYRATION*Pow2(rc2)* Outermost_core_mass(accretor));
+ (moment_of_inertia_factor2*Pow2(r2)*(accretor->mass-Outermost_core_mass(accretor))+CORE_MOMENT_OF_INERTIA_FACTOR*Pow2(rc2)* Outermost_core_mass(accretor));
}
else
{
@@ -2498,16 +2501,16 @@ end_section:
/* conserve angular momentum to estimate the spin of the
* merged star */
- star->r_star_postvrot=jtot_cgs/((k21*(donor->mass-Outermost_core_mass(donor))*r1
- +K3_CORE_GYRATION*(Outermost_core_mass(donor)))
+ star->r_star_postvrot=jtot_cgs/((moment_of_inertia_factor1*(donor->mass-Outermost_core_mass(donor))*r1
+ +CORE_MOMENT_OF_INERTIA_FACTOR*(Outermost_core_mass(donor)))
*M_SUN*R_SUN) * 1e-5;
/* Vrot2 */
/* calculate with core angular momentum conserved in the cores */
jtot_cgs=Max(0.0,jtot_cgs-djcores); // should never be <0!
- star->r_star_postvrot2=jtot_cgs/((k21*(donor->mass-Outermost_core_mass(donor))*r1
- +K3_CORE_GYRATION*(Outermost_core_mass(donor)))
+ star->r_star_postvrot2=jtot_cgs/((moment_of_inertia_factor1*(donor->mass-Outermost_core_mass(donor))*r1
+ +CORE_MOMENT_OF_INERTIA_FACTOR*(Outermost_core_mass(donor)))
*M_SUN*R_SUN) * 1e-5;
/* Vrot3 */
@@ -2515,8 +2518,8 @@ end_section:
double dm=-(donor->mass-m1init-m2init);
double dJ=dm*star->r_star_postvcrit*r1*(1e5*M_SUN*R_SUN);
jtot_cgs=Max(0.0,jtot_cgs-dJ);//sometimes < 0!
- star->r_star_postvrot3=jtot_cgs/((k21*(donor->mass-Outermost_core_mass(donor))*r1
- +K3_CORE_GYRATION*(Outermost_core_mass(donor)))
+ star->r_star_postvrot3=jtot_cgs/((moment_of_inertia_factor1*(donor->mass-Outermost_core_mass(donor))*r1
+ +CORE_MOMENT_OF_INERTIA_FACTOR*(Outermost_core_mass(donor)))
*M_SUN*R_SUN) * 1e-5;
/* save the final spin rate (but is this used any more?) */
@@ -3626,8 +3629,8 @@ void comenv_accretion(struct stardata_t * stardata,
double *mc2,
double *rc1,
double *rc2,
- double *k21,
- double *k22,
+ double *moment_of_inertia_factor1,
+ double *moment_of_inertia_factor2,
double *tn,
double common_envelope_mass,
Stellar_type *stellar_type1,
@@ -3712,7 +3715,7 @@ void comenv_accretion(struct stardata_t * stardata,
rc1,
menv,
renv,
- k21,
+ moment_of_inertia_factor1,
&(stardata->star[donor->starnum]),
stardata,
STELLAR_STRUCTURE_CALLER_common_envelope_evolution);
@@ -3735,7 +3738,7 @@ void comenv_accretion(struct stardata_t * stardata,
rc2,
menv,
renv,
- k22,
+ moment_of_inertia_factor2,
&(stardata->star[accretor->starnum]),
stardata,STELLAR_STRUCTURE_CALLER_common_envelope_evolution);
@@ -3827,7 +3830,7 @@ void comenv_accretion(struct stardata_t * stardata,
rc2,
menv,
renv,
- k22,
+ moment_of_inertia_factor2,
&(stardata->star[accretor->starnum]),
stardata,STELLAR_STRUCTURE_CALLER_common_envelope_evolution);
diff --git a/src/evolution/initialize_parameters.c b/src/evolution/initialize_parameters.c
index ed7a1dcaf1fcd3bdb54d49f3de394b42885dd6e0..ae967f8fae3be9367fde5b6a61ba8833b435b63e 100644
--- a/src/evolution/initialize_parameters.c
+++ b/src/evolution/initialize_parameters.c
@@ -367,16 +367,18 @@ void initialize_parameters(struct stardata_t * RESTRICT const stardata)
}
/* and hence set the stellar angular momentum */
- Dprint("J from omega=%g k2=%g rm=%g m=%g mc=%g rc=%g\n",
+ Dprint("J from omega=%g moment_of_inertia_factor=%g rm=%g m=%g mc=%g rc=%g\n",
star->omega,
- star->k2,
+ star->moment_of_inertia_factor,
star->radius,
star->mass,
Outermost_core_mass(star),
star->core_radius);
star->angular_momentum = star->omega *
- moment_of_inertia(star,star->radius);
+ moment_of_inertia(stardata,
+ star,
+ star->radius);
Dprint("J = %g\n",star->angular_momentum);
diff --git a/src/logging/log_every_timestep.c b/src/logging/log_every_timestep.c
index 07bbd3ecf2a6e9a885f01bb534bc2fe54aab53eb..1582348f62de298573788c6bb3c67d08be4f33fd 100644
--- a/src/logging/log_every_timestep.c
+++ b/src/logging/log_every_timestep.c
@@ -2000,7 +2000,10 @@ void log_every_timestep(struct stardata_t * RESTRICT const stardata)
{
star->v_eq = star->vrot0;
star->omega=OMEGA_FROM_VKM*star->v_eq/star->radius;
- star->angular_momentum = star->omega*moment_of_inertia(star,star->radius);
+ star->angular_momentum =
+ star->omega * moment_of_inertia(stardata,
+ star,
+ star->radius);
}
char spectype[10];
@@ -2502,8 +2505,8 @@ void log_every_timestep(struct stardata_t * RESTRICT const stardata)
Second_derivative(model,DERIVATIVE_ORBIT_ANGMOM), //22
Second_derivative(star[0],DERIVATIVE_STELLAR_ANGMOM),//23
Second_derivative(star[1],DERIVATIVE_STELLAR_ANGMOM),//24
- moment_of_inertia(&stardata->star[0],stardata->star[0].radius),//25
- moment_of_inertia(&stardata->star[1],stardata->star[1].radius),//26
+ moment_of_inertia(stardata,&stardata->star[0],stardata->star[0].radius),//25
+ moment_of_inertia(stardata,&stardata->star[1],stardata->star[1].radius),//26
stardata->model.dtm //27
);
}
@@ -2525,10 +2528,10 @@ void log_every_timestep(struct stardata_t * RESTRICT const stardata)
printf("OMEGA ω=%g J=%g I=%g : ω=%g J=%g I=%g\n",
stardata->star[0].omega,
stardata->star[0].angular_momentum,
- moment_of_inertia(&stardata->star[0],stardata->star[0].radius),
+ moment_of_inertia(stardata,&stardata->star[0],stardata->star[0].radius),
stardata->star[1].omega,
stardata->star[1].angular_momentum,
- moment_of_inertia(&stardata->star[0],stardata->star[1].radius));
+ moment_of_inertia(stardata,&stardata->star[0],stardata->star[1].radius));
show_derivatives(stardata);
*/
if(0 &&
@@ -2601,7 +2604,7 @@ void log_every_timestep(struct stardata_t * RESTRICT const stardata)
string[abs(stardata->star[N].dtlimiter)], \
Sync_timescale(N), \
_XHC(N), \
- moment_of_inertia(&stardata->star[N],stardata->star[N].radius)
+ moment_of_inertia(stardata,&stardata->star[N],stardata->star[N].radius)
#define _starheader(N) \
"st%d M%d Mc%d R%d R/RL%d J%d omega%d dMdt%d dJdt%d dt%d dtlimiter%d dtlimstring%d tsync%d XHC%d I%d ", N,N,N,N,N,N,N,N,N,N,N,N,N,N,N \
@@ -3039,7 +3042,8 @@ void log_every_timestep(struct stardata_t * RESTRICT const stardata)
#ifdef __ROTATION_TESTS
if(stardata->star[0].stellar_type <= MAIN_SEQUENCE)
{
- const double I MAYBE_UNUSED = moment_of_inertia(&stardata->star[0],
+ const double I MAYBE_UNUSED = moment_of_inertia(stardata,
+ &stardata->star[0],
stardata->star[0].radius);
const double P = YEAR_LENGTH_IN_DAYS*2.0*PI/stardata->star[0].omega;
printf("ROT %g %g %g %g %g %g\n",
diff --git a/src/logging/show_derivatives.c b/src/logging/show_derivatives.c
index 1d9a5bff2a372b8a3efa66aefd5504cba82332ff..981b7b86a2e3f3b86fa15482ca080f2d33d5380d 100644
--- a/src/logging/show_derivatives.c
+++ b/src/logging/show_derivatives.c
@@ -31,7 +31,7 @@ void show_derivatives(struct stardata_t * RESTRICT const stardata)
(STAR)->roche_radius, \
(STAR)->angular_momentum, \
(STAR)->omega, \
- moment_of_inertia((STAR),((STAR)->radius)), \
+ moment_of_inertia(stardata,(STAR),((STAR)->radius)), \
stardata->store->colours[COLOUR_RESET] \
)
diff --git a/src/rotation/breakup_angular_momentum.c b/src/rotation/breakup_angular_momentum.c
index 1f19aaccbf1d299911a11c8649118a311cc7d7f4..34d0c10d91abb00cd6ed09032349dd82a323da4f 100644
--- a/src/rotation/breakup_angular_momentum.c
+++ b/src/rotation/breakup_angular_momentum.c
@@ -2,7 +2,8 @@
No_empty_translation_unit_warning;
-double breakup_angular_momentum(struct star_t * RESTRICT const star)
+double breakup_angular_momentum(struct stardata_t * const stardata,
+ struct star_t * RESTRICT const star)
{
/*
* Return the maximum angular momentum this star can hold
@@ -13,7 +14,7 @@ double breakup_angular_momentum(struct star_t * RESTRICT const star)
double r = star->radius > star->roche_radius ?
star->effective_radius : star->radius;
calculate_rotation_variables(star,r);
- return moment_of_inertia(star,r) * star->omega_crit;
+ return moment_of_inertia(stardata,star,r) * star->omega_crit;
}
#ifdef __EXPERIMENTAL_CODE
diff --git a/src/rotation/core_moment_of_inertia.c b/src/rotation/core_moment_of_inertia.c
index 0e5fa1054c36dc9b4c384ac23782a35881eefadb..7416ba1d872f7a2879bce73b83233c869110c8d8 100644
--- a/src/rotation/core_moment_of_inertia.c
+++ b/src/rotation/core_moment_of_inertia.c
@@ -8,7 +8,7 @@ No_empty_translation_unit_warning;
* instead
*/
double Pure_function core_moment_of_inertia(const struct star_t * RESTRICT const star,
- const double rc,
+ const double core_radius,
Core_type core_type)
{
if(core_type < 0)
@@ -18,6 +18,6 @@ double Pure_function core_moment_of_inertia(const struct star_t * RESTRICT const
return
HAS_BSE_CORE(star->stellar_type) ?
- (K3_CORE_GYRATION * star->core_mass[core_type] * Pow2(rc)) :
+ (CORE_MOMENT_OF_INERTIA_FACTOR * star->core_mass[core_type] * Pow2(core_radius)) :
0.0;
}
diff --git a/src/rotation/envelope_moment_of_inertia.c b/src/rotation/envelope_moment_of_inertia.c
index 9b8ca23b111a148ebb76cc0cf7f678878ff81773..c8797d5c40fc466cc543d564a0da623aab230745 100644
--- a/src/rotation/envelope_moment_of_inertia.c
+++ b/src/rotation/envelope_moment_of_inertia.c
@@ -4,10 +4,12 @@ No_empty_translation_unit_warning;
/*
* Calculate the moment of inertia of a star's envelope
+ *
+ * This assumes the core contributes little to the moment of inertia.
*/
double envelope_moment_of_inertia(const struct star_t * RESTRICT const star,
const double r)
{
- return star->k2 * (star->mass - Outermost_core_mass(star)) * Pow2(r);
+ return star->moment_of_inertia_factor * (star->mass - Outermost_core_mass(star)) * Pow2(r);
}
diff --git a/src/rotation/moment_of_inertia.c b/src/rotation/moment_of_inertia.c
index 146184ffe7c03dd173f2e5ae2cd97d1ad8db8e10..69548607463de5289a7cb9d49c268b2126e8a112 100644
--- a/src/rotation/moment_of_inertia.c
+++ b/src/rotation/moment_of_inertia.c
@@ -13,17 +13,18 @@ static double normal_star_MOI(const struct star_t * RESTRICT const star,
const double r);
-double Pure_function moment_of_inertia(const struct star_t * RESTRICT const star,
- const double r)
+double moment_of_inertia(struct stardata_t * const stardata,
+ const struct star_t * RESTRICT const star,
+ const double r)
{
- const double momI =
+ double momI =
star->stellar_type == MASSLESS_REMNANT ? TINY :
star->stellar_type == BLACK_HOLE ? BH_MOI(star) :
star->stellar_type == NEUTRON_STAR ? neutron_star_MOI(star,r) :
normal_star_MOI(star,r);
/*
if(1)
- printf("momI [star %d type %d has core? %d, M=%g R=%g (RL=%g) L=%g J=%g] = (Menv=%g-%g) * (k2=%g) * (r=%g)^2 + (k3=%g) * (Mc=%g) * (rc=%g)^2 = %g (%g) + %g (%g) = %g\n",
+ printf("momI [star %d type %d has core? %d, M=%g R=%g (RL=%g) L=%g J=%g] = (Menv=%g-%g) * (moment_of_inertia_factor=%g) * (r=%g)^2 + (k3=%g) * (Mc=%g) * (rc=%g)^2 = %g (%g) + %g (%g) = %g\n",
star->starnum,
star->stellar_type,
HAS_BSE_CORE(star->stellar_type),
@@ -36,25 +37,67 @@ double Pure_function moment_of_inertia(const struct star_t * RESTRICT const star
star->mass,
Outermost_core_mass(star),
- star->k2,
+ star->moment_of_inertia_factor,
r,
K3,
Outermost_core_mass(star),
star->core_radius,
- (star->mass - Outermost_core_mass(star))*star->k2*Pow2(r),
+ (star->mass - Outermost_core_mass(star))*star->moment_of_inertia_factor*Pow2(r),
envelope_moment_of_inertia(star,r),
K3* Pow2(star->core_radius)*Outermost_core_mass(star),
core_moment_of_inertia(star,star->core_radius),
momI);
*/
+ if(stardata->preferences->apply_Darwin_Radau_correction == TRUE)
+ {
+ /*
+ * Darwin-Radau corretion todata into account that
+ * the star is rotating.
+ *
+ * References at
+ * https://en.wikipedia.org/wiki/Darwin%E2%80%93Radau_equation
+ */
+
+ if(Is_not_zero(star->omega_crit))
+ {
+ /* Re/Rp */
+ const double Re_Rp_ratio = r_rot(star->omega/star->omega_crit);
+
+ if(!Fequal(Re_Rp_ratio,1.0))
+ {
+ /* epsilon = 1.0 - Rp/Re */
+ double epsilon = 1.0 - 1.0/Re_Rp_ratio;
+
+ /* q = geodynamical constant */
+ double q = Pow2(star->omega / YEAR_LENGTH_IN_SECONDS) *
+ Pow3(star->radius * R_SUN) /
+ (GRAVITATIONAL_CONSTANT * star->mass * M_SUN);
+
+ /* Radau parameter */
+ const double eta = 5.0/2.0 * q / epsilon - 2.0;
+
+ const double f = (1.0 - 2.0/5.0 * sqrt(Max(0.0,1.0 + eta)));
+
+ printf("Re_Rp_ratio %30.20e \n",Re_Rp_ratio);
+ printf("epsilon %g\n",epsilon);
+ printf("q %g\n",q);
+ printf("5/2 q/epsilon %g\n",5.0/2.0*q/epsilon);
+ printf("eta %g\n",eta);
+
+ printf("multiplier %g\n",
+ f);
+ momI *= f;
+ }
+ }
+ }
+
return momI;
}
static double neutron_star_MOI(const struct star_t * RESTRICT const star,
const double r)
-{
- /*
+{ /*
* Based on Ravenhall and Pethick, 1994
* ApJ 484, 846, Eq. 6.
*
diff --git a/src/rotation/omegaspin.c b/src/rotation/omegaspin.c
index 92ac5593476e0fa302bd745a2619cadeb210f70b..4ce8538f6c87612c7d6f642def435443ae3f6b65 100644
--- a/src/rotation/omegaspin.c
+++ b/src/rotation/omegaspin.c
@@ -12,7 +12,8 @@ No_empty_translation_unit_warning;
* and the star is not a remnant, this is formally
* a bug. In this case we just return 0.0.
*/
-double Pure_function omegaspin(const struct star_t * RESTRICT const star)
+double Pure_function omegaspin(struct stardata_t * stardata,
+ const struct star_t * RESTRICT const star)
{
if(star->stellar_type==MASSLESS_REMNANT ||
Is_zero(star->effective_radius))
@@ -21,7 +22,9 @@ double Pure_function omegaspin(const struct star_t * RESTRICT const star)
}
else
{
- const double I = moment_of_inertia(star,star->effective_radius);
+ const double I = moment_of_inertia(stardata,
+ star,
+ star->effective_radius);
return Is_not_zero(I) ? (star->angular_momentum/I) : 0.0;
}
}
diff --git a/src/rotation/rotation_prototypes.h b/src/rotation/rotation_prototypes.h
index f305e59b249f271e61dd98ef5e27932531ac0c65..30fb14c8bd0834c28bce78f2f19ebc13a323c18f 100644
--- a/src/rotation/rotation_prototypes.h
+++ b/src/rotation/rotation_prototypes.h
@@ -2,7 +2,8 @@
#ifndef ROTATION_PROTOTYPES_H
#define ROTATION_PROTOTYPES_H
-double Pure_function omegaspin(const struct star_t * RESTRICT const star);
+double Pure_function omegaspin(struct stardata_t * const stardata,
+ const struct star_t * RESTRICT const star);
Constant_function double Lang_initial_rotation(const double m);
@@ -13,11 +14,12 @@ Constant_function double r_rot_inv(double f);
void calculate_rotation_variables(struct star_t * RESTRICT const star,
const double radius);
-double Pure_function moment_of_inertia(const struct star_t * RESTRICT const star,
- const double r);
+double moment_of_inertia(struct stardata_t * const stardata,
+ const struct star_t * RESTRICT const star,
+ const double r);
double Pure_function core_moment_of_inertia(const struct star_t * RESTRICT const star,
- const double rc,
+ const double core_radius,
Core_type core_type);
double envelope_moment_of_inertia(const struct star_t * RESTRICT const star,
@@ -33,12 +35,13 @@ double Constant_function bh_angular_velocity(const double M,
Boolean Pure_function stellar_structure_giant_branch_helium_ignition(struct stardata_t * const stardata,
struct star_t * const star);
-double breakup_angular_momentum(struct star_t * RESTRICT const star);
+double breakup_angular_momentum(struct stardata_t * const stardata,
+ struct star_t * RESTRICT const star);
double Pure_function stellar_angular_momentum(const struct star_t * RESTRICT const star,
const int of_what,
- const double r,
- const double rc);
+ const double radius,
+ const double core_radius);
double Rossby_number(struct stardata_t * stardata MAYBE_UNUSED,
struct star_t * const star);
diff --git a/src/rotation/stellar_angular_momentum.c b/src/rotation/stellar_angular_momentum.c
index e0ada4118558e59f4568d8172049b9f9f864c1eb..8e71fab2dff11e13e2418325cb2f5224b65b988d 100644
--- a/src/rotation/stellar_angular_momentum.c
+++ b/src/rotation/stellar_angular_momentum.c
@@ -21,8 +21,8 @@ No_empty_translation_unit_warning;
double Pure_function stellar_angular_momentum(const struct star_t * star,
const int of_what,
- const double r,
- const double rc)
+ const double radius,
+ const double core_radius)
{
double J;
if(of_what == ANGULAR_MOMENTUM_ENVELOPE)
@@ -30,7 +30,7 @@ double Pure_function stellar_angular_momentum(const struct star_t * star,
/*
* Return the angular momentum of the envelope
*/
- J = envelope_moment_of_inertia(star,r) * omega_env;
+ J = envelope_moment_of_inertia(star,radius) * omega_env;
}
else if(of_what == ANGULAR_MOMENTUM_CORE)
{
@@ -39,7 +39,7 @@ double Pure_function stellar_angular_momentum(const struct star_t * star,
*/
if(HAS_BSE_CORE(star->stellar_type))
{
- J = core_moment_of_inertia(star,rc) * omega_core;
+ J = core_moment_of_inertia(star,core_radius) * omega_core;
}
else
{
@@ -51,8 +51,8 @@ double Pure_function stellar_angular_momentum(const struct star_t * star,
/*
* Assume we want the whole star
*/
- J = envelope_moment_of_inertia(star,r) * omega_env +
- core_moment_of_inertia(star,rc) * omega_core;
+ J = envelope_moment_of_inertia(star,radius) * omega_env +
+ core_moment_of_inertia(star,core_radius) * omega_core;
}
return J;
}
diff --git a/src/setup/cmd_line_args_list.h b/src/setup/cmd_line_args_list.h
index 2566c7f8c6f834116b54990aedea31497dd3d103..1a7aa9ae08e9e199da3519777482bf66e59404b2 100644
--- a/src/setup/cmd_line_args_list.h
+++ b/src/setup/cmd_line_args_list.h
@@ -4648,3 +4648,12 @@ ARG_FLOAT_SCANF,
Var(stardata->preferences->colours),
1.0
CMD_LINE_ARG_T_OPTIONAL_ARGS },
+{
+ ARG_SECTION_STARS,
+ "apply_Darwin_Radau_correction",
+ "Apply Darwin-Radau correction to the moment of inertia to take rotation into account?",
+ ARG_BOOLEAN,
+ WTTS_USE_DEFAULT,
+ Var(stardata->preferences->apply_Darwin_Radau_correction),
+ 1.0
+ CMD_LINE_ARG_T_OPTIONAL_ARGS },
diff --git a/src/setup/set_default_preferences.c b/src/setup/set_default_preferences.c
index e69e3e9c90fcc77a8067fbf11af37eb7fb9c3980..20a949d9788248f50d46ec7482505c36db9d61b0 100644
--- a/src/setup/set_default_preferences.c
+++ b/src/setup/set_default_preferences.c
@@ -632,4 +632,5 @@ preferences->WD_accretion_rate_novae_upper_limit_other_donor = ACCRETION_RATE_NO
preferences->Teff_postAGB_min = 6000.0;
preferences->Teff_postAGB_max = 120000.0;
+ preferences->apply_Darwin_Radau_correction = FALSE;
}
diff --git a/src/setup/version.c b/src/setup/version.c
index 5fde506714940e6517e460aed58860759f541180..7a7f6f678bba4c2f5ae7cb3e33a7523b17928c25 100644
--- a/src/setup/version.c
+++ b/src/setup/version.c
@@ -1333,7 +1333,7 @@ void version(struct stardata_t * RESTRICT const stardata)
Show_string_macro(Total_mass);
Show_float_macro(K2_GYRATION);
- Show_float_macro(K3_CORE_GYRATION);
+ Show_float_macro(CORE_MOMENT_OF_INERTIA_FACTOR);
Show_float_macro(MINIMUM_SEPARATION_TO_BE_CALLED_BINARY);
Show_float_macro(MAXIMUM_SEPARATION_TO_BE_CALLED_BINARY);
Show_string_macro(System_is_single);
@@ -1571,7 +1571,7 @@ void version(struct stardata_t * RESTRICT const stardata)
#endif
Macrotest(COMENV_LOG);
Macrotest(GIANT_AGE_FATAL_ERRORS);
- Macrotest(DTR_FUDGES);
+ Macrotest(DTR_CHECKS);
Macrotest(ACQUIRE_FIX_NEGATIVE_AGES);
Macrotest(ACQUIRE_TEST_FOR_NEGATIVE_AGES);
Show_float_macro(ACCRETION_LIMIT_EDDINGTON_MULTIPLIER_STEADY_DEFAULT);
@@ -2086,7 +2086,6 @@ void version(struct stardata_t * RESTRICT const stardata)
#endif
Macrotest(HE_STAR_ACCRETING_H_FIX);
- Macrotest(K2FUDGE);
Macrotest(NOVA_BACK_ACCRETION);
Show_float_macro(POST_CE_ENVELOPE_DM_GB);
diff --git a/src/stellar_structure/convective_envelope_mass_and_radius.c b/src/stellar_structure/convective_envelope_mass_and_radius.c
index 0f162b64d7d750dd5ce597f96b3887c6a25f8ec7..2f9beabf19cb662bfa7573a8fa7d685e6e1a8e3b 100644
--- a/src/stellar_structure/convective_envelope_mass_and_radius.c
+++ b/src/stellar_structure/convective_envelope_mass_and_radius.c
@@ -2,7 +2,6 @@
No_empty_translation_unit_warning;
#ifdef BSE
-#undef K2FUDGE
void convective_envelope_mass_and_radius(struct stardata_t * const stardata,
struct star_t * const star,
@@ -31,7 +30,7 @@ void convective_envelope_mass_and_radius(struct stardata_t * const stardata,
star->v_eq,
star->v_crit_eq,
Yesno(PMS),
- star->k2
+ star->moment_of_inertia_factor
);
/*
@@ -43,17 +42,14 @@ void convective_envelope_mass_and_radius(struct stardata_t * const stardata,
{
star->menv = 0.0;
star->renv = 0.0;
- star->k2 = 0.0;
+ star->moment_of_inertia_factor = 0.0;
return;
}
- double menvg,menvt,menvz,renvg,renvt,renvz,y,u,k2bgb,k2g,k2z;
+ double menvg,menvt,menvz,renvg,renvt,renvz,y,k2bgb,k2g,k2z;
double tebgb=0.0,tauenv,tautms;
const double ajtm = star->age/star->tm;
-#ifdef K2FUDGE
- double k2e_in = star->k2;
-#endif
/*
* A function to estimate the mass and radius of the convective envelope,
* as well as the gyration radius of the envelope.
@@ -83,16 +79,17 @@ void convective_envelope_mass_and_radius(struct stardata_t * const stardata,
* fit is (presumably) not valid for very low-mass stars/planets.
*
* This was the BSE function "mrenv".
+ *
+ * Note: BSE's "k2" is k^2 = I/(MR^2) where I is the moment of inertia,
+ * M,R are the mass and radius. This was renamed
+ * moment_of_inertia_factor in binary_c to avoid confusion with the
+ * gyration constant k_2. (argh!0
*/
const double logm = log10(Max(0.1,star->phase_start_mass));
const double A = Min(0.81,Max(0.68,0.68+0.4*logm));
-
-
- //D = -0.1;
- //E = 0.025;
/*
- * Zero-age and BGB values of k^2.
+ * Zero-age and BGB values of k^2 (== moment_of_inertia_factor = I/(MR^2))
*/
k2z = Min(0.21,Max(0.09-0.27*logm,0.037+0.033*logm));
if(logm>1.3) k2z -= 0.055*Pow2(logm-1.3);
@@ -140,9 +137,7 @@ void convective_envelope_mass_and_radius(struct stardata_t * const stardata,
w = (star->luminosity - star->luminosities[L_BGB])*B;
w *= w;
y = (F - 0.033*log10(star->luminosities[L_BGB]))/k2bgb - 1.0;
-
- y = Max(0.0,y); // should never be < 0! causes fatal errors...
-
+ y = Max(0.0,y);
const double numerator = (F - 0.033*log10(star->luminosity) + 0.4*w);
const double denominator = 1.0+y*(star->luminosities[L_BGB]/star->luminosity)+w;
k2g = numerator / denominator;
@@ -159,9 +154,9 @@ void convective_envelope_mass_and_radius(struct stardata_t * const stardata,
w,
denominator);
Dprint("k2g = numerator/denominator = %g\n",k2g);
- Dprint("convective_envelope_mass_and_radius star->stellar_type=%d k2 mid star->k2=%g k2g=%g from L=%g star->luminosities[L_BGB]=%g star->mass=%g k2bgb=%g w=%g y=%g numerator=%g denominator=%g\n",
+ Dprint("convective_envelope_mass_and_radius star->stellar_type=%d k2 mid star->moment_of_inertia_factor=%g k2g=%g from L=%g star->luminosities[L_BGB]=%g star->mass=%g k2bgb=%g w=%g y=%g numerator=%g denominator=%g\n",
star->stellar_type,
- star->k2,
+ star->moment_of_inertia_factor,
k2g,
star->luminosity,
star->luminosities[L_BGB],star->mass,k2bgb,w,y,
@@ -188,7 +183,7 @@ void convective_envelope_mass_and_radius(struct stardata_t * const stardata,
k2g = k2bgb;
}
- if(effective_stellar_type<GIANT_BRANCH)
+ if(effective_stellar_type < GIANT_BRANCH)
{
menvg = 0.5;
renvg = 0.65;
@@ -242,18 +237,18 @@ void convective_envelope_mass_and_radius(struct stardata_t * const stardata,
Dprint("tau = %g from radius=%g / rzams=%g\n",
tau,star->radius,star->rzams);
const double C = Max(-2.5,Min(-1.5,-2.5+5.0*logm));
- star->k2 = (k2z - 0.025) * pow(tau,C) + 0.025 * pow(tau,-0.1);
+ star->moment_of_inertia_factor = (k2z - 0.025) * pow(tau,C) + 0.025 * pow(tau,-0.1);
Dprint("k2z = %g tau = %g C = %g \n",k2z,tau,C);
- Dprint("stellar type = %d < HeMS : star->k2 set to %g\n",
- effective_stellar_type,star->k2);
+ Dprint("stellar type = %d < HeMS : star->moment_of_inertia_factor set to %g\n",
+ effective_stellar_type,star->moment_of_inertia_factor);
}
else if(effective_stellar_type==HeMS)
{
/*
* Rough fit for naked He MS stars.
*/
- star->k2 = 0.08 - 0.03 * ajtm;
- Dprint("HeMS : star->k2 = %g\n",star->k2);
+ star->moment_of_inertia_factor = 0.08 - 0.03 * ajtm;
+ Dprint("HeMS : star->moment_of_inertia_factor = %g\n",star->moment_of_inertia_factor);
}
else if(effective_stellar_type<HeWD)
{
@@ -261,9 +256,9 @@ void convective_envelope_mass_and_radius(struct stardata_t * const stardata,
* Rough fit for HeHG stars.
*/
const double rzams = rzhef(star->phase_start_mass);
- star->k2 = 0.08 * star->rzams / star->radius;
- Dprint("<HeWD : star->k2 = %g from RZAMS = %g / R=%g (recalc %g)\n",
- star->k2,
+ star->moment_of_inertia_factor = 0.08 * star->rzams / star->radius;
+ Dprint("<HeWD : star->moment_of_inertia_factor = %g from RZAMS = %g / R=%g (recalc %g)\n",
+ star->moment_of_inertia_factor,
star->rzams,
star->radius,
rzams
@@ -274,10 +269,10 @@ void convective_envelope_mass_and_radius(struct stardata_t * const stardata,
* breaks down: don't let this happen! Really
* we need a better solution to fix this ...
*/
- star->k2 = Max(0.005,star->k2);
+ star->moment_of_inertia_factor = Max(0.005,star->moment_of_inertia_factor);
}
- compact_k2 = star->k2;
+ compact_k2 = star->moment_of_inertia_factor;
/*
* tauenv measures proximity to the Hayashi track in terms of Teff.
@@ -298,67 +293,41 @@ void convective_envelope_mass_and_radius(struct stardata_t * const stardata,
tauenv = Max(0.0,Min(1.0,iw*(tauenv-A)));
Dprint("tauenv (proximity to Hayashi track) = %g\n",tauenv);
- // NB iw=1/(1-A)
+ /* NB iw=1/(1-A) */
if(tauenv > 0.0)
{
star->menv = menvg*Pow5(tauenv);
star->renv = renvg*Pow5d4(tauenv);
- double xx=Pow3(tauenv)*(k2g - star->k2);
+ double xx=Pow3(tauenv)*(k2g - star->moment_of_inertia_factor);
Dprint("Near Hayashi track: menv=%g renv=%g\n",
star->menv,
star->renv);
if(effective_stellar_type<HERTZSPRUNG_GAP)
{
- /* Zero-age values for CE mass and radius.
+ /*
+ * Zero-age values for CE mass and radius.
*/
const double x = Max(0.0,Min(1.0,(0.1-logm)/0.55));
- u=Pow5(x);
-
- {
- double a,b;
-
-//#undef ZAMS_MENV_METALLICITY_CORRECTION
-#ifdef ZAMS_MENV_METALLICITY_CORRECTION
- /*
- * fit is valid for low-mass MS stars (m<3)
- * and anyway, m~3 has zero convective envelope
- */
- if(star->phase_start_mass < 3.0 &&
- effective_stellar_type < HERTZSPRUNG_GAP)
- {
- /*
- * Scale the zams menv with Z according to Rob's fit
- */
- z=Max(1e-4,Min(0.02,z)); // limit to range of fit
- a= (-1.24090e-4 + 6.66010)*z;
- b= (3.02190e-1 + -1.50960)*z;
- }
- else
- {
- a=0.18;
- b=0.82;
- }
-#else
- a=0.18,b=0.82; // Jarrod's fit
-#endif
- menvz = a*x + b*u;
- }
-
+ double u = Pow5(x);
+ menvz = 0.18*x + 0.82*u;
renvz = 0.4*Pow1d4(x) + 0.6*Pow2(u);
const double yy = 2.0 + 8.0*x;
- /* Values for CE mass and radius at start of the HG.
+ /*
+ * Values for CE mass and radius at start of the HG.
*/
tautms = Max(0.0,Min(1.0,(tebgb/(sqrt(sqrt(star->luminosities[L_END_MS])/star->TAMS_radius))-A)*iw));
menvt = menvg*Pow5(tautms);
renvt = renvg*Pow5d4(tautms);
- /* Modified expressions during MS evolution.
+
+ /*
+ * Modified expressions during MS evolution.
*/
- u=pow(ajtm,yy);
+ const double u2 = pow(ajtm,yy);
if(tautms>0.0)
{
- star->menv = menvz + u*(star->menv)*(1.0 - menvz/menvt);
- star->renv = renvz + u*(star->renv)*(1.0 - renvz/renvt);
+ star->menv = menvz + u2*(star->menv)*(1.0 - menvz/menvt);
+ star->renv = renvz + u2*(star->renv)*(1.0 - renvz/renvt);
}
else
{
@@ -368,8 +337,8 @@ void convective_envelope_mass_and_radius(struct stardata_t * const stardata,
xx *= u;
}//<HERTZSPRUNG_GAP
- star->k2 += xx;
- Dprint("Correction xxf = %g : star->k2 = %g\n",xx,star->k2);
+ star->moment_of_inertia_factor += xx;
+ Dprint("Correction xxf = %g : star->moment_of_inertia_factor = %g\n",xx,star->moment_of_inertia_factor);
}
else
{
@@ -384,7 +353,7 @@ void convective_envelope_mass_and_radius(struct stardata_t * const stardata,
*/
star->menv = menvg;
star->renv = renvg;
- star->k2 = k2g;
+ star->moment_of_inertia_factor = k2g;
Dprint("Set k2 = k2g = %g\n",k2g);
}
@@ -402,48 +371,26 @@ void convective_envelope_mass_and_radius(struct stardata_t * const stardata,
double f =
log10(Max(1.0,preMS_radius_factor(star->mass,1e6*star->age)));
f = Max(0.0, Min(1.0, f));
- star->k2 = k2g * f + compact_k2 * (1.0 - f);
+ star->moment_of_inertia_factor = k2g * f + compact_k2 * (1.0 - f);
}
}
#endif
- star->menv = Max(star->menv * (star->mass - Outermost_core_mass(star)),1.0E-10);
- star->renv = Max(star->renv * (star->radius - star->core_radius) ,1.0E-10);
-
-#ifdef K2FUDGE
- /*
- * For non-compact stars, do not let the moment of inertia
- * change by a fraction > K2FUDGE_MAX in any timestep.
- *
- * This is a fudge : in reality k2 should be a smooth fit!
- */
- if(effective_stellar_type<HeWD)
- {
- double k2e_change=fabs(k2e_in-star->k2)/(k2e_in);
- double k2e_wanted=star->k2;
- Dprint("k2e change %3.2f %%\n",100.0star->k2_change);
- if(Is_not_zero(k2e_in)&&(k2e_change>K2FUDGE_MAX))
- {
- star->k2 = Min(1.1star->k2_in,Max(0.9star->k2_in,star->k2));
-#ifdef K2FUDGE_WARNING
- printf("Warning : K2 changed by > %g %% : wants to be %g, was %g but capped at %g\n",K2FUDGE_MAX*100.0,k2e_wanted,k2e_in,star->k2);
-#endif
- }
- }
-#endif
+ star->menv = Max(star->menv * (star->mass - Outermost_core_mass(star)),1.0e-10);
+ star->renv = Max(star->renv * (star->radius - star->core_radius) ,1.0e-10);
Dprint("MRENV result menv=%12.12e renv=%12.12e k2e=%g R=%g\n",
star->menv,
star->renv,
- star->k2,
+ star->moment_of_inertia_factor,
star->radius);
- if(star->k2 < 0.0)
+ if(star->moment_of_inertia_factor < 0.0)
{
Backtrace;
Exit_binary_c(BINARY_C_OUT_OF_RANGE,
"k2e = %g <0 error at time %15.12e : INPUT : stellar_type=%d phase_start_mass=%g mass=%g core_mass[CORE_He]=%g (envelope mass %g, envelope radius %g) luminosity=%g radius=%g core_radius=%g age=%g tm=%g luminosities[L_END_MS]=%g luminosities[L_BGB]=%g luminosities[L_HE_IGNITION]=%g rzams=%g TAMS_radius=%g rg=%g menv=%g renv=%g k2e=%g k2g=%g\n",
- star->k2,
+ star->moment_of_inertia_factor,
stardata->model.time,
star->stellar_type,
star->phase_start_mass,
@@ -464,7 +411,7 @@ void convective_envelope_mass_and_radius(struct stardata_t * const stardata,
rg,
star->menv,
star->renv,
- star->k2,
+ star->moment_of_inertia_factor,
k2g);
}
diff --git a/src/stellar_structure/set_star_struct_from_stellar_structure.c b/src/stellar_structure/set_star_struct_from_stellar_structure.c
index 3e352d2fd0cb2dd6224efb4bfae218984d8a1500..a75ced49acb2ffc9b2d5b0a1af283b116b5fc9a0 100644
--- a/src/stellar_structure/set_star_struct_from_stellar_structure.c
+++ b/src/stellar_structure/set_star_struct_from_stellar_structure.c
@@ -23,7 +23,7 @@ void set_star_struct_from_stellar_structure(struct star_t * const newstar,
oldstar->tn = newstar->tn;
oldstar->renv = newstar->renv;
oldstar->stellar_type = newstar->stellar_type;
- oldstar->k2 = newstar->k2;
+ oldstar->moment_of_inertia_factor = newstar->moment_of_inertia_factor;
oldstar->SN_type = newstar->SN_type;
oldstar->phase_start_mass = newstar->phase_start_mass;
oldstar->phase_start_core_mass = newstar->phase_start_core_mass;
diff --git a/src/stellar_structure/set_stellar_structure_struct_from_star.c b/src/stellar_structure/set_stellar_structure_struct_from_star.c
index 4a852b8ba2d6e1932f9e05e1712e798790fea13d..165a447469bcb7871099ab7c9ba50d94482f7be5 100644
--- a/src/stellar_structure/set_stellar_structure_struct_from_star.c
+++ b/src/stellar_structure/set_stellar_structure_struct_from_star.c
@@ -19,7 +19,7 @@ void set_stellar_structure_struct_from_star(struct star_t * const newstar,
set(tm);
set(tn);
set(renv);
- set(k2);
+ set(moment_of_inertia_factor);
set(phase_start_mass);
set(phase_start_core_mass);
set(core_mass[CORE_CO]);
diff --git a/src/stellar_structure/stellar_structure_BSE_given_timescales.c b/src/stellar_structure/stellar_structure_BSE_given_timescales.c
index e62b0a1d9d58d8cd386b0561235f8506da3f0d5c..d671b1c868d46d5f047801552bd6541de8c085f0 100644
--- a/src/stellar_structure/stellar_structure_BSE_given_timescales.c
+++ b/src/stellar_structure/stellar_structure_BSE_given_timescales.c
@@ -71,7 +71,7 @@ int stellar_structure_BSE_given_timescales(struct star_t *newstar,
newstar->starnum,
RLOF_CALL);
#ifdef DEBUG
- newstar->k2 = 0.0;
+ newstar->moment_of_inertia_factor = 0.0;
#endif // DEBUG
Dprint("phase_start_mass=%g age=%g mt=%g mc=%g st=%d mc=%g k2=%g star->phase_start_mass=%g ntp=%g angular_momentum=%g\n",
newstar->phase_start_mass,
@@ -80,7 +80,7 @@ int stellar_structure_BSE_given_timescales(struct star_t *newstar,
Outermost_core_mass(newstar),
newstar->stellar_type,
Outermost_core_mass(newstar),
- newstar->k2,
+ newstar->moment_of_inertia_factor,
oldstar->phase_start_mass,
newstar->num_thermal_pulses,
newstar->angular_momentum);
@@ -462,14 +462,14 @@ int stellar_structure_BSE_given_timescales(struct star_t *newstar,
newstar->luminosities[L_BGB],
stardata->common.giant_branch_parameters);
- Dprint("convective_envelope_mass_and_radius pre k2=%g\n",newstar->k2);
+ Dprint("convective_envelope_mass_and_radius pre k2=%g\n",newstar->moment_of_inertia_factor);
convective_envelope_mass_and_radius(stardata,
newstar,
rg,
stardata->common.metallicity);
- Dprint("convective_envelope_mass_and_radius post k2=%g\n",newstar->k2);
+ Dprint("convective_envelope_mass_and_radius post k2=%g\n",newstar->moment_of_inertia_factor);
}
else
{
@@ -482,7 +482,7 @@ int stellar_structure_BSE_given_timescales(struct star_t *newstar,
/* set core radius to actual radius */
newstar->core_radius = newstar->radius;
/* generic k2 */
- newstar->k2 = 0.21;
+ newstar->moment_of_inertia_factor = 0.21;
}
#ifdef NUCSYN
@@ -496,11 +496,11 @@ int stellar_structure_BSE_given_timescales(struct star_t *newstar,
stardata);
#endif//NUCSYN
- if(newstar->k2<0.0)
+ if(newstar->moment_of_inertia_factor<0.0)
{
Exit_binary_c(BINARY_C_OUT_OF_RANGE,
"Negative k2 = %g at stellar_structure exit! (mass=%g mt=%g core_mass[CORE_He]=%g lum=%g) \n",
- newstar->k2,
+ newstar->moment_of_inertia_factor,
newstar->phase_start_mass,
newstar->mass,
Outermost_core_mass(newstar),
@@ -533,7 +533,7 @@ int stellar_structure_BSE_given_timescales(struct star_t *newstar,
newstar->core_radius,
newstar->menv,
newstar->renv,
- newstar->k2,
+ newstar->moment_of_inertia_factor,
newstar->omega,
newstar->angular_momentum,
newstar->num_thermal_pulses);
@@ -579,7 +579,7 @@ int stellar_structure_BSE_given_timescales(struct star_t *newstar,
isnan(newstar->core_mass[CORE_Si]) ||
isnan(newstar->core_mass[CORE_Fe]) ||
isnan(newstar->core_mass[CORE_NEUTRON]) ||
- isnan(newstar->k2) ||
+ isnan(newstar->moment_of_inertia_factor) ||
isnan(newstar->luminosity) ||
isnan(newstar->renv) ||
isnan(newstar->angular_momentum) ||
@@ -599,7 +599,7 @@ int stellar_structure_BSE_given_timescales(struct star_t *newstar,
newstar->core_mass[CORE_NEUTRON],
newstar->luminosity,
newstar->radius,
- newstar->k2,
+ newstar->moment_of_inertia_factor,
newstar->renv,
newstar->angular_momentum,
newstar->omega);
diff --git a/src/stellar_structure/stellar_structure_HeStar.c b/src/stellar_structure/stellar_structure_HeStar.c
index 37fefa18730e2fa778c6f4c3a5616d55b27498a0..7eb05382fd38dcc014b8cc1edebffcd08ec423c2 100644
--- a/src/stellar_structure/stellar_structure_HeStar.c
+++ b/src/stellar_structure/stellar_structure_HeStar.c
@@ -483,13 +483,13 @@ tau = newstar->age / newstar->tm;
*/
double J=
// envelope
- newstar->k2*M_SUN*(newstar->mass-newstar->core_mass[CORE_CO])*
+ newstar->moment_of_inertia_factor*M_SUN*(newstar->mass-newstar->core_mass[CORE_CO])*
Pow2(newstar->radius*R_SUN)*(
newstar->omega
/YEAR_LENGTH_IN_SECONDS)
+
// core
- K3_CORE_GYRATION*(newstar->core_mass[CORE_CO]*M_SUN)*Pow2(newstar->core_radius*R_SUN)*(
+ CORE_MOMENT_OF_INERTIA_FACTOR*(newstar->core_mass[CORE_CO]*M_SUN)*Pow2(newstar->core_radius*R_SUN)*(
newstar->omega
/YEAR_LENGTH_IN_SECONDS);
diff --git a/src/stellar_structure/stellar_structure_remnant_and_perturbations.c b/src/stellar_structure/stellar_structure_remnant_and_perturbations.c
index 778a965bef4d2473488ee30380b01f6349d81755..05fe38554162351452ebaf40ba8e4cdb0d36b4a0 100644
--- a/src/stellar_structure/stellar_structure_remnant_and_perturbations.c
+++ b/src/stellar_structure/stellar_structure_remnant_and_perturbations.c
@@ -271,15 +271,12 @@ void stellar_structure_remnant_and_perturbations(struct star_t * newstar,
newstar->core_radius = newstar->radius;
newstar->menv = 1E-10;
newstar->renv = 1E-10;
- newstar->k2 = 0.21;
+ newstar->moment_of_inertia_factor = 0.21;
}
/*
* Perturb the luminosity and radius due to small envelope mass.
- *
- * Note: only do this if we're not about to explode as a SN.
*/
- //if(newstar->SN_type == SN_NONE)
{
stellar_structure_small_envelope_perturbations(
stardata,
diff --git a/src/supernovae/supernova.c b/src/supernovae/supernova.c
index 47fd93bb39ac45fc6bd0d98476e695391abe4228..73e1d9ddb5ae5cfede56aae04fe67b9ad640fe57 100644
--- a/src/supernovae/supernova.c
+++ b/src/supernovae/supernova.c
@@ -145,7 +145,9 @@ void supernova(struct stardata_t * const stardata,
else
{
star->omega = star->angular_momentum /
- moment_of_inertia(star,star->radius);
+ moment_of_inertia(stardata,
+ star,
+ star->radius);
}
diff --git a/src/tides/tides_generic.c b/src/tides/tides_generic.c
index 6905ae2ac1043d59b300cd126fc95009505e6965..a57449b9783d0913d445d4b6e38855e30d03fc8f 100644
--- a/src/tides/tides_generic.c
+++ b/src/tides/tides_generic.c
@@ -41,7 +41,9 @@ void tides_generic(struct stardata_t * const stardata,
{
double tc,f,tcqr,rg2=0.0,ttid;
const double cq = object_mass / star->mass;
- const double starI = moment_of_inertia(star,star_radius);
+ const double starI = moment_of_inertia(stardata,
+ star,
+ star_radius);
const double ecc2 = Pow2(orbit->eccentricity);
const double omecc2 = 1.0 - ecc2; // 1-e^2
const double sqome2 = sqrt(omecc2); // sqrt(1-e^2)
@@ -120,7 +122,7 @@ void tides_generic(struct stardata_t * const stardata,
cq,
raa6,
tcqr);
- rg2 = star->k2;
+ rg2 = star->moment_of_inertia_factor;
Dprint("star%d: Zahn radiative damping tcqr=%g rg2=%g\n",
star->starnum,
tcqr,
@@ -168,14 +170,14 @@ void tides_generic(struct stardata_t * const stardata,
itc,
cq,
raa6);
- rg2 = star->k2 * Max(1e-30, menv) * im;
- Dprint("star%d: Convective (Hut) damping tcqr=%g cq=%g rg2=%g raa6=%g k2=%g menv=%g Envelope_mass=%g im=%g itc=%g\n",
+ rg2 = star->moment_of_inertia_factor * Max(1e-30, menv) * im;
+ Dprint("star%d: Convective (Hut) damping tcqr=%g cq=%g rg2=%g raa6=%g moment_of_inertia_factor=%g menv=%g Envelope_mass=%g im=%g itc=%g\n",
star->starnum,
tcqr,
cq,
raa6,
rg2,
- star->k2,
+ star->moment_of_inertia_factor,
star->menv,
menv,
im,
@@ -197,14 +199,14 @@ void tides_generic(struct stardata_t * const stardata,
star->menv,
im,
itc);
- rg2 = star->k2 * Max(1e-30,menv) * im;
- Dprint("star%d: Convective (Hut) damping tcqr=%g cq=%g rg2=%g raa6=%g k2=%g menv=%g Envelope_mass=%g im=%g itc=%g\n",
+ rg2 = star->moment_of_inertia_factor * Max(1e-30,menv) * im;
+ Dprint("star%d: Convective (Hut) damping tcqr=%g cq=%g rg2=%g raa6=%g moment_of_inertia_factor=%g menv=%g Envelope_mass=%g im=%g itc=%g\n",
star->starnum,
tcqr,
cq,
raa6,
rg2,
- star->k2,
+ star->moment_of_inertia_factor,
star->menv,
menv,
im,
@@ -226,7 +228,7 @@ void tides_generic(struct stardata_t * const stardata,
*/
f = 7.33e-09*pow(star->luminosity/star->mass,5.0/7.0);
tcqr = f * Pow2(raa2*cq) / (1.0+cq);
- rg2 = K3_CORE_GYRATION;
+ rg2 = CORE_MOMENT_OF_INERTIA_FACTOR;
Dprint("star%d: Degenerate (Campbell) damping tcqr=%g rg2=%g\n",
star->starnum,
tcqr,
diff --git a/src/timestep/timestep_TPAGB.c b/src/timestep/timestep_TPAGB.c
index 2c29d3b287f186f5c0fa1b0fb28eee805e923a8d..cc8c7dda09b37455aa83d285c91aeb068eda6ba4 100644
--- a/src/timestep/timestep_TPAGB.c
+++ b/src/timestep/timestep_TPAGB.c
@@ -73,14 +73,14 @@ void timestep_TPAGB(Timestep_prototype_args)
*time_remaining = tn - age;
-#ifdef DTR_FUDGES
+#ifdef DTR_CHECKS
// Rob's fudge for small envelopes, sometimes things get
// unstable
if(*time_remaining<1e-9)
{
*time_remaining = Max(*dt,*time_remaining);
}
-#endif //DTR_FUDGES
+#endif //DTR_CHECKS
}
else
{
diff --git a/src/timestep/timestep_hard_limits.c b/src/timestep/timestep_hard_limits.c
index 234a398e129ad2c92e8aff3db618c7803cbc9142..3699de2c3df20fe9aad1c5b32bc7a43ffd76a50c 100644
--- a/src/timestep/timestep_hard_limits.c
+++ b/src/timestep/timestep_hard_limits.c
@@ -5,7 +5,7 @@ No_empty_translation_unit_warning;
void timestep_hard_limits(Timestep_prototype_args)
{
-#ifdef DTR_FUDGES
+#ifdef DTR_CHECKS
if(*time_remaining < 0.0)
{
/*
diff --git a/src/zfuncs/zfuncs_prototypes.h b/src/zfuncs/zfuncs_prototypes.h
index cc77ddbc9908b588c82208874e28bba0d04abfd8..a944d6ab0d94e91c7467382dafbc8762cf1e2d93 100644
--- a/src/zfuncs/zfuncs_prototypes.h
+++ b/src/zfuncs/zfuncs_prototypes.h
@@ -334,10 +334,4 @@ double time_from_preMS_radius_factor(struct stardata_t * RESTRICT const stardata
const double radius_factor);
#endif // PRE_MAIN_SEQUENCE
-double Pure_function moment_of_inertia(const struct star_t * RESTRICT const star,
- const double r);
-
-
-
-
#endif /* ZFUNCS_PROTOTYPES_H */
diff --git a/tbse b/tbse
index 09b6d2017b512b9d32f6f577c2ff96b6e92043de..c5820129a90da321e59c430d2517a01546dd636d 100755
--- a/tbse
+++ b/tbse
@@ -240,17 +240,22 @@ PRE_MAIN_SEQUENCE_FIT_LOBES=False
# White dwarf cooling model
-# WHITE_DWARF_COOLING_MESTEL = 0 : Mestel
+# WHITE_DWARF_COOLING_MESTEL = 0 : Mestel (as in BSE/Hurley+2002)
# WHITE_DWARF_COOLING_MESTEL_MODIFIED = 1 : Mestel but modified
# WHITE_DWARF_COOLING_CARRASCO2014 = 2 : Carrasco 2014 tables
WHITE_DWARF_COOLING_MODEL=WHITE_DWARF_COOLING_MESTEL
# White dwarf radius model
-# WHITE_DWARF_RADIUS_NAUENBERG1972=0 Nauenberg (1972)
+# WHITE_DWARF_RADIUS_NAUENBERG1972=0 Nauenberg (1972, as in BSE/Hurley+2002)
# WHITE_DWARF_RADIUS_MU=1 mu-dependent variant of Nauenberg (1972)
# WHITE_DWARF_RADIUS_CARRASCO2014=2 Carrasco et al. (2014) tables
WHITE_DWARF_RADIUS_MODEL=WHITE_DWARF_RADIUS_NAUENBERG1972
+# Apply Darwin-Radau correction to the moment of inertia
+# to take rotation into account? Doesn't work well at
+# low rotation rates - considered experimental! (FALSE)
+APPLY_DARWIN_RADAU_CORRECTION=FALSE
+
############################################################
#
# Binary star physics
@@ -2630,6 +2635,7 @@ $ORBITING_OBJECTS \
--wind_Nieuwenhuijzen_luminosity_lower_limit $WIND_NIEUWENHUIJZEN_LUMINOSITY_LOWER_LIMT \
--wind_LBV_luminosity_lower_limit $WIND_LBV_LUMINOSITY_LOWER_LIMT \
--float_overflow_checks $FLOAT_OVERFLOW_CHECKS \
+--apply_Darwin_Radau_correction $APPLY_DARWIN_RADAU_CORRECTION \
";
# valgrind options
diff --git a/ubuntu_install_binary_c.sh b/ubuntu_install_binary_c.sh
index c0f393ebf90bf56ce371c6d5b2f9cbcf54e83560..4679da756a9c48f01cc1471ff3b19befeedee68b 100755
--- a/ubuntu_install_binary_c.sh
+++ b/ubuntu_install_binary_c.sh
@@ -116,7 +116,7 @@ git_rev ()
echo ${c}:${d}:${h}
}
-alias new="ls -lrt |head -10"
+alias new="ls -lrt | tail -10"
# git function to find the repository url
git_url ()
@@ -288,16 +288,21 @@ function install_perl
source ~/perl5/perlbrew/etc/bashrc 2>&1 || true
echo "Install cpanm"
- perlbrew -f install-cpanm
+ perlbrew -f install-cpanmu
}
function install_binary_c_python
{
# install the binary_c Python module
export LIBRARY_PATH=$LD_LIBRARY_PATH
- pip3 install -U pip setuptools
+
+ # install pip and setup tools
+ pip3 install -U pip setuptools testresources
+
+ # choose a binary_c-python install method
# pip3 install git+https://$GITLAB/binary_c-python.git
pip3 install git+https://gitlab.eps.surrey.ac.uk/ri0005/binary_c-python.git@development_version/2.1.7
+
}
function update_packages
@@ -317,8 +322,21 @@ function update_packages
echo "You should normally choose Python 3.9"
sudo update-alternatives --config python # >>> choose python3.9 <<<
sudo apt-get clean
- pip install --upgrade meson
- pip install --upgrade ninja
+
+ # update python packages required for the build
+ pip3 install --upgrade meson
+ pip3 install --upgrade ninja
+}
+
+function update_python_packages
+{
+ # update python packages and required system packages
+ sudo apt-get install libgpg-error-dev libgpgme-dev swig libcairo2-dev libgirepository1.0-dev autoconf python-dev
+
+ pip3 install pipupgrade
+ pip3 install --upgrade beautifulsoup4 html5lib
+ pipupgrade --self
+ pipupgrade --verbose --yes
}
function download_MINT_data