patch: inferred infall rates take into account gas migration

vice/src/multizone/ism.c

@@ -9,6 +9,7 @@
 #include "../utils.h"
 #include "../ism.h"
 #include "ism.h"
+#include "migration.h"
 
 
 /*
@@ -43,6 +44,7 @@ extern unsigned short update_zone_evolution(MULTIZONE *mz) {
 
 	unsigned int i;
 	double *mass_recycled = gas_recycled_in_zones(*mz);
+	double *migration_deltas = migration_gas_changes_by_zone(*mz);
 	for (i = 0; i < (*(*mz).mig).n_zones; i++) {
 		SINGLEZONE *sz = mz -> zones[i];
 
@@ -56,6 +58,12 @@ extern unsigned short update_zone_evolution(MULTIZONE *mz) {
 		 * assignment of the infall rate to NaN. This change ensures that the
 		 * infall rate will not be NaN by the time primordial_inflow is called
 		 * after one timestep has passed.
+		 *
+		 * Change Note: version 1.3.2
+		 *
+		 * The mass added to the ISM by migration is now taken into account when
+		 * computing the infall rate in gas and star formation mode. This small
+		 * correction was previously unaccounted for.
 		 */
 
 		switch (checksum((*(*sz).ism).mode)) {
@@ -67,7 +75,7 @@ extern unsigned short update_zone_evolution(MULTIZONE *mz) {
 				);
 				sz -> ism -> infall_rate = (
 					((*(*sz).ism).mass - (*(*sz).ism).specified[(*sz).timestep]
-						- mass_recycled[i]) / (*sz).dt +
+						- mass_recycled[i] - migration_deltas[i]) / (*sz).dt +
 					(*(*sz).ism).star_formation_rate + get_outflow_rate(*sz)
 				);
 				primordial_inflow(sz);
@@ -92,7 +100,7 @@ extern unsigned short update_zone_evolution(MULTIZONE *mz) {
 					*(*sz).ism).specified[(*sz).timestep + 1l];
 				double dMg = get_ism_mass_SFRmode(*sz, 0u) - (*(*sz).ism).mass;
 				sz -> ism -> infall_rate = (
-					(dMg - mass_recycled[i]) / (*sz).dt +
+					(dMg - mass_recycled[i] - migration_deltas[i]) / (*sz).dt +
 					(*(*sz).ism).star_formation_rate + get_outflow_rate(*sz)
 				);
 				sz -> ism -> mass += dMg;

vice/src/multizone/migration.c

@@ -213,6 +213,46 @@ extern void migrate(MULTIZONE *mz) {
 }
 
 
+/*
+ * Compute the net change in the ISM mass for all zones due solely to gas
+ * migration.
+ *
+ * Parameters
+ * ==========
+ * mz:		The multizone object for the current simulation.
+ *
+ * Returns
+ * =======
+ * deltas: a pointer to the array of doubles describing, component-wise, the
+ * net change in each zone's mass.
+ *
+ * header: migration.h
+ */
+extern double *migration_gas_changes_by_zone(MULTIZONE mz) {
+
+	double **changes = get_changes(mz, -1);
+	double *deltas = (double *) malloc ((*mz.mig).n_zones * sizeof(double));
+	unsigned int i;
+	for (i = 0u; i < (*mz.mig).n_zones; i++) {
+		/*
+		 * The changes array is indexed such that the ij'th component holds the
+		 * amount of mass that migrates from zone i into zone j. In other words,
+		 * the i'th column (i.e. the first axis of indexing) describes the
+		 * amount that zone i loses, while the i'th row (the second axis of
+		 * indexing) describes the amount that zone i gains.
+		 */
+		deltas[i] = 0;
+		unsigned int j;
+		/* changes[i][i] = 0 for all i's (see get_changes below). */
+		for (j = 0u; j < (*mz.mig).n_zones; j++) deltas[i] += changes[i][j];
+		for (j = 0u; j < (*mz.mig).n_zones; j++) deltas[i] -= changes[j][i];
+	}
+
+	return deltas;
+
+}
+
+
 /*
  * Updates a tracer particle's current zone number based on the zone_history
  * array at the next timestep.

vice/src/multizone/migration.h

@@ -21,6 +21,23 @@ extern "C" {
  */
 extern void migrate(MULTIZONE *mz);
 
+/*
+ * Compute the net change in the ISM mass for all zones due solely to gas
+ * migration.
+ *
+ * Parameters
+ * ==========
+ * mz:		The multizone object for the current simulation.
+ *
+ * Returns
+ * =======
+ * deltas: a pointer to the array of doubles describing, component-wise, the
+ * net change in each zone's mass.
+ *
+ * source: migration.c
+ */
+extern double *migration_gas_changes_by_zone(MULTIZONE mz);
+
 /*
  * Performs a sanity check on a given migration matrix by making sure the sum
  * of migration probabilities out of a given zone at all times is <= 1.