diff --git a/src/binary_c_macros.h b/src/binary_c_macros.h
index 2eb41d9d4a8e7e8fb0fdad92394501a66ef3004f..194feb8a9824e03c25f8ab2b3191fb9d7db2abe9 100644
--- a/src/binary_c_macros.h
+++ b/src/binary_c_macros.h
@@ -42,6 +42,7 @@
#include "supernovae/sn.h"
#include "file/file_macros.h"
#include "disc/disc_macros.h"
+#include "disc/disc_function_macros.h"
#include "nucsyn/nucsyn_macros.h"
#include "evolution/rejection_macros.h"
@@ -907,7 +908,7 @@
#define L_SOLAR "L☉"
#define M_SOLAR "M☉"
#define R_SOLAR "R☉"
-
+#define J_CGS "g cm^2 s^-1"
/*
* Magnetic braking algorithms
diff --git a/src/binary_c_maths.h b/src/binary_c_maths.h
index af7ab3bfb637bb86e112529ea248ca46ffff332e..aa148f12c1c5b0be6da59980eedaea7b7708b827 100644
--- a/src/binary_c_maths.h
+++ b/src/binary_c_maths.h
@@ -524,6 +524,12 @@ Power_function func_pow4d3(const double x) {return _Pow4d3_macro(x);}
*/
#define Safepow10(X) (isnan(X) ? (X) : Pow10(X)))
+/*
+ * Macro to check if X is a non-zero number
+ */
+#define Is_nonzero_number(X) \
+ ((Is_really_not_zero(X) && \
+ !isnan(X) && !isinf(X)))
#endif // BINARY_MATHS_H
diff --git a/src/binary_c_structures.h b/src/binary_c_structures.h
index 28e553a382964290646a5792ceb6d9d22ea6f035..89e8cc0d64403c85777a627b504d7842c71c1001 100644
--- a/src/binary_c_structures.h
+++ b/src/binary_c_structures.h
@@ -53,6 +53,7 @@
#ifdef DISCS
#include "disc/disc_parameters.h"
#include "disc/disc_macros.h"
+#include "disc/disc_function_macros.h"
#include "disc/disc_power_laws.h"
#include "disc/disc_constraints.h"
#endif//DISCS
diff --git a/src/disc/disc.h b/src/disc/disc.h
index 4f6ed84c2105e37c134ebf1f4be8ec80046fdcc2..940c1acbf52f93c3c1d1459caf4ff9c6672c48b9 100644
--- a/src/disc/disc.h
+++ b/src/disc/disc.h
@@ -11,6 +11,7 @@
#include "disc_parameters.h"
#include "disc_macros.h"
+#include "disc_function_macros.h"
#include "disc_power_laws.h"
#include "disc_prototypes.h"
#include "disc_thermal_zones.h"
diff --git a/src/disc/disc_bisection_rooter.c b/src/disc/disc_bisection_rooter.c
index 7717c8da293ad3beba36cafef6992e0d7e02364d..bb3c1f59ab08d52ce3f54f04966b9a0ef8a17826 100644
--- a/src/disc/disc_bisection_rooter.c
+++ b/src/disc/disc_bisection_rooter.c
@@ -14,8 +14,6 @@ static void disc_nsector_wrapper2(const int n,
/*
* Hard-coded bisection rooter
*/
-#undef DISC_DEBUG
-#define DISC_DEBUG 2
int disc_bisection_rooter(const int n,
const disc_parameter * parameters,
@@ -32,7 +30,7 @@ int disc_bisection_rooter(const int n,
int error;
unsigned int i,j;
- printf("BISECTOR at t = %g y\n",disc->lifetime / YEAR_LENGTH_IN_SECONDS);
+ //printf("BISECTOR at t = %g y\n",disc->lifetime / YEAR_LENGTH_IN_SECONDS);
while(converged == FALSE && failed == FALSE)
{
@@ -117,23 +115,28 @@ int disc_bisection_rooter(const int n,
"Bisection gave error %d \"%s\"\n",
error,
Bisect_error_string(error));
-
- if(disc->lifetime > 1000.0 * YEAR_LENGTH_IN_SECONDS)
+
+#ifdef ____DEBUG
+ if(0 && disc->lifetime > 1000.0 * YEAR_LENGTH_IN_SECONDS)
{
+ /*
printf("Want disc->M = %g have %g\n",
disc->M,
disc_total_mass(disc));
-
+ */
double Tvwas = disc->Tvisc0;
- double dTv = 0.001 * disc->Tvisc0;
+ double dTv = 1e-4 * disc->Tvisc0;
struct disc_t * discwas = New_disc_from(disc);
- double Tv;
+ double Tv = 1.0819e17;
+ double Tvmax = 1.09e17;
double residual[n];
int failure;
- for(Tv = 0.9 * Tvwas;
+
+ while(Tv < Tvmax)
+/* for(Tv = 0.9 * Tvwas;
Tv <= 1.1 * Tvwas;
Tv += dTv)
- {
+*/ {
Copy_disc(discwas,disc);
disc->Tvisc0 = Tv;
disc_nsector2(stardata,
@@ -143,20 +146,28 @@ int disc_bisection_rooter(const int n,
constraints,
residual,
&failure);
+ const double eps = (disc_total_mass(disc) - disc->M)/disc->M;
printf("CHECKNEAR %30.15e %30.15g eps %g nzones %u zone 0 type %d Rin %g Rout %g\n",
disc->Tvisc0,
disc_total_mass(disc),
- (disc_total_mass(disc) - disc->M)/disc->M,
+ eps,
disc->n_thermal_zones,
disc->thermal_zones[0].type,
disc->Rin,
disc->Rout
);
-
+ disc_show_thermal_zones(disc,NULL,-1);
+ if(eps > 0.0)
+ {
+ Flexit;
+ }
+ Tv += dTv;
}
- disc_monotonic_check(stardata,disc,binary,n);
+
+ //disc_monotonic_check(stardata,disc,binary,n);
Flexit;
}
+#endif
if(error != BINARY_C_BISECT_ERROR_NONE)
{
failed = i + 1;
@@ -333,8 +344,6 @@ static double disc_nsector_wrapper(const double x,
const double parameter_value[1] = { x };
double residual[1];
-
- printf("parameter value %g\n",x);
disc_nsector_wrapper2(n,
parameter_value,
@@ -348,14 +357,6 @@ static double disc_nsector_wrapper(const double x,
failed
);
- printf("disc mass %g = %g Msun : want %g = %g Msun -> residual %g\n",
- disc_total_mass(disc),
- disc_total_mass(disc)/M_SUN,
- disc->M,
- disc->M/M_SUN,
- residual[0]
- );
-
((struct GSL_args_t *)p)->error = *failed;
Discdebug(2,
diff --git a/src/disc/disc_calc_disc_structure.c b/src/disc/disc_calc_disc_structure.c
index c54a20eddc1f49a41ce7caaeef55cabbb42c693a..7e54a5e2ceaa956403b476b3db1aaa6404e3987a 100644
--- a/src/disc/disc_calc_disc_structure.c
+++ b/src/disc/disc_calc_disc_structure.c
@@ -12,9 +12,6 @@
#include "disc_constraints.h"
#include "disc_boundary_conditions.h"
-#undef DISC_DEBUG
-#define DISC_DEBUG 2
-
/************************************************************
* disc_calc_disc_structure
*
@@ -144,8 +141,7 @@ int disc_calc_disc_structure(const struct binary_system_t * const binary,
/*
* Converge, if possible, with the default inner edge torque.
*/
-printf("call %d\n",n);
-converge_with_torque(n,
+ converge_with_torque(n,
parameter_types,
constraints,
binary,
@@ -998,11 +994,6 @@ static int disc_GSL_multiroot(const int solver,
"Failed to allocated memory for GSL multiroot fsolver\n");
}
#endif//ALLOC_CHECKS
-
- /*
- * use logs?
- */
-#define USELOG
/*
* Set up the test function
@@ -1033,11 +1024,11 @@ static int disc_GSL_multiroot(const int solver,
gsl_vector_set(
guess_vector,
(size_t)i,
-#ifdef USELOG
+#ifdef DISC_NSECTION_USELOG
log(x)
#else
x
-#endif//USELOG
+#endif//DISC_NSECTION_USELOG
);
}
@@ -1149,17 +1140,6 @@ static int disc_GSL_multiroot(const int solver,
status,
gsl_strerror(status),
iter);
-
- /* debugging */
- if(status==27 && 0)
- {
- printf("no progress error\n");
- disc_monotonic_check(stardata,
- disc,
- binary,
- n);
- Flexit;
- }
break;
}
@@ -1244,12 +1224,9 @@ static int GSL_multiroot_testfunc(const gsl_vector * x,
{
parameter_values[i] =
gsl_vector_get(x,i);
-#ifdef USELOG
+#ifdef DISC_NSECTION_USELOG
parameter_values[i] = exp(parameter_values[i]);
#endif
- printf("parameter %d -> %g\n",
- i,
- parameter_values[i]);
}
double residuals[n];
diff --git a/src/disc/disc_function_macros.h b/src/disc/disc_function_macros.h
new file mode 100644
index 0000000000000000000000000000000000000000..a8dad9c73fa7d0d06935a7f2f92b07854660fbc5
--- /dev/null
+++ b/src/disc/disc_function_macros.h
@@ -0,0 +1,154 @@
+#pragma once
+#ifndef DISC_FUNCTION_MACROS_H
+#define DISC_FUNCTION_MACROS_H
+
+
+/* convert position to disc index for an existing disc */
+#define Disc_pos_to_n(OBJ,POS) \
+ ((POS) == DISC_INNER_EDGE ? 0 : \
+ (POS) == DISC_OUTER_EDGE ? (OBJ)->ndiscs-1 : \
+ (POS))
+
+/* convert position to disc index for a new disc */
+#define Disc_newpos_to_n(OBJ,POS) \
+ ((POS) == DISC_INNER_EDGE ? 0 : \
+ (POS) == DISC_OUTER_EDGE ? (OBJ)->ndiscs : \
+ (POS))
+
+#define Discloop(OBJ,A) \
+ for((A)=0;(A)<(OBJ)->ndiscs;(A)++)
+
+
+
+#define Convergence_status(S) \
+ ( \
+ (S) == CONVERGENCE_FAILED ? "Failed" : \
+ (S) == CONVERGENCE_SUCCEEDED ? "Success" : \
+ (S) == CONVERGENCE_FLUX_NEGATIVE ? "Flux negative" : \
+ "Unknown" \
+ )
+
+
+
+/*
+ * Macro to determine whether a disc "really exists".
+ *
+ * This can often be used to determine whether other
+ * functions, like integrals over the disc, will work.
+ */
+#define Disc_is_disc(D) \
+ ( \
+ Is_nonzero_number((D)->M) \
+ && \
+ Is_nonzero_number((D)->J) \
+ && \
+ (D)->lifetime < DISC_MAX_LIFETIME \
+ )
+
+
+
+#define Disc_removal_string(S) \
+ ( \
+ (S)==DISC_REMOVE_FORMED_RING ? "Formed ring" : \
+ (S)==DISC_REMOVE_REPLACE ? "Replaced" : \
+ (S)==DISC_REMOVE_EVAPORATE ? "Evaporated" : \
+ (S)==DISC_REMOVE_ALL ? "Remove all" : \
+ "Unknown" \
+ )
+
+
+/* solvers */
+#define Solver_string(N) ( \
+ (N)==(DISC_SOLVER_BISECTION) ? "Bi" : \
+ (N)==(GSL_MULTIROOT_FSOLVER_HYBRIDS) ? "Hs" : \
+ (N)==(GSL_MULTIROOT_FSOLVER_HYBRID) ? "H" : \
+ "?" \
+ )
+
+
+#define Show_disc_derivative(N) \
+ Discdebug(1, \
+ "d(%20s)/dt = M % 20g (Msun/y), J % 20g (cgs), e % 20g (s^-1)\n", \
+ Disc_derivative_string(N), \
+ disc->loss[(N)].mdot / M_SUN * YEAR_LENGTH_IN_SECONDS, \
+ disc->loss[(N)].jdot, \
+ disc->loss[(N)].edot);
+
+
+#define Clear_disc_feedback(D) \
+ { \
+ (D)->dM_ejected = \
+ (D)->dM_binary = \
+ (D)->dJ_binary = \
+ (D)->de_binary = 0.0; \
+ }
+
+#define Disc_derivative_string(A) ( \
+ (A) == DISC_LOSS_GLOBAL ? "Global" : \
+ (A) == DISC_LOSS_INNER_VISCOUS ? "Viscous" : \
+ (A) == DISC_LOSS_INNER_L2_CROSSING ? "L2 crossing": \
+ (A) == DISC_LOSS_FUV ? "FUV" : \
+ (A) == DISC_LOSS_XRAY ? "X-ray" : \
+ (A) == DISC_LOSS_BINARY_TORQUE ? "Binary torque" : \
+ (A) == DISC_LOSS_ISM ? "ISM ram stripping" : \
+ (A) == DISC_LOSS_RESONANCES ? "Resonances" : \
+ (A) == DISC_LOSS_NEGATIVE_FLUX ? "Negative flux" : \
+ (A) == DISC_LOSS_INNER_EDGE ? "Inner edge" : \
+ (A) == DISC_LOSS_OUTER_EDGE ? "Outer edge" : \
+ "Unknown" )
+
+
+/* macro to define mass loss from the inner or outer edge */
+#define Disc_loss_from_edge(I) \
+ ((I)==DISC_LOSS_INNER_EDGE || \
+ (I)==DISC_LOSS_OUTER_EDGE)
+
+
+/*
+ * ... and a string to log them.
+ */
+#define Disc_logger(A) ( \
+ (A) == DISC_LOG_EVOLVE_DISC_STRUCTURE ? "disc_evolve_disc_structure" : \
+ (A) == DISC_LOG_EVOLVE_DISC_STRUCTURE2 ? "disc_evolve_disc_structure (2)" : \
+ (A) == DISC_LOG_EVERY_TIMESTEP ? "log_every_timestep" : \
+ "Unknown")
+
+
+
+/*
+ * Copy disc from source S to destination D
+ */
+#define Copy_disc(S,D) memcpy((D),(S),sizeof(struct disc_t));
+
+/*
+ * New disc copied from another
+ */
+#define New_disc_from(D) Copy_disc((D),New_disc)
+
+
+
+#define Evaporate_disc(D,REASON) \
+ (D)->M = 0.0; \
+ (D)->J = 0.0; \
+ (D)->dt = LONG_TIMESTEP; \
+ (D)->converged = FALSE; \
+ if(DISC_DEBUG) \
+ printf("Disc evaporated at line %d in file %s because %s\n", \
+ __LINE__, \
+ __FILE__, \
+ REASON); \
+ /* disc evaporated */
+
+#define Disc_timestep_limit_string(X) \
+ ( \
+ (X) == DISC_TIMESTEP_LIMIT_M ? "M" : \
+ (X) == DISC_TIMESTEP_LIMIT_J ? "J" : \
+ (X) == DISC_TIMESTEP_LIMIT_F ? "F" : \
+ (X) == DISC_TIMESTEP_LIMIT_MIN ? "Min" : \
+ (X) == DISC_TIMESTEP_LIMIT_MAX ? "Max" : \
+ (X) == DISC_TIMESTEP_LIMIT_TOO_FAST ? "Too fast" : \
+ "Unknown" \
+ )
+
+
+#endif // DISC_FUNCTION_MACROS_H
diff --git a/src/disc/disc_initialize_disc.c b/src/disc/disc_initialize_disc.c
index b8db91fdee634f32614068af9685ff5650ead3a9..ffa918939e7c396a9d59c714eff61ab229a7e6f0 100644
--- a/src/disc/disc_initialize_disc.c
+++ b/src/disc/disc_initialize_disc.c
@@ -70,10 +70,8 @@ void disc_initialize_disc(struct stardata_t * const stardata,
disc->prevzone = DISC_NO_PREVZONE_YET;
#endif
#ifdef MEMOIZE
- printf("disc memoize init was %p\n",disc->memo);
disc->memo = NULL;
memoize_initialize(&disc->memo);
- printf("disc memoize init %p\n",disc->memo);
#endif
/*
diff --git a/src/disc/disc_macros.h b/src/disc/disc_macros.h
index d9052ce9ec0d9c1decc4927b9e42ab8ac3181b80..670ba75a37589acd3b0300f79cb5d8273d339b2f 100644
--- a/src/disc/disc_macros.h
+++ b/src/disc/disc_macros.h
@@ -3,8 +3,7 @@
#define DISC_MACROS_H
/*
- * Macros related to the binary_c treatment of circumstellar and
- * circumbinary discs
+ * Macros related to the binary_c treatment of discs
*/
/* disc types (unsigned int) */
@@ -17,21 +16,6 @@
#define DISC_INNER_EDGE -1
#define DISC_OUTER_EDGE -2
-/* convert position to disc index for an existing disc */
-#define DISC_POS_TO_N(OBJ,POS) \
- ((POS) == DISC_INNER_EDGE ? 0 : \
- (POS) == DISC_OUTER_EDGE ? (OBJ)->ndiscs-1 : \
- (POS))
-
-/* convert position to disc index for a new disc */
-#define DISC_NEWPOS_TO_N(OBJ,POS) \
- ((POS) == DISC_INNER_EDGE ? 0 : \
- (POS) == DISC_OUTER_EDGE ? (OBJ)->ndiscs : \
- (POS))
-
-#define Discloop(OBJ,A) \
- for((A)=0;(A)<(OBJ)->ndiscs;(A)++)
-
#define BISECT_M 1
#define BISECT_J 2
@@ -41,39 +25,6 @@
#define CONVERGENCE_SUCCEEDED 1
#define CONVERGENCE_FLUX_NEGATIVE 2
-#define Convergence_status(S) \
- ( \
- (S) == CONVERGENCE_FAILED ? "Failed" : \
- (S) == CONVERGENCE_SUCCEEDED ? "Success" : \
- (S) == CONVERGENCE_FLUX_NEGATIVE ? "Flux negative" : \
- "Unknown" \
- )
-
-
-/*
- * Macro to determine whether a disc is actually a
- * disc. It cannot be if either M or J is zero,
- * or if its lifetime exceeds DISC_MAX_LIFETIME.
- */
-
-#define IS_NONZERO_NUMBER(X) \
- ((Is_really_not_zero(X) && \
- !isnan(X) && !isinf(X)))
-
-/*
- * Macro to determine whether a disc "really exists".
- *
- * This can often be used to determine whether other
- * functions, like integrals over the disc, will work.
- */
-#define Disc_is_disc(D) \
- ( \
- IS_NONZERO_NUMBER((D)->M) \
- && \
- IS_NONZERO_NUMBER((D)->J) \
- && \
- (D)->lifetime < DISC_MAX_LIFETIME \
- )
/*
* Reasons for disc removal
@@ -83,24 +34,6 @@
#define DISC_REMOVE_EVAPORATE 2
#define DISC_REMOVE_ALL 3
-#define Disc_removal_string(S) \
- ( \
- (S)==DISC_REMOVE_FORMED_RING ? "Formed ring" : \
- (S)==DISC_REMOVE_REPLACE ? "Replaced" : \
- (S)==DISC_REMOVE_EVAPORATE ? "Evaporated" : \
- (S)==DISC_REMOVE_ALL ? "Remove all" : \
- "Unknown" \
- )
-
-
-
-/* solvers */
-#define Solver_string(N) ( \
- (N)==(DISC_SOLVER_BISECTION) ? "Bi" : \
- (N)==(GSL_MULTIROOT_FSOLVER_HYBRIDS) ? "Hs" : \
- (N)==(GSL_MULTIROOT_FSOLVER_HYBRID) ? "H" : \
- "?" \
- )
#define DISC_SOLVER_ORDERED_LIST \
GSL_MULTIROOT_FSOLVER_HYBRIDS, \
@@ -113,23 +46,6 @@
#define DISC_SOLVER_NONE 0
#define DISC_SOLVER_BISECTION -1
-#define Show_disc_derivative(N) \
- Discdebug(1, \
- "d(%20s)/dt = M % 20g (Msun/y), J % 20g (cgs), e % 20g (s^-1)\n", \
- Disc_derivative_string(N), \
- disc->loss[(N)].mdot / M_SUN * YEAR_LENGTH_IN_SECONDS, \
- disc->loss[(N)].jdot, \
- disc->loss[(N)].edot);
-
-
-#define Clear_disc_feedback(D) \
- { \
- (D)->dM_ejected = \
- (D)->dM_binary = \
- (D)->dJ_binary = \
- (D)->de_binary = 0.0; \
- }
-
/*
* Circumbinary disc eccentricity pumping algorithms
@@ -191,26 +107,6 @@
"Inner edge", \
"Outer edge"
-#define Disc_derivative_string(A) ( \
- (A) == DISC_LOSS_GLOBAL ? "Global" : \
- (A) == DISC_LOSS_INNER_VISCOUS ? "Viscous" : \
- (A) == DISC_LOSS_INNER_L2_CROSSING ? "L2 crossing": \
- (A) == DISC_LOSS_FUV ? "FUV" : \
- (A) == DISC_LOSS_XRAY ? "X-ray" : \
- (A) == DISC_LOSS_BINARY_TORQUE ? "Binary torque" : \
- (A) == DISC_LOSS_ISM ? "ISM ram stripping" : \
- (A) == DISC_LOSS_RESONANCES ? "Resonances" : \
- (A) == DISC_LOSS_NEGATIVE_FLUX ? "Negative flux" : \
- (A) == DISC_LOSS_INNER_EDGE ? "Inner edge" : \
- (A) == DISC_LOSS_OUTER_EDGE ? "Outer edge" : \
- "Unknown" )
-
-/* macro to define mass loss from the inner or outer edge */
-#define Disc_loss_from_edge(I) \
- ((I)==DISC_LOSS_INNER_EDGE || \
- (I)==DISC_LOSS_OUTER_EDGE)
-
-
#define DISC_ZONES_FAILED 0
#define DISC_ZONES_OK 1
#define DISC_ZONES_NEW_ZONE_LIST_RIN_EQUALS_ROUT 2
@@ -229,14 +125,6 @@
#define DISC_LOG_EVOLVE_DISC_STRUCTURE2 1
#define DISC_LOG_EVERY_TIMESTEP 2
-/*
- * ... and a string to log them.
- */
-#define Disc_logger(A) ( \
- (A) == DISC_LOG_EVOLVE_DISC_STRUCTURE ? "disc_evolve_disc_structure" : \
- (A) == DISC_LOG_EVOLVE_DISC_STRUCTURE2 ? "disc_evolve_disc_structure (2)" : \
- (A) == DISC_LOG_EVERY_TIMESTEP ? "log_every_timestep" : \
- "Unknown")
/*
@@ -283,35 +171,11 @@
#define DISC_EDGE_STRIPPING_BEFORE_CONVERGENCE_STEP 2
-/*
- * Copy disc from source S to destination D
- */
-#define Copy_disc(S,D) memcpy((D),(S),sizeof(struct disc_t));
-
/*
* Space for a new disc
*/
#define New_disc Malloc(sizeof(struct disc_t))
-/*
- * New disc copied from another
- */
-#define New_disc_from(D) Copy_disc((D),New_disc)
-
-
-
-
-#define Evaporate_disc(D,REASON) \
- (D)->M = 0.0; \
- (D)->J = 0.0; \
- (D)->dt = LONG_TIMESTEP; \
- (D)->converged = FALSE; \
- if(DISC_DEBUG) \
- printf("Disc evaporated at line %d in file %s because %s\n", \
- __LINE__, \
- __FILE__, \
- REASON); \
- /* disc evaporated */
/*
* Disc timestep limiters
@@ -322,16 +186,7 @@
#define DISC_TIMESTEP_LIMIT_MIN 3
#define DISC_TIMESTEP_LIMIT_MAX 4
#define DISC_TIMESTEP_LIMIT_TOO_FAST 5
-#define Disc_timestep_limit_string(X) \
- ( \
- (X) == DISC_TIMESTEP_LIMIT_M ? "M" : \
- (X) == DISC_TIMESTEP_LIMIT_J ? "J" : \
- (X) == DISC_TIMESTEP_LIMIT_F ? "F" : \
- (X) == DISC_TIMESTEP_LIMIT_MIN ? "Min" : \
- (X) == DISC_TIMESTEP_LIMIT_MAX ? "Max" : \
- (X) == DISC_TIMESTEP_LIMIT_TOO_FAST ? "Too fast" : \
- "Unknown" \
- )
+
/*
* disc log levels
*/
diff --git a/src/disc/disc_monotonic_check.c b/src/disc/disc_monotonic_check.c
index a6a855efae61af440606a63039c9356dd66cb068..d04ea42d4435dfdb668fb30f8dcaa3aaeb641d6b 100644
--- a/src/disc/disc_monotonic_check.c
+++ b/src/disc/disc_monotonic_check.c
@@ -9,14 +9,14 @@ Boolean disc_monotonic_check(struct stardata_t * const stardata,
)
{
/*
- * Check that the disc M,J,F change monotonically in
- * the disc parameters
+ * Check that the disc's M,J,F change monotonically in
+ * the disc parameters.
*/
- double unknown MAYBE_UNUSED;
+ double unknown MAYBE_UNUSED = 0.0;
unsigned int i;
/*
- * Allocate memory
+ * Allocate memory for the checks
*/
double *ref = Malloc(sizeof(double)*n);
disc_parameter * parameters = Malloc(sizeof(disc_parameter)*n);
diff --git a/src/disc/disc_new_zone_radii.c b/src/disc/disc_new_zone_radii.c
index 21398d80a58b87198d426cc2aa4298d3840b6cb4..32db3f7b297e6917fd7170d81ece0f63668b2d0c 100644
--- a/src/disc/disc_new_zone_radii.c
+++ b/src/disc/disc_new_zone_radii.c
@@ -9,12 +9,12 @@
Boolean disc_new_zone_radii(struct disc_t * const disc)
{
/*
- * Determine disc zone locations and their radii
+ * Determine disc zone locations and their radii.
*
* Note that this function assumes that the temperature
- * functions are of the form:
+ * functions are of the form,
*
- * T(R) = T_0 R^k
+ * T(R) = T_0 R^k,
*
* where k<0.
*
@@ -25,13 +25,41 @@ Boolean disc_new_zone_radii(struct disc_t * const disc)
* On input to this function, disc->thermal_zones contains
* all possible zones.
*
- * First we make a copy of these zones, in an array P,
- * and set all their inner radii to the disc reference length (1cm),
- * and make them all "valid".
+ * We use two arrays of thermal zones, P and Q.
+ *
+ * First we copy the disc's zones to P, set all their inner
+ * radii to the disc reference length (1cm), and make them
+ * all "valid".
+ *
+ * We then set the first zone in Q[0] : this is the hottest
+ * at the reference length. Then we find the next crossing
+ * point which defines zone Q[1], then continue to Q[2].
+ *
+ * Next the zones are trimmed according to the disc size,
+ * and finally set into the disc.
+ */
+ Zprintf("\n\n\nNew zones %g %s %g %s; %g %s %g %s\n\n\n",
+ disc->Rin/R_SUN,
+ R_SOLAR,
+ disc->Rout/R_SUN,
+ R_SOLAR,
+ disc->M/M_SUN,
+ M_SOLAR,
+ disc->J,
+ J_CGS);
+
+ /* allocate memory for P and Q */
+ struct disc_thermal_zone_t * RESTRICT P =
+ Malloc(sizeof(struct disc_thermal_zone_t)*
+ (DISCS_MAX_N_ZONES));
+ struct disc_thermal_zone_t * RESTRICT Q =
+ Malloc(sizeof(struct disc_thermal_zone_t)*
+ (DISCS_MAX_N_ZONES));
+
+ /*
+ * Copy the disc's zones to P and set all P's
+ * zones to be valid.
*/
- Zprintf("\n\n\nNew zones %g %g; %g %g\n\n\n",disc->Rin,disc->Rout,disc->M,disc->J);
- struct disc_thermal_zone_t * P = Malloc(sizeof(struct disc_thermal_zone_t)*
- (DISCS_MAX_N_ZONES));
Disc_zone_counter i;
for(i=0;i<DISCS_MAX_N_ZONES;i++)
{
@@ -43,8 +71,6 @@ Boolean disc_new_zone_radii(struct disc_t * const disc)
}
Zprintf("\n\n\n");
- struct disc_thermal_zone_t * Q = Malloc(sizeof(struct disc_thermal_zone_t)*
- (DISCS_MAX_N_ZONES));
/*
* Now we use the variable r as the inner edge of the "next"
@@ -54,14 +80,36 @@ Boolean disc_new_zone_radii(struct disc_t * const disc)
*/
double r = DISC_REFERENCE_LENGTH;
-
/*
* Counters: the number of remaining and allocated zones
*/
Disc_zone_counter zones_remaining = DISCS_MAX_N_ZONES;
Disc_zone_counter zones_allocated = 0;
-
+#ifdef ___DEBUGGING
+ if(0)
+ {
+ /*
+ * Debugging:
+ * Show the crossing radius matrix
+ */
+ int i,j;
+ for(i=0;i<DISCS_MAX_N_ZONES;i++)
+ {
+ for(j=0;j<=i;j++)
+ {
+ const double Rx = disc_Rcross(&P[i],
+ &P[j]);
+ printf("%d x %d : %g %s\n",
+ i,
+ j,
+ Rx/R_SUN,
+ R_SOLAR);
+ }
+ }
+ }
+#endif// ___DEBUGGING
+
/*
* Loop until all zones have been allocated
*/
@@ -72,10 +120,9 @@ Boolean disc_new_zone_radii(struct disc_t * const disc)
* Find the hottest valid zone at this radius.
*
* This is zone ihottest with temperature
- * Thottest
+ * Thottest.
*
- * Zones that have already been used are already
- * invalid so cannot be repeated.
+ * Zones that have already been used are ignored.
*/
double Thottest = 0.0;
Disc_zone_counter ihottest = -1;
@@ -83,38 +130,46 @@ Boolean disc_new_zone_radii(struct disc_t * const disc)
{
if(P[i].valid == TRUE)
{
- double T = power_law(&P[i].Tlaw,r);
+ const double T = power_law(&P[i].Tlaw,r);
if(T > Thottest)
{
ihottest = i;
Thottest = T;
}
- Zprintf("Zone %u : T = %g (A=%g exp=%g) : current hottest %u\n",
- i,
- T,
- P[i].Tlaw.A0,
- P[i].Tlaw.exponent,
- ihottest);
+ Zprintf("At r = %30.15e %s, Zone %u : T = %g (A=%g exp=%g) : current hottest %u\n",
+ r/R_SUN,
+ R_SOLAR,
+ i,
+ T,
+ P[i].Tlaw.A0,
+ P[i].Tlaw.exponent,
+ ihottest);
}
}
- Zprintf("Hottest %u\n",ihottest);
+ Zprintf("Hottest at r = %30.15e %s is zone %u\n",
+ r/R_SUN,
+ R_SOLAR,
+ ihottest);
/*
- * If ihottest is -1 we're in trouble
+ * If ihottest is -1 we're in trouble:
+ * this should not happen!
*/
if(ihottest == -1)
{
Safe_free(P);
Safe_free(Q);
- printf("ihottest -1 error\n");Flexit;
+ printf("ihottest -1 error\n");
return FALSE;
}
/*
- * Zone ihottest, in P, is the hottest at r
+ * Zone P[ihottest] is the hottest zone at radius r
*
- * So copy this onto the new zone list, Q.
+ * So copy this onto the new zone list, Q, at
+ * its next position, and set P[ihottest] to be invalid
+ * so it cannot be set into Q twice.
*/
Zprintf("Set zone P[%u] in Q[%u]\n",ihottest,zones_allocated);
Copy_zone(&P[ihottest],
@@ -128,42 +183,56 @@ Boolean disc_new_zone_radii(struct disc_t * const disc)
* we cross into the next zone.
*
* If there are no zones remaining, just set
- * smallest_Rx to be DISC_LARGE_RADIUS,
+ * next_crossing_radius to be DISC_LARGE_RADIUS,
* and leave found_boundary as FALSE.
*/
-
Boolean found_boundary = FALSE;
- double smallest_Rx = DISC_LARGE_RADIUS;
- Disc_zone_counter MAYBE_UNUSED smallest_ix;
+ double next_crossing_radius = DISC_LARGE_RADIUS;
+ Disc_zone_counter MAYBE_UNUSED next_crossing_radius_index;
if(zones_remaining > 0)
- {
+ {
for(i=0;i<DISCS_MAX_N_ZONES;i++)
{
if(i != ihottest &&
P[i].valid == TRUE)
{
- /*
- Zprintf("Call RCross %u (%g %g %g) %u (%g %g %g)\n",
- ihottest,
- P[ihottest].Tlaw.A0,
- P[ihottest].Tlaw.exponent,
- P[ihottest].Tlaw.R0,
- P[i].Tlaw.A0,
- P[i].Tlaw.exponent,
- P[i].Tlaw.R0
+ Zprintf("Call RCross %u (%g %g %g %s) %u (%g %g %g %s)\n",
+ ihottest,
+ P[ihottest].Tlaw.A0,
+ P[ihottest].Tlaw.exponent,
+ P[ihottest].Tlaw.R0/R_SUN,
+ R_SOLAR,
+ i,
+ P[i].Tlaw.A0,
+ P[i].Tlaw.exponent,
+ P[i].Tlaw.R0/R_SUN,
+ R_SOLAR
);
- */
- double Rx = disc_Rcross(&P[ihottest],
- &P[i]);
- if(Rx < smallest_Rx)
+
+ /*
+ * Find the crossing point of the hottest
+ * zone at r with the other zones
+ */
+ const double Rx = disc_Rcross(&P[ihottest],
+ &P[i]);
+
+ /*
+ * If this crossing radius, Rx, is outside the
+ * current radius but smaller than next_crossing_radius,
+ * replace next_crossing_radius with Rx.
+ */
+ if(Rx > r &&
+ Rx < next_crossing_radius)
{
- smallest_Rx = Rx;
- smallest_ix = i;
+ next_crossing_radius = Rx;
+ next_crossing_radius_index = i;
found_boundary = TRUE;
}
- Zprintf("Zone %u : Rx at %g, current smallest %u\n",
- i,Rx,smallest_ix);
+ Zprintf("Zone %u : Rx at %30.15e %s\n",
+ i,
+ Rx/R_SUN,
+ R_SOLAR);
}
}
}
@@ -171,21 +240,26 @@ Boolean disc_new_zone_radii(struct disc_t * const disc)
/*
* The end of the zone is at the crossing point.
*
- * We should also set r here because this is the new test point.
+ * We should also set r here because this is the new test point
+ * on the next iteration of the main loop.
*/
- r = Q[zones_allocated].rend = smallest_Rx;
- Zprintf("Set zone %u rend = %g\n",zones_allocated,Q[zones_allocated].rend);
+ r = Q[zones_allocated].rend = next_crossing_radius;
+ Zprintf("This zone (%u) ends at next crossing radius = %30.15e %s\n",
+ zones_allocated,
+ Q[zones_allocated].rend/R_SUN,
+ R_SOLAR);
/*
- * Make zone valid
+ * Make new zone in Q valid
*/
Q[zones_allocated].valid = TRUE;
if(found_boundary == TRUE)
{
- Zprintf("Next crossing with zone %u, R = %g\n",
- smallest_ix,
- smallest_Rx);
+ Zprintf("Next crossing with zone %u, R = %30.15e %s\n",
+ next_crossing_radius_index,
+ next_crossing_radius/R_SUN,
+ R_SOLAR);
}
else
{
@@ -199,26 +273,33 @@ Boolean disc_new_zone_radii(struct disc_t * const disc)
}
/*
- * See starting radii
+ * Set starting radii
*/
-
Q[0].rstart = DISC_REFERENCE_LENGTH;
for(i=1;i<DISCS_MAX_N_ZONES;i++)
{
Q[i].rstart = Q[i-1].rend;
- Zprintf("Set zone %u rstart = %g\n",i,Q[i].rstart);
+ Zprintf("Set zone %u rstart = %30.15e %s\n",
+ i,
+ Q[i].rstart/R_SUN,
+ R_SOLAR);
}
for(i=0;i<DISCS_MAX_N_ZONES;i++)
{
- Zprintf("ZONE %u : rstart = %g to rend = %g\n",
+ const struct disc_thermal_zone_t * MAYBE_UNUSED z = Q + i;
+ Zprintf("ZONE %u %s %s : rstart = %30.15e %s to rend = %30.15e %s\n",
i,
- Q[i].rstart,
- Q[i].rend);
+ DISC_ZONE_TYPE_STRING(z),
+ DISC_ZONE_VALID_STRING(z),
+ z->rstart/R_SUN,
+ R_SOLAR,
+ z->rend/R_SUN,
+ R_SOLAR);
}
/*
- * make sure all radii are in the disc
+ * Trim the zones in Q to make sure they are inside the disc
*/
for(i=0;i<DISCS_MAX_N_ZONES;i++)
{
@@ -242,18 +323,40 @@ Boolean disc_new_zone_radii(struct disc_t * const disc)
nzones--;
}
}
-
+
/*
- * Copy valid thermal zones from the new, ordered zone list Q
+ * Copy valid thermal zones from Q
* back into the disc
- */
-
+ */
for(i=0;i<nzones;i++)
{
Copy_zone(Q+i,&disc->thermal_zones[i]);
}
disc->n_thermal_zones = nzones;
-
+
+ /*
+ * Check zone radii
+ */
+ for(i=0;i<disc->n_thermal_zones;i++)
+ {
+ const struct disc_thermal_zone_t * const z = disc->thermal_zones + i;
+ if(z->rstart > z->rend)
+ {
+ fprintf(stderr,
+ "Disc zone %u has rstart = %g %s > rend = %g %s\n",
+ i,
+ z->rstart/R_SUN,
+ R_SOLAR,
+ z->rend/R_SUN,
+ R_SOLAR
+ );
+ Safe_free(P);
+ Safe_free(Q);
+ return FALSE;
+ }
+ }
+
+
/*
* Set up the temperature power laws
*/
diff --git a/src/disc/disc_nsector.c b/src/disc/disc_nsector.c
index c036b1b6328024cb27ae402c709caefab8a5fbaf..b8b026644d432e9251cb914d0d44b6e485f413e2 100644
--- a/src/disc/disc_nsector.c
+++ b/src/disc/disc_nsector.c
@@ -51,11 +51,6 @@ void disc_nsector(const unsigned int n,
double unknown MAYBE_UNUSED;
struct disc_t * discwas;
unsigned int i,j;
-
- printf("conv in disc_nsector Tvisc0 = %g Rin = %g Rout = %g\n",
- disc->Tvisc0,
- disc->Rin,
- disc->Rout);
/*
* Back up disc
diff --git a/src/disc/disc_parameters.h b/src/disc/disc_parameters.h
index 6a4f07f0438536b55b7348d9c4d8194e9eb35575..60c78dc59ece28afab8eaae04231e2c041a93bfb 100644
--- a/src/disc/disc_parameters.h
+++ b/src/disc/disc_parameters.h
@@ -384,6 +384,14 @@
*/
#define DISC_MONTE_CARLO_GUESSES_USE_LOG
+
+/*
+ * Use log(x) instead of (x) in nsection? This is
+ * slightly quicker and prevents negative values
+ * being used.
+ */
+#define DISC_NSECTION_USELOG
+
/*
* Use logarithm(x) instead of (x) in bisection?
*
diff --git a/src/disc/disc_setup_convergence_parameter_types_and_constraints.c b/src/disc/disc_setup_convergence_parameter_types_and_constraints.c
index 7a3cec802bf9088d146353727ccbc03853ff515d..b6766dca2ea315ff6b0c4d6561f52e63c183c351 100644
--- a/src/disc/disc_setup_convergence_parameter_types_and_constraints.c
+++ b/src/disc/disc_setup_convergence_parameter_types_and_constraints.c
@@ -4,8 +4,7 @@
#include "disc_constraints.h"
#include "disc_boundary_conditions.h"
-#undef DISC_DEBUG
-#define DISC_DEBUG 2
+
/*
* Macro to set up of parameter_types, constraints and boundary
* conditions, and update the _p pointers
diff --git a/src/disc/disc_show_thermal_zones.c b/src/disc/disc_show_thermal_zones.c
index 0480de6cbc298ae12ed3aa5de468ad900404e8c4..d94ce8300bb74b3370cedf76d62360e992a724df 100644
--- a/src/disc/disc_show_thermal_zones.c
+++ b/src/disc/disc_show_thermal_zones.c
@@ -8,6 +8,11 @@ void disc_show_thermal_zones(struct disc_t * disc,
const struct disc_thermal_zone_t * zones,
const int n)
{
+ if(zones == NULL)
+ {
+ zones = disc->thermal_zones;
+ }
+
const int max = n==-1 ? DISCS_MAX_N_ZONES : n;
int i;
const Boolean show_invalid = FALSE;
diff --git a/src/disc/disc_thermal_zones.h b/src/disc/disc_thermal_zones.h
index eea365f40eea83952cddce9e296bdb2482148083..26d30661de174d9e73fbfe86fdea6beb7f97cd7c 100644
--- a/src/disc/disc_thermal_zones.h
+++ b/src/disc/disc_thermal_zones.h
@@ -66,10 +66,16 @@
#define DISC_ZONE_TYPE_RADIATIVE_OUTER 2
#define DISC_ZONE_TYPE_UNDEFINED -1
-#define DISC_ZONE_TYPE_STRING(Z) ((Z->type)==DISC_ZONE_TYPE_VISCOUS ? "Viscous" : \
- ((Z->type)==DISC_ZONE_TYPE_RADIATIVE_INNER ? "Rad in " : "Rad out"))
+#define DISC_ZONE_TYPE_STRING(Z) \
+ ( \
+ (Z->type)==DISC_ZONE_TYPE_VISCOUS ? "Viscous" : \
+ (Z->type)==DISC_ZONE_TYPE_RADIATIVE_INNER ? "Rad in " : \
+ (Z->type)==DISC_ZONE_TYPE_RADIATIVE_OUTER ? "Rad out" : \
+ "undefined" \
+ )
-#define DISC_ZONE_VALID_STRING(Z) ((Z)->valid ? " valid" : "invalid")
+#define DISC_ZONE_VALID_STRING(Z) \
+ ((Z)->valid ? " valid" : "invalid")
/*
* Temperature power law definitions
diff --git a/src/disc/inner_disc.c b/src/disc/inner_disc.c
index 186fe00bb90268a172d505f120fee029da7e3928..d307e592da63b5051bc41e847e1ee9090d56b6d5 100644
--- a/src/disc/inner_disc.c
+++ b/src/disc/inner_disc.c
@@ -13,12 +13,12 @@ struct disc_t * inner_disc(const Disc_type type,
if(type == DISC_CIRCUMSTELLAR)
{
struct star_t * star = (struct star_t *) object;
- disc = star->discs[DISC_POS_TO_N(star,DISC_INNER_EDGE)];
+ disc = star->discs[Disc_pos_to_n(star,DISC_INNER_EDGE)];
}
else if(type == DISC_CIRCUMBINARY)
{
struct stardata_t * stardata = (struct stardata_t *) object;
- disc = &stardata->common.discs[DISC_POS_TO_N(&(stardata->common),
+ disc = &stardata->common.discs[Disc_pos_to_n(&(stardata->common),
DISC_INNER_EDGE)];
}
else
diff --git a/src/disc/new_disc.c b/src/disc/new_disc.c
index 87b2bd134e90817307290caa14dacfc48d112197..53e88a85300b1dc9ff26d91e3b1d1fb95c1393e1 100644
--- a/src/disc/new_disc.c
+++ b/src/disc/new_disc.c
@@ -42,7 +42,7 @@ struct disc_t * new_disc(
BINARY_C_ALLOC_FAILED,
"Memory allocation (Realloc) of a new circumstellar disc pointer failed");
}
- int n = DISC_NEWPOS_TO_N(star,pos);
+ int n = Disc_newpos_to_n(star,pos);
int i;
for(i=star->ndiscs-1; i>=n; i--)
{
@@ -64,7 +64,7 @@ struct disc_t * new_disc(
{
struct stardata_t * stardata = object;
struct common_t * common = &(stardata->common);
- int n = DISC_NEWPOS_TO_N(common,pos);
+ int n = Disc_newpos_to_n(common,pos);
if(n < NDISCS)
{
struct disc_t * disc = &common->discs[n];
diff --git a/src/disc/outer_disc.c b/src/disc/outer_disc.c
index 2d36945984c7d179d86dce556ef4a0271e3b087a..0116977bdac33400283da163eef3708c44f6a566 100644
--- a/src/disc/outer_disc.c
+++ b/src/disc/outer_disc.c
@@ -13,12 +13,12 @@ struct disc_t * outer_disc(const Disc_type type,
if(type == DISC_CIRCUMSTELLAR)
{
struct star_t * star = (struct star_t *) object;
- disc = star->discs[DISC_POS_TO_N(star,DISC_OUTER_EDGE)];
+ disc = star->discs[Disc_pos_to_n(star,DISC_OUTER_EDGE)];
}
else if(type == DISC_CIRCUMBINARY)
{
struct stardata_t * stardata = (struct stardata_t *) object;
- disc = &stardata->common.discs[DISC_POS_TO_N(&(stardata->common),DISC_OUTER_EDGE)];
+ disc = &stardata->common.discs[Disc_pos_to_n(&(stardata->common),DISC_OUTER_EDGE)];
}
else
{
diff --git a/src/disc/remove_disc.c b/src/disc/remove_disc.c
index b64bf660125808fb19434b4490c9ecfee5b18e0e..8022b121dea7976a1db731780151a0bd5307331d 100644
--- a/src/disc/remove_disc.c
+++ b/src/disc/remove_disc.c
@@ -29,11 +29,9 @@ void remove_disc(
BINARY_C_OUT_OF_RANGE,
"No circumstellar disc to remove in remove_disc\n");
}
- int n = DISC_POS_TO_N(star,pos);
+ int n = Disc_pos_to_n(star,pos);
#ifdef MEMOIZE
- printf("remove_disc memoize_free %p\n",star->discs[n]->memo);fflush(NULL);
memoize_free(&star->discs[n]->memo);
- printf("remove_disc now %p\n",star->discs[n]->memo);fflush(NULL);
#endif // MEMOIZE
Safe_free(star->discs[n]);
int i;
@@ -72,7 +70,7 @@ void remove_disc(
}
else
{
- int n = DISC_POS_TO_N(common,pos);
+ int n = Disc_pos_to_n(common,pos);
if(DISC_DEBUG)
{
@@ -109,9 +107,7 @@ void remove_disc(
if(disc->memo != NULL)
{
Dprint("call memoize_free\n");
- printf("remove_disc memoize_free %p\n",disc->memo);fflush(NULL);
memoize_free(&disc->memo);
- printf("remove_disc now %p\n",disc->memo);fflush(NULL);
}
#endif//MEMOIZE
memset(disc,0,sizeof(struct disc_t));
diff --git a/src/maths/generic_bisect.c b/src/maths/generic_bisect.c
index 764301a69d025f85c4ee53eadecbab2b9fc9adc6..a82d54b9bc199c57fdfb85eb14961f6d8d95b3a8 100644
--- a/src/maths/generic_bisect.c
+++ b/src/maths/generic_bisect.c
@@ -90,7 +90,7 @@ double generic_bisect(int * RESTRICT error,
max = log10(Max(REALLY_TINY,max));
guess = log10(guess);
}
- int vb = 1;
+ Boolean vb = FALSE;
*error = BINARY_C_BISECT_ERROR_NONE;
va_list args_master,args;
double top,bottom;
@@ -108,11 +108,11 @@ double generic_bisect(int * RESTRICT error,
* we should use the equivalent positive identifier
* and turn verbosity on.
*/
- vb = 1;
+ vb = TRUE;
n = -n;
}
- //if(n==BISECTOR_DISC_BISECTION_ROOTER) vb = 1;
+ //if(n==BISECTOR_DISC_BISECTION_ROOTER) vb = TRUE;
/*
* call va_start to get args_master
diff --git a/src/maths/generic_bisect.h b/src/maths/generic_bisect.h
index f9374a2b346e00f1bcb3c1ac2930f0a6ada9dc5a..15612bea8efadee4b4cacd19cde2a7117c81d0cb 100644
--- a/src/maths/generic_bisect.h
+++ b/src/maths/generic_bisect.h
@@ -63,7 +63,7 @@
/* allow verbosity */
#define VBPRINT(...) \
- if(vb) \
+ if(vb == TRUE) \
{ \
generic_bisect_vbprint(n,__VA_ARGS__); \
}
diff --git a/tbse b/tbse
index 51f36d68a5da0fd6205f8baec1ff3ca6bca1e9b0..44485787d85003167450dacc257c84403176685a 100755
--- a/tbse
+++ b/tbse
@@ -594,7 +594,7 @@ DISC_N_MONTE_CARLO_GUESSES=0
CBDISC_GAMMA=1.4 # 1.4 (molecular hydrogen) or 1.6666666 (monatomic gas)
CBDISC_ALPHA=1e-3 # 1e-2 or 1e-3
-CBDISC_KAPPA=1e-3 # cgs (cm^2 g^-1) : 1e-2 for cool gas, 1e5 for dust in the visible cf https://ned.ipac.caltech.edu/level5/Sept05/Li2/Li4.html
+CBDISC_KAPPA=1e-2 # cgs (cm^2 g^-1) : 1e-2 for cool gas, 1e5 for dust in the visible cf https://ned.ipac.caltech.edu/level5/Sept05/Li2/Li4.html
CBDISC_TORQUEF=1e-4 # 1e-3 tends to keep discs outside L2
##############################
@@ -678,7 +678,7 @@ CBDISC_MASS_LOSS_INNER_VISCOUS_ACCRETION_METHOD=CBDISC_MASS_LOSS_INNER_VISCOUS_A
CBDISC_MASS_LOSS_INNER_L2_CROSS_MULTIPLIER=1.0
# couple L2 and viscous losses at the inner edge (True)
-CBDISC_VISCOUS_L2_COUPLING=False
+CBDISC_VISCOUS_L2_COUPLING=True
# ISM pressure in units of the Boltzmann constant/K (3000) e.g.
# http://www.astronomy.ohio-state.edu/~pogge/Ast871/Notes/Intro.pdf page 15
@@ -698,8 +698,8 @@ CBDISC_MASS_LOSS_XRAY_MULTIPLIER=1.0
CBDISC_VISCOUS_PHOTOEVAPORATION_COUPLING=True
# edge stripping (True/False , both true by default)
-CBDISC_INNER_EDGE_STRIPPING=F
-CBDISC_OUTER_EDGE_STRIPPING=F
+CBDISC_INNER_EDGE_STRIPPING=True
+CBDISC_OUTER_EDGE_STRIPPING=True
# timescales for edge stripping:
# DISC_STRIPPING_TIMESCALE_INSTANT = 0 = instant (default)