Use a habitability factor when deciding on planet populations

data/systems/00_sol.lua

@@ -16,7 +16,7 @@ local mercury = CustomSystemBody:new('Mercury', 'PLANET_TERRESTRIAL')
 	:seed(6)
 	:radius(f(38,100))
 	:mass(f(55,1000))
-	:temp(340)
+	:temp(375)
 	:semi_major_axis(f(387,1000))
 	:eccentricity(f(205,1000))
 	:inclination(math.deg2rad(7.0))
@@ -66,7 +66,7 @@ local earth = CustomSystemBody:new('Earth', 'PLANET_TERRESTRIAL')
 	:atmos_density(f(1,1))
 	:atmos_oxidizing(f(99,100))
 	:ocean_cover(f(7,10))
-	:ice_cover(f(5,10))
+	:ice_cover(f(3,100))
 	:life(f(9,10))
 	:orbital_phase_at_start(fixed.deg2rad(f(336,1)))
 	:rings(false)
@@ -117,12 +117,14 @@ local moon = {
 		:semi_major_axis(f(257,100000))
 		:eccentricity(f(549,10000))
 		:height_map('moon.hmap',1)
+		:metallicity(f(1,2))
 		:inclination(math.deg2rad(5.145))
 		:rotation_period(f(273,10))
 		:axial_tilt(fixed.deg2rad(f(668,100)))
 		:orbital_phase_at_start(fixed.deg2rad(f(0,1)))
 		:rotational_phase_at_start(fixed.deg2rad(f(0,1)))
 		:volcanicity(f(0,1))
+		:ice_cover(f(5,1000))
 		:rings(false),
 	{
 		CustomSystemBody:new('Tranquility Base', 'STARPORT_SURFACE')
@@ -147,7 +149,7 @@ local mars = CustomSystemBody:new('Mars', 'PLANET_TERRESTRIAL')
 	:rotation_period(f(1027,1000))
 	:axial_tilt(fixed.deg2rad(f(2519,100)))
 	-- XXX composition copied from earth until there's a way to indicate terraformed
-	:metallicity(f(14.5,5))
+	:metallicity(f(4,5))
 	:volcanicity(f(2,10))
 	:atmos_density(f(489,1000))
 	:atmos_oxidizing(f(950001,1000000))

data/systems/01_epsilon_eridani.lua

@@ -48,7 +48,7 @@ local atlantica = CustomSystemBody:new('Atlantica', 'PLANET_TERRESTRIAL')
 	:axial_tilt(fixed.deg2rad(f(26,10)))
 	:metallicity(f(5,6))
 	:volcanicity(f(6,10))
-	:atmos_density(f(9,1))
+	:atmos_density(f(2,1))
 	:atmos_oxidizing(f(1,1))
 	:ocean_cover(f(8,10))
 	:ice_cover(f(0,1))
@@ -65,7 +65,7 @@ local newhope = CustomSystemBody:new('New Hope', 'PLANET_TERRESTRIAL')
 	:rotation_period(f(4,6))
 	:axial_tilt(fixed.deg2rad(f(1741,100)))
 	:metallicity(f(5,6))
-	:volcanicity(f(68,100))
+	:volcanicity(f(28,100))
 	:atmos_density(f(15,10))
 	:atmos_oxidizing(f(7,10))
 	:ocean_cover(f(45,100))

data/systems/custom/00_sol.json

@@ -62,7 +62,7 @@
 			],
 			"atmosDensity": "f0/0",
 			"atmosOxidizing": "f0/0",
-			"averageTemp": 340,
+			"averageTemp": 375,
 			"axialTilt": "f0/749614",
 			"eccentricity": "f0/880468295",
 			"hasRings": false,
@@ -168,7 +168,7 @@
 			"seed": 1857401659,
 			"semiMajorAxis": "f1/0",
 			"type": "PLANET_TERRESTRIAL",
-			"volatileIces": "f0/2147483648",
+			"volatileIces": "f0/128849019",
 			"volatileLiquid": "f0/3006477107",
 			"volcanicity": "f0/429496729"
 		},
@@ -523,7 +523,7 @@
 			"inclination": "f0/385675994",
 			"life": "f0/0",
 			"mass": "f0/51539607",
-			"metallicity": "f0/0",
+			"metallicity": "f0/2147483648",
 			"name": "Moon",
 			"orbitalOffset": "f1/1244416494",
 			"orbitalPhase": "f0/0",
@@ -535,7 +535,7 @@
 			"seed": 4294967291,
 			"semiMajorAxis": "f0/11038065",
 			"type": "PLANET_TERRESTRIAL",
-			"volatileIces": "f0/0",
+			"volatileIces": "f0/214748365",
 			"volatileLiquid": "f0/0",
 			"volcanicity": "f0/0"
 		},
@@ -642,7 +642,7 @@
 			"inclination": "f0/138678443",
 			"life": "f0/429496729",
 			"mass": "f0/459561500",
-			"metallicity": "f2/3435973836",
+			"metallicity": "f0/3435973836",
 			"name": "Mars",
 			"orbitalOffset": "f4/2960333524",
 			"orbitalPhase": "f0/899537940",
@@ -2759,4 +2759,4 @@
 	"stars": [
 		"STAR_G"
 	]
-}
+}

src/galaxy/StarSystemGenerator.cpp

@@ -37,7 +37,8 @@ static const fixedf<48> EARTH_MASS_TO_VOL_MM3 = fixedf<48>(25946, 100); // 259.4
 // Convert a distance in megameters to earth radii
 static const fixedf<48> MM_TO_EARTH_RAD = fixedf<48>(15678, 100000);
 
-// max surface gravity for a permanent human settlement
+// min and max surface gravity for a permanent human settlement
+static const double MIN_SETTLEMENT_SURFACE_GRAVITY = 1; // m/s2 .. roughly 0.1 g - can't settle a planet if nothing sticks to the surface.
 static const double MAX_SETTLEMENT_SURFACE_GRAVITY = 50; // m/s2 .. roughly 5 g
 
 static const Uint32 POLIT_SEED = 0x1234abcd;
@@ -1454,6 +1455,82 @@ void PopulateStarSystemGenerator::PositionSettlementOnPlanet(SystemBody *sbody,
 	sbody->SetOrbitFromParameters();
 }
 
+static fixed FalloffFunction(const fixed x, const fixed mean, const fixed standard_dev)
+{
+	// Doing an actual normalized gaussian function (normal distribution curve) in
+	// fixed point is rather a pain, so instead this sort of approximates a normal
+	// distribution, if you squint.
+
+	fixed z = (x - mean) / (fixed(4, 1) * standard_dev);
+	fixed lorentz = fixed(1, 1) / (fixed(1, 1) + (z * z));
+	return lorentz * lorentz * lorentz * lorentz * lorentz * lorentz * lorentz * lorentz * lorentz;
+}
+
+static fixed CalculateHabitability(SystemBody *sbody)
+{
+	PROFILE_SCOPED()
+
+	// The habitability score is 1.0 for something exactly like Earth, and 0 for completely uninhabitable.
+
+	// Assume everything is like Earth to start with. We then calculate normal distribution-like curves so that
+	// the further conditions deviate from Earth-like, the worse the habitability gets. However, we have to be a
+	// little lenient - if we model things too closely to reality then nobody wants to live anywhere except on
+	// M-class planets, and that isn't so fun.
+	fixed habitability = fixed(1, 1);
+
+	// Surface gravity - mean at 7.5 with standard deviation of 7.5, so slightly less than Earth gravity is a
+	// pleasure to work in, zero gravity is still habitable, but has a penalty, while heavier than Earth gravity
+	// gets progressively worse. Multiply the result a small amount so that Earth itself is back up to 1.0
+	const fixed gravity = sbody->CalcSurfaceGravityAsFixed();
+	habitability *= FalloffFunction(gravity, fixed(15, 2), fixed(15, 2)) * fixed(106, 100);
+
+	// If the entire planet is covered in ice (ice coverage = 1) then that's hugely inconvenient for industry and
+	// for living, while zero ice wouldn't cause any problems at all. So set mean to 0 and have a rapid fall off towards 1.
+	const fixed ice_coverage = sbody->GetVolatileIcesAsFixed();
+	habitability *= FalloffFunction(ice_coverage, fixed(0, 1), fixed(4, 10));
+
+	// On the other hand, having a mostly oceanic planet would be great for temperature regulation, Earth-like clouds, etc.
+	// So set mean to 0.7. However, no water at all would be a huge problem. So we still want a rapid fall-off.
+	const fixed water_coverage = sbody->GetVolatileLiquidAsFixed();
+	habitability *= FalloffFunction(water_coverage, fixed(7, 10), fixed(3, 10));
+
+	// Generally humans do not want to play a live action game of "the floor is lava". Treat volcanicity like ice, but worse.
+	const fixed volcanoes = sbody->GetVolcanicityAsFixed();
+	habitability *= FalloffFunction(volcanoes, fixed(0, 1), fixed(3, 10));
+
+	// We don't treat air pressure, temperature, and atmospheric composition as independent, to avoid penalising moons
+	// too heavily. Fundamentally, if there is very low air pressure, then the temperature and composition don't matter.
+	// So first check air pressure, and if it's below a threshold then skip over the other two.
+
+	// For atmospheric pressure, we'll use the gas density, since it's already in fixed point and is good enough for this
+	// rough habitability calculation. 1.225 is Earth atmosphere density, so that's our baseline, with medium fall-off.
+	const fixed density = sbody->GetVolatileGasAsFixed();
+	habitability *= FalloffFunction(density, fixed(1225, 1000), fixed(11, 10));
+
+	if (density > fixed(6, 10)) {
+		// Temperature is in Kelvin. Average temperature on Earth's service is 15 Celsius, so 273+15 is our mean.
+		// Humans live at extremes on Earth, plus we're dealing with the future, so 50 for a standard deviation seems
+		// reasonable, though narrow in astrophysical terms
+		const fixed temperature = fixed(sbody->GetAverageTemp(), 1);
+		habitability *= FalloffFunction(temperature, fixed(288, 1), fixed(50, 1));
+
+		// Breathable air is super convenient, while air that actively kills you is even worse than a vacuum. So let's heavily
+		// favour oxidising atmospheres over reducing onces.
+		const fixed air_type = sbody->GetAtmosOxidizingAsFixed();
+		habitability *= FalloffFunction(air_type, fixed(1, 1), fixed(3, 10));
+	}
+
+	// In practice a system body that doesn't have a regular day/night cycle would presumably be roasting hot on one side and
+	// freezing on the other. Not impossible to live on, but if the boundary is slowly shifting that would be pretty annoying.
+	// On the other hand, an extremely rapid day/night cycle would be very disorientating, assuming that everything isn't
+	// flung off from the planet's surface anyway. EXCEPT - if the planet has a higher gravity, then these two things could
+	// cancel out to "feels like 1g" at the equator. Since the game doesn't currently do that with rotational speed, let's
+	// not complicate it here either for now.
+
+	return habitability;
+}
+
+
 /*
  * Set natural resources, tech level, industry strengths and population levels
  */
@@ -1487,8 +1564,16 @@ void PopulateStarSystemGenerator::PopulateStage1(SystemBody *sbody, StarSystem::
 
 	sbody->m_population = fixed();
 
+	// See if there are any surface starports already created. If yes then it's a bit late to say there is no population.
+	fixed surfaceStarportPop = fixed();
+	for (const auto &child : sbody->m_children) {
+		if (child->GetType() == SystemBody::TYPE_STARPORT_SURFACE)
+			surfaceStarportPop += child->GetPopulationAsFixed();
+	}
+
 	/* Bad type of planet for settlement */
-	if ((sbody->GetAverageTemp() > CELSIUS + 100) || (sbody->GetAverageTemp() < 100) || (sbody->CalcSurfaceGravity() > MAX_SETTLEMENT_SURFACE_GRAVITY) ||
+	if ((sbody->GetAverageTemp() > CELSIUS + 100) || (sbody->GetAverageTemp() < 100) ||
+		(sbody->CalcSurfaceGravity() < MIN_SETTLEMENT_SURFACE_GRAVITY) || (sbody->CalcSurfaceGravity() > MAX_SETTLEMENT_SURFACE_GRAVITY) ||
 		(sbody->GetType() != SystemBody::TYPE_PLANET_TERRESTRIAL && sbody->GetType() != SystemBody::TYPE_PLANET_ASTEROID)) {
 		Random starportPopRand;
 		starportPopRand.seed(_init, 6);
@@ -1509,7 +1594,8 @@ void PopulateStarSystemGenerator::PopulateStage1(SystemBody *sbody, StarSystem::
 			}
 		}
 
-		return;
+		if (surfaceStarportPop == 0)
+			return;
 	}
 
 	sbody->m_agricultural = fixed();
@@ -1523,9 +1609,15 @@ void PopulateStarSystemGenerator::PopulateStage1(SystemBody *sbody, StarSystem::
 	} else {
 		// don't bother populating crap planets
 		if (sbody->GetMetallicityAsFixed() < fixed(5, 10) &&
-			sbody->GetMetallicityAsFixed() < (fixed(1, 1) - system->GetHumanProx())) return;
+			sbody->GetMetallicityAsFixed() < (fixed(1, 1) - system->GetHumanProx()) &&
+			(surfaceStarportPop == 0)) return;
 	}
 
+	// Work out the general habitability of the system body, so we can scale the economy according to that.
+	// In general, Earth-like places will have a habitability of 1, and tiny rocky moons will tend towards 0.
+	// This ensures we don't get tiny rocky moons with absurdly high populations and industrial output.
+	const fixed habitability = CalculateHabitability(sbody);
+
 	/* Commodities we produce (mining and agriculture) */
 
 	fixed workforce = fixed();
@@ -1554,6 +1646,8 @@ void PopulateStarSystemGenerator::PopulateStage1(SystemBody *sbody, StarSystem::
 		else
 			affinity = fixed(1, 1);
 
+		affinity *= habitability;
+
 		affinity *= rand.Fixed();
 
 		if (GalacticEconomy::Consumables().count(commodity.id)) {
@@ -1586,8 +1680,15 @@ void PopulateStarSystemGenerator::PopulateStage1(SystemBody *sbody, StarSystem::
 		}
 	}
 
+	// At this point we have the workforce, and m_population is still zero. Now calculate what the population would be.
+
+	// The population is at least as large as the workforce.
 	sbody->m_population += workforce;
 
+	// The more prosperous an economy is, the more non-working people there will be, e.g. children, retirees.
+	// We don't really have a good proxy for that at the moment, so use distance-from-Sol instead.
+	sbody->m_population += sbody->m_population * system->GetHumanProx();
+
 	if (!system->HasCustomBodies() && sbody->GetPopulationAsFixed() > 0)
 		sbody->m_name = Pi::luaNameGen->BodyName(sbody, namerand);
 
@@ -1611,6 +1712,9 @@ void PopulateStarSystemGenerator::PopulateStage1(SystemBody *sbody, StarSystem::
 	// well, outdoor worlds should have way more people
 	sbody->m_population = fixed(1, 10) * sbody->m_population + sbody->m_population * sbody->GetAgriculturalAsFixed();
 
+	// Add in any surface starport populations x2. This ensures that the planet pop is always larger.
+	sbody->m_population += surfaceStarportPop * fixed(2, 1);
+
 	//	Output("%s: pop %.3f billion\n", name.c_str(), sbody->m_population.ToFloat());
 
 	outTotalPop += sbody->GetPopulationAsFixed();

src/galaxy/SystemBody.cpp

@@ -969,6 +969,17 @@ double SystemBody::CalcSurfaceGravity() const
 	}
 }
 
+fixed SystemBody::CalcSurfaceGravityAsFixed() const
+{
+	if (m_radius <= fixed(0, 1))
+		return fixed(0, 1);
+	// m_radius and m_mass are already in native body units
+	// (Earth radius/mass for planets, Sol radius/mass for stars)
+	// so we don't need big G here - just multiply by the relevant value of g.
+	const fixed unitG = (GetSuperType() == SUPERTYPE_STAR) ? G_SOL : G_EARTH;
+	return unitG * ((m_mass / m_radius) / m_radius);
+}
+
 double SystemBody::CalcEscapeVelocity() const
 {
 	double r = GetRadius();

src/galaxy/SystemBody.h

@@ -14,6 +14,9 @@
 #include "galaxy/SystemPath.h"
 #include "gameconsts.h"
 
+static const fixed G_EARTH = fixed(980665, 100000); // 9.80665 m/s^2
+static const fixed G_SOL = fixed(2742, 10);	  // 274.2 m/s^2
+
 class StarSystem;
 
 struct AtmosphereParameters;
@@ -261,6 +264,7 @@ class SystemBody : public RefCounted, public SystemBodyType, protected SystemBod
 	double GetVolatileIces() const { return m_volatileIces.ToDouble(); }
 	fixed GetVolcanicityAsFixed() const { return m_volcanicity; }
 	double GetVolcanicity() const { return m_volcanicity.ToDouble(); }
+	fixed GetAtmosOxidizingAsFixed() const { return m_atmosOxidizing; }
 	double GetAtmosOxidizing() const { return m_atmosOxidizing.ToDouble(); }
 	fixed GetLifeAsFixed() const { return m_life; }
 	double GetLife() const { return m_life.ToDouble(); }
@@ -270,6 +274,7 @@ class SystemBody : public RefCounted, public SystemBodyType, protected SystemBod
 	fixed GetPopulationAsFixed() const { return m_population; }
 	double GetPopulation() const { return m_population.ToDouble(); }
 
+	fixed CalcSurfaceGravityAsFixed() const;
 	double CalcSurfaceGravity() const;
 	double CalcEscapeVelocity() const;
 	double CalcMeanDensity() const;