Merge pull request #81 from sirus20x6/feat/cold-giant-metallicity-gate

feat(accrete): metallicity-dependent giant-formation gate (fix cold-giant over-production)

Author: sirus20x6

accrete.cpp

@@ -885,6 +885,31 @@ auto accrete::dist_planetary_masses(sun &the_sun, long double inner_dust,
 
   // std::cout << planet_inner_bound << " " << planet_outer_bound << std::endl;
 
+  // --- Metallicity-dependent giant-formation gate (default-on; see const.h) -----
+  // Decide ONCE whether this star can form gas giants, from a per-star [Fe/H] draw
+  // and the observed giant-frequency-metallicity relation. Three fixed draws in a
+  // fixed order (Box-Muller [Fe/H] = 2 draws, then the acceptance draw) -> no
+  // rejection loops, serial == parallel byte-identical. If not giant-capable, gas
+  // runaway is suppressed below (crit_mass made unreachable) so bodies stay cores.
+  bool giant_capable = true;
+  {
+    const long double u1 = random_ctx->randDouble(1.0e-12L, 1.0L);  // avoid log(0)
+    const long double u2 = random_ctx->randDouble(0.0L, 1.0L);
+    const long double feh =
+        GIANT_METALLICITY_MEAN +
+        GIANT_METALLICITY_SIGMA * sqrt(-2.0L * log(u1)) * cos(2.0L * PI * u2);
+    long double p_giant = GIANT_FORMATION_NORM * stell_mass_ratio *
+                          std::pow(10.0L, GIANT_FORMATION_FEH_SLOPE * feh);
+    if (stell_mass_ratio > GIANT_FORMATION_MASS_TURNOVER_MSUN) {  // steep drop for A-stars
+      p_giant *= GIANT_FORMATION_MASS_TURNOVER_MSUN / stell_mass_ratio;
+    }
+    if (p_giant > GIANT_FORMATION_PROB_CAP) {
+      p_giant = GIANT_FORMATION_PROB_CAP;
+    }
+    const long double u_giant = random_ctx->randDouble(0.0L, 1.0L);
+    giant_capable = (u_giant < p_giant);
+  }
+
   // while there's still dust left...
   while (dust_left) {
     if (seeds != nullptr) {
@@ -952,6 +977,9 @@ auto accrete::dist_planetary_masses(sun &the_sun, long double inner_dust,
       dust_density = dust_density_coeff * sqrt(stell_mass_ratio) *
                      exp(-ALPHA * std::pow(a, 1.0 / NDENSITY));
       crit_mass = critical_limit(a, e, stell_luminosity_ratio);
+      if (!giant_capable) {
+        crit_mass = INCREDIBLY_LARGE_NUMBER;  // metallicity gate: suppress gas runaway -> rocky/ice core
+      }
       if (total_mass == PROTOPLANET_MASS && is_seed == false) {
         // std::cout << "test1\n";
         // std::cout << total_mass << " " << dust_mass << " " << gas_mass << " " << a
@@ -1009,6 +1037,9 @@ auto accrete::dist_planetary_masses(sun &the_sun, long double inner_dust,
       r_inner = inner_effect_limit(a, e, reduced_mass);
       r_outer = outer_effect_limit(a, e, reduced_mass);
       crit_mass = critical_limit(a, e, stell_luminosity_ratio);
+      if (!giant_capable) {
+        crit_mass = INCREDIBLY_LARGE_NUMBER;  // metallicity gate: suppress gas runaway -> rocky/ice core
+      }
       if (planet_inner_bound > a) {
         planet_inner_bound = a;
       }

const.h

@@ -201,12 +201,29 @@ constexpr double MIGRATION_MIN_MASS_MJUP = 0.1;     // only bodies > 0.1 M_Jup m
 constexpr double MIGRATION_GAP_MASS_MJUP = 0.5;     // Type II gap-opening threshold (per M_star)
 constexpr double MIGRATION_TYPE1_NORM_YR = 1.1E5;   // Tanaka 2002 Type I normalization (yr)
 constexpr double MIGRATION_TYPE2_NORM_YR = 7.0E5;   // Type II viscous normalization at 1 AU (yr)
-constexpr double MIGRATION_EFFICIENCY    = 0.0001;     // overall rate multiplier (primary calibration knob)
+constexpr double MIGRATION_EFFICIENCY    = 0.0005;     // overall rate multiplier (primary calibration knob)
 constexpr double MIGRATION_INNER_EDGE_AU = 0.05;    // inner trap a_stop (disk edge / cavity)
 constexpr double MIGRATION_CIRC_AU       = 0.05;    // tidal circularization locus (hot -> e=0)
 constexpr double GIANT_ECC_SIGMA         = 0.215;   // Rayleigh sigma -> <e>~0.27 for cold giants
 constexpr double GIANT_ECC_MAX           = 0.80;    // eccentricity cap
 
+/* Metallicity-dependent giant-formation gate (default-on). StarGen's accretion
+ * forms a gas giant in ~84% of systems (crit_mass is low everywhere); the observed
+ * giant frequency is ~10-20% and rises steeply with stellar [Fe/H]. We draw a
+ * per-star [Fe/H] ~ N(mean,sigma) -- solar-neighborhood FGK metallicity distribution
+ * is ~Gaussian, mean ~0, sigma ~0.20 dex (Casagrande 2011; Nordstrom 2004) -- and
+ * accept the star as giant-capable with probability
+ *   P = NORM * (Mstar/Msun) * 10^(SLOPE*[Fe/H])   (Fischer & Valenti 2005 slope ~2.0),
+ * NORM calibrated so the population integrates to ~10-20% (closed form
+ * <P> = NORM*10^(SLOPE*mean)*exp((SLOPE*ln10)^2*sigma^2/2)). A non-giant-capable
+ * star has gas runaway suppressed (crit_mass made unreachable) -> rocky/ice cores. */
+constexpr double GIANT_METALLICITY_MEAN  = 0.0;    // dex (solar-centered FGK MDF)
+constexpr double GIANT_METALLICITY_SIGMA = 0.20;   // dex (Casagrande 2011 / Nordstrom 2004)
+constexpr double GIANT_FORMATION_NORM    = 0.10;   // calibration knob (-> ~12% integrated)
+constexpr double GIANT_FORMATION_FEH_SLOPE = 2.0;  // Fischer & Valenti 2005
+constexpr double GIANT_FORMATION_PROB_CAP = 0.90;  // cap (relation validated to [Fe/H]~+0.5)
+constexpr double GIANT_FORMATION_MASS_TURNOVER_MSUN = 2.0; // occurrence drops above ~2 Msun (Reffert 2015)
+
 // SI-unit physical constants used by the enviro.cpp acceleration helpers. These
 // are the EXACT values those functions previously redefined locally, kept
 // verbatim so generated output stays byte-identical. They intentionally differ