diff --git a/src/disc/disc_bisection_rooter.c b/src/disc/disc_bisection_rooter.c
index bb3c1f59ab08d52ce3f54f04966b9a0ef8a17826..777cef8072725b337d665df4ef833d6985b9004b 100644
--- a/src/disc/disc_bisection_rooter.c
+++ b/src/disc/disc_bisection_rooter.c
@@ -192,7 +192,7 @@ int disc_bisection_rooter(const int n,
}
else if(local_constraints[0] == DISC_CONSTRAINT_J_TO_DISC_J)
{
- DISC_OUTER_ZONE(disc)->rend = disc->Rout;
+ Disc_outer_zone(disc)->rend = disc->Rout;
}
diff --git a/src/disc/disc_extend_powerlaw_from_reference_point.c b/src/disc/disc_extend_powerlaw_from_reference_point.c
index f92941e72e7a2c286cc18d935516fcfb54ceecad..c9e1fac449195e9e7c38fe603c12ef7b2ac9580a 100644
--- a/src/disc/disc_extend_powerlaw_from_reference_point.c
+++ b/src/disc/disc_extend_powerlaw_from_reference_point.c
@@ -26,9 +26,9 @@ void disc_extend_powerlaw_from_reference_point(struct disc_t * const disc,
* Note that this zone IS a valid zone and should be the outermost
* zone.
*/
- int iref = DISC_NEAREST_ZONE_N(Rref,disc);
- struct disc_thermal_zone_t * zref = DISC_NEAREST_ZONE(Rref,disc);
- int n_nearest = DISC_NEAREST_ZONE_N(Rref,disc);
+ int iref = Disc_nearest_zone_n(Rref,disc);
+ struct disc_thermal_zone_t * zref = Disc_nearest_zone(Rref,disc);
+ int n_nearest = Disc_nearest_zone_n(Rref,disc);
/*
* Hence A0 and A1 in this zone
diff --git a/src/disc/disc_function_macros.h b/src/disc/disc_function_macros.h
index a8dad9c73fa7d0d06935a7f2f92b07854660fbc5..2062024c9a2da8b62190653f8d59aa460fcb7677 100644
--- a/src/disc/disc_function_macros.h
+++ b/src/disc/disc_function_macros.h
@@ -151,4 +151,8 @@
)
+#define Disc_Mdot(X) (Max(0.0,-disc->loss[(X)].mdot / M_SUN * YEAR_LENGTH_IN_SECONDS))
+#define Disc_Jdot(X) (Max(0.0,-disc->loss[(X)].jdot))
+
+
#endif // DISC_FUNCTION_MACROS_H
diff --git a/src/disc/disc_hottest_zone_at_radius.c b/src/disc/disc_hottest_zone_at_radius.c
new file mode 100644
index 0000000000000000000000000000000000000000..ea141b7698267b975053d4183db937dd1eb30282
--- /dev/null
+++ b/src/disc/disc_hottest_zone_at_radius.c
@@ -0,0 +1,42 @@
+#include "../binary_c.h"
+
+#ifdef DISCS
+
+Disc_zone_counter disc_hottest_zone_at_radius(struct disc_thermal_zone_t * zones,
+ const double radius,
+ const Boolean all,
+ double * Thottest_p)
+{
+ /*
+ * Return the counter of the hottest zone at
+ * given radius, and set *Thottest_p to this temperature
+ * if Thottest_p != NULL.
+ *
+ * Returns -1 on failure.
+ *
+ * Invalid zones are skipped, unless all == TRUE.
+ */
+ double Thottest = 0.0;
+ Disc_zone_counter i,ihottest = -1;
+
+ for(i=0;i<DISCS_MAX_N_ZONES;i++)
+ {
+ if(all == TRUE ||
+ zones[i].valid == TRUE)
+ {
+ const double T = power_law(&zones[i].Tlaw,radius);
+ if(T > Thottest)
+ {
+ ihottest = i;
+ Thottest = T;
+ }
+ }
+ }
+ if(ihottest!=-1 && Thottest_p != NULL)
+ {
+ *Thottest_p = Thottest;
+ }
+ return ihottest;
+}
+
+#endif //DISCS
diff --git a/src/disc/disc_logging.c b/src/disc/disc_logging.c
index 005b26aa18f970a67f4ce57e923ad71605d84c0e..f6a3324fa84aec49c2c3662e8590cc055b029a24 100644
--- a/src/disc/disc_logging.c
+++ b/src/disc/disc_logging.c
@@ -21,9 +21,6 @@ static void close_disc_logs(struct stardata_t * stardata);
* binary struct, thus saving us the effort.
*/
-#define Mdot(X) (Max(0.0,-disc->loss[(X)].mdot / M_SUN * YEAR_LENGTH_IN_SECONDS))
-#define Jdot(X) (Max(0.0,-disc->loss[(X)].jdot))
-
void disc_logging(struct stardata_t * const stardata,
const struct binary_system_t * const binary,
@@ -400,7 +397,7 @@ void disc_logging(struct stardata_t * const stardata,
Disc_log_filepointer2d,
"%g %d %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g\n",
disc->lifetime/YEAR_LENGTH_IN_SECONDS,
- DISC_ZONE_TYPE(R,disc),
+ Disc_zone_type(R,disc),
R/R_SUN,
T,
Sigma,
@@ -651,20 +648,20 @@ void disc_logging(struct stardata_t * const stardata,
const double RMhalf = disc_half_mass_radius(disc,b);
const double RJhalf = disc_half_angular_momentum_radius(disc,b);
const double Mdotwind =
- Mdot(DISC_LOSS_GLOBAL)+
- Mdot(DISC_LOSS_ISM)+
- Mdot(DISC_LOSS_FUV)+
- Mdot(DISC_LOSS_XRAY)+
- Mdot(DISC_LOSS_INNER_EDGE)+
- Mdot(DISC_LOSS_OUTER_EDGE);
+ Disc_Mdot(DISC_LOSS_GLOBAL)+
+ Disc_Mdot(DISC_LOSS_ISM)+
+ Disc_Mdot(DISC_LOSS_FUV)+
+ Disc_Mdot(DISC_LOSS_XRAY)+
+ Disc_Mdot(DISC_LOSS_INNER_EDGE)+
+ Disc_Mdot(DISC_LOSS_OUTER_EDGE);
const double Jdotwind =
- Jdot(DISC_LOSS_GLOBAL)+
- Jdot(DISC_LOSS_ISM)+
- Jdot(DISC_LOSS_FUV)+
- Jdot(DISC_LOSS_XRAY)+
- Jdot(DISC_LOSS_INNER_EDGE)+
- Jdot(DISC_LOSS_OUTER_EDGE);
+ Disc_Jdot(DISC_LOSS_GLOBAL)+
+ Disc_Jdot(DISC_LOSS_ISM)+
+ Disc_Jdot(DISC_LOSS_FUV)+
+ Disc_Jdot(DISC_LOSS_XRAY)+
+ Disc_Jdot(DISC_LOSS_INNER_EDGE)+
+ Disc_Jdot(DISC_LOSS_OUTER_EDGE);
const double hwind = Is_really_not_zero(Mdotwind) ?
(Jdotwind / Mdotwind) : 0.0;
@@ -672,10 +669,10 @@ void disc_logging(struct stardata_t * const stardata,
const double min_Toomre_Q = disc_minimum_Toomre_Q(disc,b);
- double RZ0,RZ1; /* zone crossing radii */
- //disc_zone_crossing_radii(disc,binary,&RZ0,&RZ1);
- RZ0 = RZ1 = 0.0;
-
+ /* zone crossing radii */
+ double * RZcross;
+ disc_zone_crossing_radii(disc,b,&RZcross);
+
/* output */
binary_c_filter_fprintf(
Disc_log_filepointer,
@@ -689,10 +686,10 @@ void disc_logging(struct stardata_t * const stardata,
// thermal zone information
disc->n_thermal_zones,
- (disc->thermal_zones[0].valid==TRUE ? disc->thermal_zones[0].type : DISC_ZONE_TYPE_UNDEFINED),
- (disc->thermal_zones[1].valid==TRUE ? disc->thermal_zones[1].type : DISC_ZONE_TYPE_UNDEFINED),
- (disc->thermal_zones[2].valid==TRUE ? disc->thermal_zones[2].type : DISC_ZONE_TYPE_UNDEFINED),
-
+ Disc_valid_zone_type(disc,0),
+ Disc_valid_zone_type(disc,1),
+ Disc_valid_zone_type(disc,2),
+
// radii
disc->Rin/R_SUN,
disc->Rout/R_SUN,
@@ -712,9 +709,9 @@ void disc_logging(struct stardata_t * const stardata,
stardata->star[0].roche_radius,
stardata->star[1].roche_radius,
- // zone crossing radii
- RZ0/R_SUN,
- RZ1/R_SUN,
+ // zone crossing radii, or negative
+ RZcross != NULL ? (RZcross[0]/R_SUN) : -1.0,
+ RZcross != NULL ? (RZcross[1]/R_SUN) : -1.0,
// misc.
disc->torqueF,
@@ -735,25 +732,25 @@ void disc_logging(struct stardata_t * const stardata,
// disc mass loss rates
- Mdot(DISC_LOSS_INNER_VISCOUS),
- Mdot(DISC_LOSS_GLOBAL),
- Mdot(DISC_LOSS_INNER_L2_CROSSING),
- Mdot(DISC_LOSS_ISM),
- Mdot(DISC_LOSS_FUV),
- Mdot(DISC_LOSS_XRAY),
- Mdot(DISC_LOSS_INNER_EDGE),
- Mdot(DISC_LOSS_OUTER_EDGE),
+ Disc_Mdot(DISC_LOSS_INNER_VISCOUS),
+ Disc_Mdot(DISC_LOSS_GLOBAL),
+ Disc_Mdot(DISC_LOSS_INNER_L2_CROSSING),
+ Disc_Mdot(DISC_LOSS_ISM),
+ Disc_Mdot(DISC_LOSS_FUV),
+ Disc_Mdot(DISC_LOSS_XRAY),
+ Disc_Mdot(DISC_LOSS_INNER_EDGE),
+ Disc_Mdot(DISC_LOSS_OUTER_EDGE),
// angular momentum rates
- Jdot(DISC_LOSS_INNER_VISCOUS),
- Jdot(DISC_LOSS_GLOBAL),
- Jdot(DISC_LOSS_INNER_L2_CROSSING),
- Jdot(DISC_LOSS_ISM),
- Jdot(DISC_LOSS_FUV),
- Jdot(DISC_LOSS_XRAY),
- Jdot(DISC_LOSS_INNER_EDGE),
- Jdot(DISC_LOSS_OUTER_EDGE),
- Jdot(DISC_LOSS_BINARY_TORQUE),
+ Disc_Jdot(DISC_LOSS_INNER_VISCOUS),
+ Disc_Jdot(DISC_LOSS_GLOBAL),
+ Disc_Jdot(DISC_LOSS_INNER_L2_CROSSING),
+ Disc_Jdot(DISC_LOSS_ISM),
+ Disc_Jdot(DISC_LOSS_FUV),
+ Disc_Jdot(DISC_LOSS_XRAY),
+ Disc_Jdot(DISC_LOSS_INNER_EDGE),
+ Disc_Jdot(DISC_LOSS_OUTER_EDGE),
+ Disc_Jdot(DISC_LOSS_BINARY_TORQUE),
// specific angular momentum loss rate in winds
hwind,
@@ -828,6 +825,7 @@ void disc_logging(struct stardata_t * const stardata,
stardata->star[1].derivative[DERIVATIVE_STELLAR_CORE_MASS]
);
binary_c_fflush(Disc_log_filepointer);
+ Safe_free(RZcross);
}
#endif // DISC_LOG
diff --git a/src/disc/disc_new_zone_list.c b/src/disc/disc_new_zone_list.c
index 129b083dd3322b757e4c86621e2620b3bc0fb4fe..a19c82e3c202836cb8bd75f54b06ba6e2070cdfa 100644
--- a/src/disc/disc_new_zone_list.c
+++ b/src/disc/disc_new_zone_list.c
@@ -59,9 +59,12 @@ Boolean disc_new_zone_list(struct disc_t * const disc,
* Create pointers to the zones
* (for convenience and ease of reading in English)
*/
- struct disc_thermal_zone_t * radiative_in = &(thermal_zones[DISC_ZONE_TYPE_RADIATIVE_INNER]);
- struct disc_thermal_zone_t * radiative_out = &(thermal_zones[DISC_ZONE_TYPE_RADIATIVE_OUTER]);
- struct disc_thermal_zone_t * viscous = &(thermal_zones[DISC_ZONE_TYPE_VISCOUS]);
+ struct disc_thermal_zone_t * const radiative_in =
+ &(thermal_zones[DISC_ZONE_TYPE_RADIATIVE_INNER]);
+ struct disc_thermal_zone_t * const radiative_out =
+ &(thermal_zones[DISC_ZONE_TYPE_RADIATIVE_OUTER]);
+ struct disc_thermal_zone_t * const viscous =
+ &(thermal_zones[DISC_ZONE_TYPE_VISCOUS]);
/*
* Set up the outer radiative zone, in which Rref should lie
@@ -196,7 +199,7 @@ Boolean disc_new_zone_list(struct disc_t * const disc,
* Default all thermal zones to be invalid
*/
{
- int i;
+ Disc_zone_counter i;
for(i=0;i<DISCS_MAX_N_ZONES;i++)
{
thermal_zones[i].valid = FALSE;
diff --git a/src/disc/disc_new_zone_radii.c b/src/disc/disc_new_zone_radii.c
index 32db3f7b297e6917fd7170d81ece0f63668b2d0c..991fdeab27865c6fcaef2b26174a7d5bde186808 100644
--- a/src/disc/disc_new_zone_radii.c
+++ b/src/disc/disc_new_zone_radii.c
@@ -93,7 +93,7 @@ Boolean disc_new_zone_radii(struct disc_t * const disc)
* Debugging:
* Show the crossing radius matrix
*/
- int i,j;
+ Disc_zone_counter i,j;
for(i=0;i<DISCS_MAX_N_ZONES;i++)
{
for(j=0;j<=i;j++)
@@ -119,34 +119,10 @@ 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.
- *
* Zones that have already been used are ignored.
*/
- double Thottest = 0.0;
- Disc_zone_counter ihottest = -1;
- for(i=0;i<DISCS_MAX_N_ZONES;i++)
- {
- if(P[i].valid == TRUE)
- {
- const double T = power_law(&P[i].Tlaw,r);
- if(T > Thottest)
- {
- ihottest = i;
- Thottest = T;
- }
- 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);
- }
- }
-
+ Disc_zone_counter ihottest = disc_hottest_zone_at_radius(P,r,FALSE,NULL);
+
Zprintf("Hottest at r = %30.15e %s is zone %u\n",
r/R_SUN,
R_SOLAR,
@@ -290,8 +266,8 @@ Boolean disc_new_zone_radii(struct disc_t * const disc)
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,
- DISC_ZONE_TYPE_STRING(z),
- DISC_ZONE_VALID_STRING(z),
+ Disc_zone_type_string(z),
+ Disc_zone_valid_string(z),
z->rstart/R_SUN,
R_SOLAR,
z->rend/R_SUN,
diff --git a/src/disc/disc_parameters.h b/src/disc/disc_parameters.h
index 60c78dc59ece28afab8eaae04231e2c041a93bfb..deb608bd40c6626a79e8f3a087b792ea845d7d01 100644
--- a/src/disc/disc_parameters.h
+++ b/src/disc/disc_parameters.h
@@ -184,7 +184,7 @@
* and these work fine). Defaults: 1e-6 and 1000.
*/
-#define DISC_TOLERANCE 1e-4
+#define DISC_TOLERANCE 1e-6
#define DISC_BISECTION_MAX_ITERATIONS 1000
/*
diff --git a/src/disc/disc_prototypes.h b/src/disc/disc_prototypes.h
index 15193660c11ed9bb663f00fd16be2d0ebf1afb9e..1c7b06c81d3cd9cc52f787d9a2f0249a2383a9e5 100644
--- a/src/disc/disc_prototypes.h
+++ b/src/disc/disc_prototypes.h
@@ -479,11 +479,10 @@ void disc_equation_checker(const double r,
int disc_n_valid_zones(const struct disc_t * const disc);
-void disc_zone_crossing_radii(struct disc_t * const disc,
- const struct binary_system_t * const binary,
- double * const r1,
- double * const r2);
+void disc_zone_crossing_radii(struct disc_t * const disc,
+ const struct binary_system_t * const binary,
+ double ** radii_p);
Boolean disc_new_zone_radii(struct disc_t * const disc);
@@ -494,6 +493,12 @@ Boolean disc_monotonic_check(struct stardata_t * const stardata,
const int n
);
+
+Disc_zone_counter disc_hottest_zone_at_radius(struct disc_thermal_zone_t * zones,
+ const double radius,
+ const Boolean all,
+ double * Thottest_p);
+
#endif // DISCS
#endif // DISC_PROTOTYPES_H
diff --git a/src/disc/disc_rezone.c b/src/disc/disc_rezone.c
index 4de093e374efb5dfadbd6dfa683f0a26f24132c7..8f56981de042381db6929e2663c36ad689e5d65b 100644
--- a/src/disc/disc_rezone.c
+++ b/src/disc/disc_rezone.c
@@ -116,7 +116,7 @@ Boolean disc_rezone(struct disc_t * const disc,
* to be correct at the new rstart
*/
Rezone_debug("set T power laws\n");
- zone = DISC_INNER_ZONE(disc);
+ zone = Disc_inner_zone(disc);
zone->Tlaw.A0 = power_law(&zone->Tlaw,
zone->rstart);
zone->Tlaw.R0 = zone->rstart;
diff --git a/src/disc/disc_show_thermal_zones.c b/src/disc/disc_show_thermal_zones.c
index d94ce8300bb74b3370cedf76d62360e992a724df..3c023d32670fdebbaa4462b2d38896c3e04f4895 100644
--- a/src/disc/disc_show_thermal_zones.c
+++ b/src/disc/disc_show_thermal_zones.c
@@ -39,8 +39,8 @@ void disc_show_thermal_zones(struct disc_t * disc,
printf("Zone %d at %p :",i,z);
fflush(stdout);
printf("type %s : %s %s: ",
- DISC_ZONE_TYPE_STRING(z),
- DISC_ZONE_VALID_STRING(z),
+ Disc_zone_type_string(z),
+ Disc_zone_valid_string(z),
NORMAL_COLOUR(z)
);
diff --git a/src/disc/disc_structure_functions.c b/src/disc/disc_structure_functions.c
index 5a088001e8b2cc20b8682d24f8280f38273e7038..003af166ac58604c7d160e62244111ae7d811d5a 100644
--- a/src/disc/disc_structure_functions.c
+++ b/src/disc/disc_structure_functions.c
@@ -35,7 +35,7 @@ double Pure_function disc_inner_edge_temperature(const struct disc_t * const dis
/*
* Temperature at the disc's inner edge (K)
*/
- const struct disc_thermal_zone_t * zone = DISC_INNER_ZONE(disc);
+ const struct disc_thermal_zone_t * zone = Disc_inner_zone(disc);
return Max(0.0,power_law(&(zone->Tlaw),zone->rstart));
}
@@ -44,7 +44,7 @@ double Pure_function disc_outer_edge_temperature(const struct disc_t * const dis
/*
* Temperature at the disc's outer edge (K)
*/
- const struct disc_thermal_zone_t * zone = DISC_OUTER_ZONE(disc);
+ const struct disc_thermal_zone_t * zone = Disc_outer_zone(disc);
return Max(0.0,power_law(&(zone->Tlaw),zone->rend));
}
@@ -66,7 +66,7 @@ double generic_power_law(double radius,
* if radius is negative, we seek only the nearest zone
*/
radius = - radius;
- z = DISC_NEAREST_ZONE(radius,disc);
+ z = Disc_nearest_zone(radius,disc);
}
else
{
@@ -74,7 +74,7 @@ double generic_power_law(double radius,
* But by default, we seek a zone INSIDE the disc,
* i.e. there is some error checking.
*/
- z = DISC_ZONE(radius,disc);
+ z = Disc_zone(radius,disc);
}
return power_law(n==POWER_LAW_TEMPERATURE ?
diff --git a/src/disc/disc_thermal_zones.h b/src/disc/disc_thermal_zones.h
index 26d30661de174d9e73fbfe86fdea6beb7f97cd7c..d6a4207a0e96c4f83f5c5187b8867d8f19bd29a3 100644
--- a/src/disc/disc_thermal_zones.h
+++ b/src/disc/disc_thermal_zones.h
@@ -15,44 +15,44 @@
/*
* Given a radius, return the zone number
*/
-//#define DISC_ZONE_N(R,D) ((D)->n_thermal_zones==1 ? 0 : (R) < (D)->thermal_zones[1].rstart ? 0 : 1)
+//#define Disc_zone_n(R,D) ((D)->n_thermal_zones==1 ? 0 : (R) < (D)->thermal_zones[1].rstart ? 0 : 1)
// was 2.38s
-#define DISC_ZONE_N(R,D) (disc_zone_n((R),(D)))
+#define Disc_zone_n(R,D) (disc_zone_n((R),(D)))
-#define DISC_NEAREST_ZONE_N(R,D) (disc_nearest_zone_n((R),(D)))
+#define Disc_nearest_zone_n(R,D) (disc_nearest_zone_n((R),(D)))
/*
* Disc's inner and outer zone numbers
*/
-#define DISC_INNER_ZONE_N(D) (0)
-#define DISC_OUTER_ZONE_N(D) ((D)->n_thermal_zones>0 ? ((D)->n_thermal_zones-1) : 0)
+#define Disc_inner_zone_n(D) (0)
+#define Disc_outer_zone_n(D) ((D)->n_thermal_zones>0 ? ((D)->n_thermal_zones-1) : 0)
/*
* Given a radius, return the zone type
*/
-#define DISC_ZONE_TYPE(R,D) ((D)->thermal_zones[DISC_ZONE_N((R),(D))].type)
+#define Disc_zone_type(R,D) ((D)->thermal_zones[Disc_zone_n((R),(D))].type)
/*
* Given a radius, return (a pointer to) the current zone structure
*/
-#define DISC_ZONE(R,D) (&((D)->thermal_zones[DISC_ZONE_N((R),(D))]))
-#define DISC_NEAREST_ZONE(R,D) (&((D)->thermal_zones[DISC_NEAREST_ZONE_N((R),(D))]))
+#define Disc_zone(R,D) (&((D)->thermal_zones[Disc_zone_n((R),(D))]))
+#define Disc_nearest_zone(R,D) (&((D)->thermal_zones[Disc_nearest_zone_n((R),(D))]))
/*
* Pointers to the inner and outer zone, respectively
*/
-#define DISC_INNER_ZONE(D) (&((D)->thermal_zones[DISC_INNER_ZONE_N(D)]))
-#define DISC_OUTER_ZONE(D) (&((D)->thermal_zones[DISC_OUTER_ZONE_N(D)]))
+#define Disc_inner_zone(D) (&((D)->thermal_zones[Disc_inner_zone_n(D)]))
+#define Disc_outer_zone(D) (&((D)->thermal_zones[Disc_outer_zone_n(D)]))
/*
* Pointer to the previous and next zone,
* or NULL if these do not exist.
*/
-#define DISC_PREVIOUS_ZONE(D,Z) ( \
+#define Disc_previous_zone(D,Z) ( \
(Z)->type==0 ? NULL \
: &((D)->thermal_zones[(Z)->type-1]))
-#define DISC_NEXT_ZONE(D,Z) ( \
+#define Disc_next_zone(D,Z) ( \
(Z)->type==DISCS_MAX_N_ZONES ? NULL \
: &((D)->thermal_zones[(Z)->type+1]))
@@ -66,7 +66,19 @@
#define DISC_ZONE_TYPE_RADIATIVE_OUTER 2
#define DISC_ZONE_TYPE_UNDEFINED -1
-#define DISC_ZONE_TYPE_STRING(Z) \
+#define Disc_valid_zone_type_z(Z) \
+ ( \
+ ((Z)->valid == TRUE) ? \
+ ((Z)->type) : \
+ DISC_ZONE_TYPE_UNDEFINED \
+ )
+
+#define Disc_valid_zone_type(DISC,N) \
+ ( \
+ Disc_valid_zone_type_z(&(DISC)->thermal_zones[(N)]) \
+ )
+
+#define Disc_zone_type_string(Z) \
( \
(Z->type)==DISC_ZONE_TYPE_VISCOUS ? "Viscous" : \
(Z->type)==DISC_ZONE_TYPE_RADIATIVE_INNER ? "Rad in " : \
@@ -74,7 +86,7 @@
"undefined" \
)
-#define DISC_ZONE_VALID_STRING(Z) \
+#define Disc_zone_valid_string(Z) \
((Z)->valid ? " valid" : "invalid")
/*
diff --git a/src/disc/disc_trisector.c b/src/disc/disc_trisector.c
index 77dc2f099508d56304473c02138b23bc6a6eef95..dde3eb6e857c48970361b6f78cd40f5c8076a180 100644
--- a/src/disc/disc_trisector.c
+++ b/src/disc/disc_trisector.c
@@ -115,8 +115,8 @@ void disc_trisector2(struct stardata_t * stardata,
M/M_SUN,disc->M/M_SUN,
J,disc->J,
F,disc->F,
- DISC_INNER_ZONE(disc)->rstart,
- DISC_INNER_ZONE(disc)->rend
+ Disc_inner_zone(disc)->rstart,
+ Disc_inner_zone(disc)->rend
);
residual[0] = M/disc->M - 1.0;
residual[1] = J/disc->J - 1.0;
diff --git a/src/disc/disc_zone_crossing_radii.c b/src/disc/disc_zone_crossing_radii.c
index 6fa19dd299437ef512a792a4168541ce3c021fe4..5bd666cfc985c483e197756e5f97b31bc8529cac 100644
--- a/src/disc/disc_zone_crossing_radii.c
+++ b/src/disc/disc_zone_crossing_radii.c
@@ -2,21 +2,143 @@
#ifdef DISCS
-void disc_zone_crossing_radii(struct disc_t * const disc,
- const struct binary_system_t * const binary,
- double * const r1,
- double * const r2)
+void disc_zone_crossing_radii(
+ struct disc_t * const disc,
+ const struct binary_system_t * const binary,
+ double ** radii_p
+ )
{
/*
- * Set r1,r2 to the zone crossing radii
- */
- struct disc_thermal_zone_t W[DISCS_MAX_N_ZONES+1];
- disc_new_zone_list(disc,binary,W);
+ * Set the crossing radii of disc in **radii.
+ *
+ * If the new zone list cannot be made, set *radii to NULL.
+ */
+ struct disc_thermal_zone_t * P =
+ Malloc(sizeof(struct disc_thermal_zone_t)*
+ (DISCS_MAX_N_ZONES));
- *r1 = disc_Rcross(W+0,W+1); // VR
- *r2 = disc_Rcross(W+1,W+2); // IO
+ if(disc_new_zone_list(disc,binary,P) == DISC_ZONES_OK)
+ {
+ Disc_zone_counter i,j,current=0;
- printf("r1 = %g, r2 = %g\n",*r1,*r2);
+ /*
+ * Allocate space for radii list: this should
+ * be freed by the calling function
+ */
+ *radii_p = Malloc(sizeof(double)*DISCS_MAX_N_ZONES);
+ double * const radii = *radii_p;
+
+ /*
+ * Default all radii to unphysical -1.0
+ * and make all zones valid to start with
+ */
+ for(i=0;i<DISCS_MAX_N_ZONES;i++)
+ {
+ radii[i] = -1.0;
+ P[i].valid = TRUE;
+ }
+
+ /*
+ * The first radius should be the shortest
+ */
+ double r = DISC_LARGE_RADIUS;
+ for(i=0;i<DISCS_MAX_N_ZONES;i++)
+ {
+ for(j=0;j<DISCS_MAX_N_ZONES;j++)
+ {
+ const double Rx = disc_Rcross(&P[i],&P[j]);
+ if(Rx < r)
+ {
+ r = Rx;
+ current = i;
+ }
+ }
+ }
+ P[current].valid = FALSE;
+ //printf("Zone at 1cm %u\n",current);
+
+ /*
+ * Loop to find the next
+ */
+ Boolean valid = TRUE;
+ Disc_zone_counter icount = 0;
+ while(valid == TRUE)
+ {
+ double next_crossing_radius = DISC_LARGE_RADIUS;
+ Disc_zone_counter next_crossing_index = current;
+ Boolean found = FALSE;
+
+ /*
+ * Find the next crossing radius
+ */
+ for(i=0;i<DISCS_MAX_N_ZONES;i++)
+ {
+ /*
+ printf("check %u : current? %s : valid? %s\n",
+ i,
+ Yesno(i==current),
+ Yesno(P[i].valid));
+ */
+ if(i != current &&
+ P[i].valid == TRUE)
+ {
+ const double Rx = disc_Rcross(P + current,
+ P + i);
+ if(Rx > r &&
+ Rx < next_crossing_radius)
+ {
+ next_crossing_radius = Rx;
+ next_crossing_index = i;
+ found = TRUE;
+ }
+ }
+ }
+ /*
+ printf("found? %s : next X at Rx = %g with zone i = %u, icount = %u\n",
+ Yesno(found),
+ next_crossing_radius,
+ next_crossing_index,
+ icount);
+ */
+ if(found == TRUE)
+ {
+ /*
+ * Save to radii list
+ */
+ radii[icount] = next_crossing_radius;
+ //printf("set radius %u to %g\n",icount,radii[icount]);
+ icount++;
+
+ /*
+ * Make current zone invalid
+ */
+ current = next_crossing_index;
+ P[current].valid = FALSE;
+
+ /*
+ * Stop looping if all zones are invalid
+ */
+ valid = FALSE;
+ for(i=0;i<DISCS_MAX_N_ZONES;i++)
+ {
+ if(P[i].valid == TRUE)
+ {
+ valid = TRUE;
+ }
+ }
+ }
+ else
+ {
+ valid = FALSE;
+ }
+ }
+ }
+ else
+ {
+ *radii_p = NULL;
+ }
+
+ Safe_free(P);
}
#endif
diff --git a/src/disc/disc_zone_n.c b/src/disc/disc_zone_n.c
index 881028f22208ff099f10d5d89264540b220fa9b1..8bf29f19d64632c49cb18ecc67fe18b9faa95b28 100644
--- a/src/disc/disc_zone_n.c
+++ b/src/disc/disc_zone_n.c
@@ -15,18 +15,18 @@ Disc_zone_counter disc_zone_n(const double radius,
*/
#if DISC_DEBUG >= 3
- double d = radius-DISC_INNER_ZONE(disc)->rstart;
+ double d = radius-Disc_inner_zone(disc)->rstart;
Discdebug(3,
"R = %g, Disc Rin = %g, Rout = %g, Zone R inner edge = %g, R outer edge = %g, diff = %g, diff / R = %g, diff / INNER = %g, Abs_diff %g\n",
radius,
disc->Rin,
disc->Rout,
- DISC_INNER_ZONE(disc)->rstart,
- DISC_OUTER_ZONE(disc)->rend,
+ Disc_inner_zone(disc)->rstart,
+ Disc_outer_zone(disc)->rend,
d,
d/radius,
- d/DISC_INNER_ZONE(disc)->rstart,
- Abs_diff(radius,DISC_INNER_ZONE(disc)->rstart)
+ d/Disc_inner_zone(disc)->rstart,
+ Abs_diff(radius,Disc_inner_zone(disc)->rstart)
);
#endif //DISC_DEBUG>=3
@@ -54,17 +54,17 @@ Disc_zone_counter disc_zone_n(const double radius,
}
/* compare to inner and outer edges */
- else if(unlikely(Float_same_within_eps(radius/DISC_OUTER_ZONE(disc)->rend,
+ else if(unlikely(Float_same_within_eps(radius/Disc_outer_zone(disc)->rend,
1.0,
DISC_EDGE_EPS)))
{
- return DISC_OUTER_ZONE_N(disc);
+ return Disc_outer_zone_n(disc);
}
- else if(unlikely(Float_same_within_eps(radius/DISC_INNER_ZONE(disc)->rstart,
+ else if(unlikely(Float_same_within_eps(radius/Disc_inner_zone(disc)->rstart,
1.0,
DISC_EDGE_EPS)))
{
- return DISC_INNER_ZONE_N(disc);
+ return Disc_inner_zone_n(disc);
}
else
{
@@ -153,20 +153,20 @@ Disc_zone_counter disc_zone_n(const double radius,
fprintf(stderr,
"eps in cf %g %g = %g\n",
radius,
- DISC_INNER_ZONE(disc)->rstart,
- Abs_diff(radius,DISC_INNER_ZONE(disc)->rstart));
+ Disc_inner_zone(disc)->rstart,
+ Abs_diff(radius,Disc_inner_zone(disc)->rstart));
fprintf(stderr,
"eps out cf %g %g = %g\n",
radius,
- DISC_OUTER_ZONE(disc)->rend,
- Abs_diff(radius,DISC_OUTER_ZONE(disc)->rend));
+ Disc_outer_zone(disc)->rend,
+ Abs_diff(radius,Disc_outer_zone(disc)->rend));
// failure
Backtrace;
#if DISC_DEBUG<3
// not defined above...
- double d = radius-DISC_INNER_ZONE(disc)->rstart;
+ double d = radius-Disc_inner_zone(disc)->rstart;
#endif
fprintf(stderr,
@@ -174,12 +174,12 @@ Disc_zone_counter disc_zone_n(const double radius,
radius,
disc->Rin,
disc->Rout,
- DISC_INNER_ZONE(disc)->rstart,
- DISC_OUTER_ZONE(disc)->rend,
+ Disc_inner_zone(disc)->rstart,
+ Disc_outer_zone(disc)->rend,
d,
d/radius,
- d/DISC_INNER_ZONE(disc)->rstart,
- Abs_diff(radius,DISC_INNER_ZONE(disc)->rstart)
+ d/Disc_inner_zone(disc)->rstart,
+ Abs_diff(radius,Disc_inner_zone(disc)->rstart)
);
Exit_binary_c_no_stardata(
@@ -188,12 +188,12 @@ Disc_zone_counter disc_zone_n(const double radius,
radius,
disc->Rin,
disc->Rout,
- Abs_diff(radius,DISC_INNER_ZONE(disc)->rstart),
- Abs_diff(radius,DISC_OUTER_ZONE(disc)->rend),
- Fequal(radius,DISC_INNER_ZONE(disc)->rstart),
- Fequal(radius,DISC_OUTER_ZONE(disc)->rend),
- DISC_INNER_ZONE(disc)->valid,
- DISC_OUTER_ZONE(disc)->valid
+ Abs_diff(radius,Disc_inner_zone(disc)->rstart),
+ Abs_diff(radius,Disc_outer_zone(disc)->rend),
+ Fequal(radius,Disc_inner_zone(disc)->rstart),
+ Fequal(radius,Disc_outer_zone(disc)->rend),
+ Disc_inner_zone(disc)->valid,
+ Disc_outer_zone(disc)->valid
);
diff --git a/src/nucsyn/nucsyn_parameters.h b/src/nucsyn/nucsyn_parameters.h
index c53527665d7ed88e43bd2fcb418a69eb38462366..0b67abc82d463c04c7c6059f572219758b8ce0dd 100644
--- a/src/nucsyn/nucsyn_parameters.h
+++ b/src/nucsyn/nucsyn_parameters.h
@@ -740,7 +740,7 @@ undances with s-process
/* define NUCSYN_GCE to modify yield output for gce.pl
* NB this just activates other defines, and turns some off.
*/
-#define NUCSYN_GCE
+//#define NUCSYN_GCE
#ifdef NUCSYN_GCE
diff --git a/tbse b/tbse
index 44485787d85003167450dacc257c84403176685a..a1a4ab17c1a7f57b97cdf812db8c15de4b1e4db6 100755
--- a/tbse
+++ b/tbse
@@ -575,7 +575,7 @@ DISC_TIMESTEP_FACTOR=1.0
# DISC_LOG_LEVEL_SUBTIMESTEP 2 = on with sub-timestep logging
# DISC_LOG_LEVEL_NORMAL_FIRST_DISC_ONLY = as 1 for first disc only
# DISC_LOG_LEVEL_SUBTIMESTEP_FIRST_DISC_ONLY = as 2 for first disc only
-DISC_LOG=DISC_LOG_LEVEL_SUBTIMESTEP_FIRST_DISC_ONLY
+DISC_LOG=DISC_LOG_LEVEL_NORMAL_FIRST_DISC_ONLY
# soft 2D disc log switch, same options as DISC_LOG
DISC_LOG2D=DISC_LOG_LEVEL_NONE