aissia/docs/02-systems/GEOLOGICAL_SIMULATION_SYSTEM.md

Geological Simulation System

Status: Designed - Ready for Implementation Scope: Complete planetary formation from accretion to industrial-ready geology Duration: 4.65 billion years simulation in 95 cycles

System Overview

Revolutionary geological simulation system that generates realistic planetary geology through scientifically-inspired processes. Creates diverse, coherent geology with proper resource distribution for industrial gameplay.

Key Innovations

• RegionalInfluence framework adapted for geological processes
• Tectonic regions as simple circles with forces (not complex mesh)
• Carbon region system for realistic coal/oil formation
• Dynamic sea level affecting all geological processes
• Optimized 24-byte tiles supporting geological temperature ranges

Simulation Phases

Phase 1: Planetary Accretion (30 cycles × 100M years = 3.0 Ga)

Initial Conditions:

• Temperature: -100°C (space cold)
• Elevation: -30,000m (no surface)
• Empty planetary core

Process - Meteorite Bombardment:

For each wave (100M years):
    1. Generate 100-500 random meteorite impacts
    2. Each impact adds:
       - Kinetic energy → heat (up to +2000°C)
       - Mass → elevation change (+crater formation)
       - Metal composition → resource deposits
    3. Heavy metals sink to planetary core
    4. Core temperature drives volcanic eruptions
    5. Eruptions redistribute core metals to surface
    6. Gradual cooling between waves

Expected Results:

• Surface elevation: -30,000m → -15,000m
• Core composition: 80% iron, 15% nickel, 5% precious metals
• Crater-based geology with metal deposits
• Temperature stabilization around 1000-1500°C

Phase 2: Tectonic Formation (25 cycles × 100M years = 2.5 Ga)

Reduced Meteorite Activity:

// Phase 2 meteorite parameters (reduced from Phase 1)
impacts_per_wave = rng.randInt(10, 50);        // 10x reduction
meteorite_mass = rng.randFloat(10^12, 10^15);  // Smaller impacts
impact_probability = 0.3f;                     // 30% chance per cycle

Tectonic Region System:

struct TectonicRegion {
    uint32_t region_id;
    float center_x, center_y;           // Mobile center (world coordinates)
    float radius;                       // Current radius (tiles)
    float velocity_x, velocity_y;       // Movement (tiles per cycle)
    float mass;                         // For repulsive force calculations
    RegionType type;                    // OCEANIC, CONTINENTAL, VOLCANIC

    // Evolution parameters
    float split_probability;            // Chance to divide (0.01 = 1% per cycle)
    float growth_rate;                  // Radius change per cycle (+/- tiles)
    float stability;                    // Resistance to external forces

    // Formation tracking
    int formation_cycle;                // When this region was created
    uint32_t parent_region_id;          // If split from another region
};

enum class RegionType {
    OCEANIC,        // Dense, subducts under continental
    CONTINENTAL,    // Light, forms mountains when compressed
    VOLCANIC,       // Active volcanic zone (temporary)
};

Tectonic Physics Engine:

void updateTectonicRegions(float time_step) {
    // 1. Calculate repulsive forces between all regions
    for (auto& region1 : tectonic_regions) {
        Vector2 total_force = {0.0f, 0.0f};

        for (auto& region2 : tectonic_regions) {
            if (region1.region_id == region2.region_id) continue;

            float distance = calculateDistance(region1.center, region2.center);
            float overlap = (region1.radius + region2.radius) - distance;

            if (overlap > 0) {
                // Collision detected!
                Vector2 repulsion_direction = normalize(region1.center - region2.center);

                // Force magnitude: F = k * overlap^2 / (mass1 + mass2)
                float force_magnitude = REPULSION_CONSTANT * overlap * overlap /
                                      (region1.mass + region2.mass);

                Vector2 repulsion_force = repulsion_direction * force_magnitude;
                total_force += repulsion_force;

                // Create volcanic zone at collision point
                createVolcanicZone(region1, region2, overlap);
            }
        }

        // Apply acceleration: a = F/m
        Vector2 acceleration = total_force / region1.mass;
        region1.velocity_x += acceleration.x * time_step;
        region1.velocity_y += acceleration.y * time_step;

        // Apply velocity damping (friction)
        region1.velocity_x *= VELOCITY_DAMPING;
        region1.velocity_y *= VELOCITY_DAMPING;
    }

    // 2. Update positions and sizes
    for (auto& region : tectonic_regions) {
        // Move regions
        region.center_x += region.velocity_x * time_step;
        region.center_y += region.velocity_y * time_step;

        // Boundary wrapping (world is spherical)
        wrapCoordinates(region.center_x, region.center_y);

        // Evolve size
        region.radius += region.growth_rate * time_step;
        region.radius = std::clamp(region.radius, MIN_REGION_RADIUS, MAX_REGION_RADIUS);

        // Check for splitting
        if (region.radius > SPLIT_THRESHOLD &&
            rng.probability(region.split_probability * time_step)) {
            splitTectonicRegion(region);
        }
    }
}

Volcanic Zone Creation:

void createVolcanicZone(TectonicRegion& r1, TectonicRegion& r2, float collision_intensity) {
    // Volcanic zone at intersection point
    Vector2 intersection = calculateIntersectionCenter(r1, r2);
    float volcanic_radius = collision_intensity * VOLCANIC_RADIUS_FACTOR;

    // Create temporary volcanic region
    TectonicRegion volcanic_zone = {
        .center_x = intersection.x,
        .center_y = intersection.y,
        .radius = volcanic_radius,
        .type = RegionType::VOLCANIC,
        .formation_cycle = current_cycle,
        .split_probability = 0.05f,  // High chance to break apart
        .growth_rate = -0.1f         // Shrinks over time
    };

    // Apply volcanic influence to affected tiles
    RegionalInfluence volcanic_influence = {
        .type = "volcanic_mountain_formation",
        .intensity = collision_intensity,
        .properties = json{
            {"elevation_bonus", collision_intensity * 500.0f},  // +500m per intensity
            {"temperature_bonus", collision_intensity * 300.0f}, // +300°C
            {"resource_deposits", json::array({"iron", "sulfur", "rare_metals"})},
            {"volcanic_features", json::array({"volcano", "hot_springs", "lava_tubes"})}
        }
    };

    applyRegionalInfluence(volcanic_influence, intersection.x, intersection.y, volcanic_radius);
    tectonic_regions.push_back(volcanic_zone);
}

Region Splitting Algorithm:

void splitTectonicRegion(TectonicRegion& parent) {
    // Create two child regions
    float split_distance = parent.radius * 0.4f;
    Vector2 split_direction = randomUnitVector();

    TectonicRegion child1 = parent;
    TectonicRegion child2 = parent;

    // Assign new IDs
    child1.region_id = next_region_id++;
    child2.region_id = next_region_id++;
    child1.parent_region_id = parent.region_id;
    child2.parent_region_id = parent.region_id;

    // Separate positions
    child1.center_x += split_direction.x * split_distance;
    child1.center_y += split_direction.y * split_distance;
    child2.center_x -= split_direction.x * split_distance;
    child2.center_y -= split_direction.y * split_distance;

    // Reduce sizes
    child1.radius *= 0.7f;
    child2.radius *= 0.7f;

    // Opposite velocities (moving apart)
    child1.velocity_x = split_direction.x * SPLIT_VELOCITY;
    child1.velocity_y = split_direction.y * SPLIT_VELOCITY;
    child2.velocity_x = -split_direction.x * SPLIT_VELOCITY;
    child2.velocity_y = -split_direction.y * SPLIT_VELOCITY;

    // Replace parent with children
    removeRegion(parent);
    addRegion(child1);
    addRegion(child2);

    // Create rift valley between children
    createRiftValley(child1, child2);
}

Expected Results:

• Surface elevation: -15,000m → -5,000m (+10km crustal thickening)
• Formation of 15-25 stable tectonic regions
• Mountain ranges at collision zones (+1000-3000m elevation)
• Rift valleys where regions separate (-500m elevation)
• Distinct continental vs oceanic regions

Phase 3: Hydrological Cycles (25 cycles × 20M years = 0.5 Ga)

Dynamic Sea Level System:

struct OceanLevel {
    float current_sea_level;        // Global reference (meters)
    float ice_volume;               // Polar ice volume (km³)
    float thermal_expansion_factor; // Ocean thermal expansion coefficient
    float tectonic_basin_volume;    // Total oceanic basin volume

    void updateSeaLevel(float global_temperature, float volcanic_co2) {
        // 1. Ice volume changes (glaciation/melting)
        float target_ice_volume = calculateIceVolume(global_temperature);
        float ice_change = (target_ice_volume - ice_volume) * 0.1f; // 10% change per cycle
        ice_volume += ice_change;

        // Sea level change: 1 km³ ice = +2.7mm sea level
        float ice_effect = -ice_change * 0.0027f;

        // 2. Thermal expansion
        float thermal_effect = (global_temperature - 15.0f) * thermal_expansion_factor;

        // 3. Tectonic effects (basin formation/destruction)
        float tectonic_effect = (tectonic_basin_volume - baseline_basin_volume) * 0.001f;

        // 4. Update sea level
        current_sea_level += ice_effect + thermal_effect + tectonic_effect;
        current_sea_level = std::clamp(current_sea_level, -200.0f, +100.0f); // Realistic bounds
    }

    float calculateIceVolume(float global_temp) const {
        if (global_temp < -5.0f) return MAX_ICE_VOLUME;      // Ice age
        if (global_temp > 25.0f) return 0.0f;               // No ice
        return MAX_ICE_VOLUME * (25.0f - global_temp) / 30.0f; // Linear interpolation
    }
};

Hydraulic Erosion System:

void applyHydraulicErosion(float cycle_duration_years) {
    const float HYDRAULIC_EFFICIENCY = 10.0f; // 10x more effective than tectonic

    for (int y = 0; y < map_height; ++y) {
        for (int x = 0; x < map_width; ++x) {
            WorldTile tile = world_map->getTile(x, y);
            float elevation = tile.getElevation();

            if (elevation > current_sea_level) {
                // CONTINENTAL EROSION

                // 1. Calculate water flow (slope-based)
                float average_neighbor_elevation = calculateAverageNeighborElevation(x, y);
                float slope = elevation - average_neighbor_elevation;
                float water_flow = std::max(0.0f, slope * 0.01f); // Flow intensity

                // 2. Erosion rate based on elevation and flow
                float elevation_factor = (elevation - current_sea_level) / 1000.0f; // Higher = more erosion
                float erosion_rate = water_flow * elevation_factor * HYDRAULIC_EFFICIENCY;

                // 3. Apply erosion
                float erosion_amount = erosion_rate * cycle_duration_years / 1000000.0f; // Per million years
                tile.setElevation(elevation - erosion_amount);

                // 4. Transport sediments downstream
                transportSediments(tile, x, y, erosion_amount);

            } else if (elevation > current_sea_level - 200.0f) {
                // COASTAL EROSION

                float depth_below_sea = current_sea_level - elevation;
                float coastal_erosion_rate = 50.0f; // Intense coastal erosion

                if (depth_below_sea < 50.0f) { // Shallow coastal zone
                    float erosion_amount = coastal_erosion_rate * cycle_duration_years / 1000000.0f;
                    tile.setElevation(elevation - erosion_amount);

                    // Form beaches and deltas
                    if (hasRiverInput(x, y)) {
                        formDelta(tile, x, y);
                    }
                }

            } else {
                // DEEP OCEAN - sediment deposition

                float sedimentation_rate = calculateSedimentInput(x, y);
                float deposition_amount = sedimentation_rate * cycle_duration_years / 1000000.0f;
                tile.setElevation(elevation + deposition_amount);
            }
        }
    }
}

void transportSediments(WorldTile& source_tile, int x, int y, float sediment_amount) {
    // Find steepest downhill direction
    auto [target_x, target_y] = findSteepestDescentDirection(x, y);

    if (target_x != x || target_y != y) {
        WorldTile target_tile = world_map->getTile(target_x, target_y);

        // Deposit sediments at target
        float current_elevation = target_tile.getElevation();
        target_tile.setElevation(current_elevation + sediment_amount * 0.5f); // 50% deposition

        // Create river network
        createRiverSegment(x, y, target_x, target_y);

        // Continue transport if target is above sea level
        if (target_tile.getElevation() > current_sea_level) {
            transportSediments(target_tile, target_x, target_y, sediment_amount * 0.5f);
        }
    }
}

River Network Formation:

struct RiverNetwork {
    std::unordered_map<uint32_t, std::vector<uint32_t>> river_connections; // tile -> downstream tiles
    std::unordered_map<uint32_t, float> flow_volumes; // tile -> water volume

    void createRiverSegment(int from_x, int from_y, int to_x, int to_y) {
        uint32_t from_index = from_y * map_width + from_x;
        uint32_t to_index = to_y * map_width + to_x;

        // Add connection
        river_connections[from_index].push_back(to_index);

        // Increase flow volume
        flow_volumes[to_index] += flow_volumes[from_index] + 1.0f;

        // Mark tiles as having river features
        WorldTile from_tile = world_map->getTile(from_x, from_y);
        WorldTile to_tile = world_map->getTile(to_x, to_y);

        from_tile.setFlag(TileFlags::HAS_RIVER, true);
        to_tile.setFlag(TileFlags::HAS_RIVER, true);
    }

    void formDelta(WorldTile& tile, int x, int y) {
        // Delta formation where river meets ocean
        RegionalInfluence delta_influence = {
            .type = "river_delta",
            .intensity = flow_volumes[y * map_width + x] / 100.0f, // Intensity based on river size
            .properties = json{
                {"elevation_bonus", 10.0f}, // Slight elevation increase
                {"sediment_richness", 1.5f}, // Rich sediments
                {"features", json::array({"wetlands", "fertile_soil", "natural_harbors"})}
            }
        };

        applyRegionalInfluence(delta_influence, x, y, 5.0f); // 5-tile radius
    }
};

Climate Integration:

void updateGlobalClimate(float cycle_duration) {
    // 1. Volcanic CO2 emissions
    float volcanic_co2 = calculateVolcanicCO2Emissions();
    atmospheric_co2 += volcanic_co2;

    // 2. Weathering CO2 absorption
    float weathering_absorption = calculateWeatheringRate() * cycle_duration;
    atmospheric_co2 -= weathering_absorption;

    // 3. Temperature from CO2 (greenhouse effect)
    global_temperature = BASE_TEMPERATURE + log(atmospheric_co2 / BASELINE_CO2) * CLIMATE_SENSITIVITY;

    // 4. Update sea level based on new temperature
    ocean_level.updateSeaLevel(global_temperature, volcanic_co2);

    // 5. Regional climate effects
    updateRegionalClimate();
}

void updateRegionalClimate() {
    for (int y = 0; y < map_height; ++y) {
        for (int x = 0; x < map_width; ++x) {
            WorldTile tile = world_map->getTile(x, y);
            float elevation = tile.getElevation();

            // Temperature varies with elevation (lapse rate: -6.5°C/km)
            float local_temperature = global_temperature - (elevation / 1000.0f) * 6.5f;

            // Distance from ocean affects temperature (continental effect)
            float ocean_distance = calculateDistanceToOcean(x, y);
            float continental_effect = ocean_distance / 1000.0f * 2.0f; // +2°C per 1000km inland
            local_temperature += continental_effect;

            tile.setTemperature(local_temperature);

            // Humidity based on proximity to water and temperature
            float humidity = calculateHumidity(ocean_distance, local_temperature, elevation);
            tile.setHumidity(humidity);
        }
    }
}

Expected Results:

• Surface elevation: -5,000m → 0m (modern sea level)
• Formation of permanent river networks draining to oceans
• Coastal features: deltas, beaches, fjords, cliffs
• Climate zones established based on elevation and ocean proximity
• Sedimentary basins filled with eroded material

Phase 4: Carboniferous Period (35 cycles × 10M years = 0.35 Ga)

Dynamic Forest System:

// Forest seeds are continuously created and destroyed by geological events
void manageDynamicForests(GMap& map, int cycle) {
    // Always seed new forests with base probability
    seedNewForests(map);

    // Existing forests attempt to expand
    expandExistingForests(map);

    // Geological events destroy forests
    destroyForestsByGeologicalEvents(map, cycle);
}

void seedNewForests(GMap& map) {
    for (each_tile) {
        float seed_probability = base_forest_seed_rate;

        // Environmental modifiers
        if (tile.getTemperature() >= 15.0f && tile.getTemperature() <= 35.0f) {
            seed_probability *= temperature_bonus;  // Optimal temperature
        }
        if (tile.getHumidity() > 0.4f) {
            seed_probability *= humidity_bonus;     // Sufficient moisture
        }
        if (tile.getElevation() < 2000.0f) {
            seed_probability *= elevation_bonus;    // Below treeline
        }

        if (random() < seed_probability) {
            createForestSeed(tile);
        }
    }
}

void expandExistingForests(GMap& map) {
    for (each_forest_tile) {
        for (each_neighbor) {
            if (isViableForForest(neighbor) && random() < forest_expansion_rate) {
                expandForestTo(neighbor);
            }
        }
    }
}

void destroyForestsByGeologicalEvents(GMap& map, int cycle) {
    // Volcanic eruptions destroy forests
    for (auto& volcanic_event : getCurrentVolcanicEvents(cycle)) {
        destroyForestsInRadius(volcanic_event.center, volcanic_event.destruction_radius);
    }

    // Major floods destroy forests
    for (each_tile_with_extreme_water_level) {
        if (tile.getWaterLevel() > flood_destruction_threshold) {
            destroyForest(tile);
        }
    }

    // Extreme erosion destroys forests
    for (each_tile) {
        if (getCurrentErosionRate(tile) > erosion_destruction_threshold) {
            destroyForest(tile);
        }
    }

    // Climate extremes destroy forests
    if (tile.getTemperature() < -10.0f || tile.getTemperature() > 50.0f) {
        destroyForest(tile);
    }
}

Carbon Region System:

struct CarbonRegion {
    uint32_t region_id;
    float center_x, center_y;           // Position (world coordinates)
    float radius;                       // Affected area (tiles)
    float carbon_mass;                  // Total carbon content (megatons)
    int formation_cycle;                // When this region was created
    CarbonType type;                    // Current state: COAL, OIL, GAS

    // Tectonic attachment
    uint32_t attached_tectonic_region_id;    // Moves with tectonic plate
    float attachment_strength;               // 0.8-1.0 (strong attachment)

    // Conversion tracking
    float original_coal_mass;               // Initial coal amount
    float conversion_progress;              // 0.0-1.0 (coal→oil conversion)
    int cycles_underwater;                  // How long submerged
};

enum class CarbonType {
    COAL,           // Terrestrial formation
    OIL,            // Marine conversion
    NATURAL_GAS     // Late-stage oil maturation
};

Forest Growth and Coal Formation:

void processForestGrowth(int cycle) {
    for (int y = 0; y < map_height; ++y) {
        for (int x = 0; x < map_width; ++x) {
            WorldTile tile = world_map->getTile(x, y);

            if (isForestSuitable(tile, x, y)) {
                float biomass_potential = calculateBiomassPotential(tile);
                float forest_density = std::min(1.0f, biomass_potential);

                if (forest_density > 0.3f) { // Minimum density for coal formation
                    // Create or enhance carbon region
                    CarbonRegion* existing_region = findNearestCarbonRegion(x, y, 5.0f);

                    if (existing_region) {
                        // Add to existing region
                        enhanceCarbonRegion(*existing_region, forest_density);
                    } else {
                        // Create new carbon region
                        CarbonRegion new_region = createCarbonRegion(x, y, forest_density, cycle);
                        carbon_regions.push_back(new_region);
                    }
                }
            }
        }
    }

    // Merge nearby carbon regions
    mergeCarbonRegions();
}

bool isForestSuitable(const WorldTile& tile, int x, int y) {
    float elevation = tile.getElevation();
    float temperature = tile.getTemperature();
    float humidity = tile.getHumidity();

    return (elevation > current_sea_level + 10.0f) &&     // Above sea level
           (elevation < current_sea_level + 1000.0f) &&   // Not too mountainous
           (temperature >= 10.0f && temperature <= 30.0f) && // Temperate range
           (humidity > 0.4f) &&                           // Sufficient moisture
           (!tile.getFlag(TileFlags::VOLCANIC_ACTIVE));   // Not volcanically active
}

CarbonRegion createCarbonRegion(int x, int y, float forest_density, int cycle) {
    CarbonRegion region;
    region.region_id = next_carbon_region_id++;
    region.center_x = x;
    region.center_y = y;
    region.radius = forest_density * 8.0f; // Dense forests = larger regions
    region.carbon_mass = forest_density * 100.0f; // Base carbon mass (megatons)
    region.formation_cycle = cycle;
    region.type = CarbonType::COAL;

    // Attach to nearest tectonic region
    TectonicRegion* nearest_tectonic = findNearestTectonicRegion(x, y);
    if (nearest_tectonic) {
        region.attached_tectonic_region_id = nearest_tectonic->region_id;
        region.attachment_strength = 0.9f; // Strong attachment
    }

    return region;
}

Carbon Region Movement and Merging:

void updateCarbonRegionMovement(float time_step) {
    for (auto& carbon_region : carbon_regions) {
        TectonicRegion* tectonic = getTectonicRegion(carbon_region.attached_tectonic_region_id);

        if (tectonic) {
            // Move with tectonic plate
            carbon_region.center_x += tectonic->velocity_x * carbon_region.attachment_strength * time_step;
            carbon_region.center_y += tectonic->velocity_y * carbon_region.attachment_strength * time_step;

            // Handle tectonic collisions
            if (tectonic->isInCollision()) {
                float collision_stress = tectonic->getCollisionIntensity();

                if (collision_stress > 2.0f) {
                    // High stress can redistribute carbon deposits
                    redistributeCarbonDeposit(carbon_region, collision_stress);
                }
            }
        }
    }
}

void mergeCarbonRegions() {
    for (size_t i = 0; i < carbon_regions.size(); ++i) {
        for (size_t j = i + 1; j < carbon_regions.size(); ++j) {
            float distance = calculateDistance(carbon_regions[i].center, carbon_regions[j].center);
            float merge_threshold = (carbon_regions[i].radius + carbon_regions[j].radius) * 0.8f;

            if (distance < merge_threshold) {
                // Merge regions
                CarbonRegion merged;
                merged.region_id = next_carbon_region_id++;
                merged.center_x = (carbon_regions[i].center_x + carbon_regions[j].center_x) / 2.0f;
                merged.center_y = (carbon_regions[i].center_y + carbon_regions[j].center_y) / 2.0f;
                merged.radius = sqrt(carbon_regions[i].radius * carbon_regions[i].radius +
                                   carbon_regions[j].radius * carbon_regions[j].radius);
                merged.carbon_mass = carbon_regions[i].carbon_mass + carbon_regions[j].carbon_mass;
                merged.formation_cycle = std::min(carbon_regions[i].formation_cycle, carbon_regions[j].formation_cycle);
                merged.type = CarbonType::COAL; // Combined regions start as coal

                // Remove original regions and add merged
                carbon_regions.erase(carbon_regions.begin() + j);
                carbon_regions.erase(carbon_regions.begin() + i);
                carbon_regions.push_back(merged);

                break; // Start over after modification
            }
        }
    }
}

Coal to Oil Conversion:

void processCoalToOilConversion(int current_cycle) {
    const float CONVERSION_RATE_PER_CYCLE = 0.05f; // 5% per cycle (10M years)
    const float NATURAL_GAS_RATIO = 0.15f;         // 15% of oil becomes gas

    for (auto& coal_region : carbon_regions) {
        if (coal_region.type == CarbonType::COAL && isRegionUnderwater(coal_region)) {
            coal_region.cycles_underwater++;

            // Only convert if underwater for multiple cycles (pressure + time)
            if (coal_region.cycles_underwater >= 3) {
                float conversion_amount = coal_region.carbon_mass * CONVERSION_RATE_PER_CYCLE;

                if (conversion_amount > 1.0f) { // Minimum threshold
                    // Create oil region at same location
                    CarbonRegion oil_region = coal_region;
                    oil_region.region_id = next_carbon_region_id++;
                    oil_region.type = CarbonType::OIL;
                    oil_region.carbon_mass = conversion_amount;
                    oil_region.radius *= 0.8f; // Oil regions are more concentrated
                    oil_region.formation_cycle = current_cycle;

                    // Create natural gas as byproduct
                    if (conversion_amount > 10.0f) {
                        CarbonRegion gas_region = oil_region;
                        gas_region.region_id = next_carbon_region_id++;
                        gas_region.type = CarbonType::NATURAL_GAS;
                        gas_region.carbon_mass = conversion_amount * NATURAL_GAS_RATIO;
                        gas_region.radius *= 1.2f; // Gas spreads more widely

                        carbon_regions.push_back(gas_region);
                    }

                    // Reduce original coal region
                    coal_region.carbon_mass -= conversion_amount;
                    coal_region.conversion_progress += CONVERSION_RATE_PER_CYCLE;

                    // Add oil region
                    carbon_regions.push_back(oil_region);

                    // Remove depleted coal regions
                    if (coal_region.carbon_mass < 1.0f) {
                        // Mark for removal
                        coal_region.carbon_mass = 0.0f;
                    }
                }
            }
        }
    }

    // Remove depleted regions
    carbon_regions.erase(
        std::remove_if(carbon_regions.begin(), carbon_regions.end(),
                      [](const CarbonRegion& region) { return region.carbon_mass < 1.0f; }),
        carbon_regions.end()
    );
}

bool isRegionUnderwater(const CarbonRegion& region) {
    // Check if majority of region is below sea level
    int underwater_tiles = 0;
    int total_tiles = 0;

    int radius_int = static_cast<int>(region.radius);
    for (int dy = -radius_int; dy <= radius_int; ++dy) {
        for (int dx = -radius_int; dx <= radius_int; ++dx) {
            float distance = sqrt(dx*dx + dy*dy);
            if (distance <= region.radius) {
                int x = static_cast<int>(region.center_x) + dx;
                int y = static_cast<int>(region.center_y) + dy;

                if (world_map->isValidCoordinate(x, y)) {
                    WorldTile tile = world_map->getTile(x, y);
                    total_tiles++;

                    if (tile.getElevation() < current_sea_level - 50.0f) { // 50m underwater minimum
                        underwater_tiles++;
                    }
                }
            }
        }
    }

    return (static_cast<float>(underwater_tiles) / total_tiles) > 0.6f; // 60% underwater
}

Resource Application to Tiles:

void applyCarbonRegionsToTiles() {
    for (const auto& region : carbon_regions) {
        int radius_int = static_cast<int>(region.radius);

        for (int dy = -radius_int; dy <= radius_int; ++dy) {
            for (int dx = -radius_int; dx <= radius_int; ++dx) {
                float distance = sqrt(dx*dx + dy*dy);
                if (distance <= region.radius) {
                    int x = static_cast<int>(region.center_x) + dx;
                    int y = static_cast<int>(region.center_y) + dy;

                    if (world_map->isValidCoordinate(x, y)) {
                        WorldTile tile = world_map->getTile(x, y);
                        float influence = 1.0f - (distance / region.radius); // Gradient from center

                        // Apply resource based on carbon type
                        switch (region.type) {
                            case CarbonType::COAL:
                                addResourceToTile(tile, "coal", region.carbon_mass * influence);
                                break;
                            case CarbonType::OIL:
                                addResourceToTile(tile, "petroleum", region.carbon_mass * influence);
                                break;
                            case CarbonType::NATURAL_GAS:
                                addResourceToTile(tile, "natural_gas", region.carbon_mass * influence);
                                break;
                        }
                    }
                }
            }
        }
    }
}

Expected Results:

• Dynamic forest evolution: Forests continuously seed, expand, and get destroyed by geological events
• Multiple forest generations: Layers of carbon deposits from different geological periods
• Realistic forest patterns: Forests avoid active volcanic zones, extreme climates, and flood-prone areas
• Coal deposit diversity: Various ages and qualities of coal from different forest cycles
• Oil formation: Submerged forest areas naturally convert to petroleum over time
• Geological storytelling: Each coal seam represents a specific period of forest growth and burial

Phase 5: Pre-Faunal Stabilization (15 cycles × 10M years = 0.15 Ga)

Maturation-Only Phase - No New Formation:

void processStabilizationCycle(int cycle) {
    // 1. HYDROCARBON MATURATION (no new formation)
    continueHydrocarbonMaturation();

    // 2. PETROLEUM MIGRATION TO TRAPS
    migratePetroleumToGeologicalTraps();

    // 3. FINAL EROSION PHASE
    applyFinalErosion();

    // 4. SOIL LAYER DEVELOPMENT
    developSoilLayers();

    // 5. CLIMATE STABILIZATION
    stabilizeClimate();

    // 6. GEOLOGICAL STRUCTURE FINALIZATION
    finalizeGeologicalStructures();
}

Continued Hydrocarbon Maturation:

void continueHydrocarbonMaturation() {
    const float REDUCED_CONVERSION_RATE = 0.03f; // Slower than Carboniferous (3% per cycle)

    for (auto& coal_region : carbon_regions) {
        if (coal_region.type == CarbonType::COAL && isRegionUnderwater(coal_region)) {
            // Continue coal→oil conversion at reduced rate
            float conversion_amount = coal_region.carbon_mass * REDUCED_CONVERSION_RATE;

            if (conversion_amount > 0.5f) { // Lower threshold for final conversion
                processCoalToOilConversion(coal_region, conversion_amount);
            }
        }
    }

    // Oil→Gas maturation for old oil deposits
    for (auto& oil_region : carbon_regions) {
        if (oil_region.type == CarbonType::OIL) {
            int oil_age = current_cycle - oil_region.formation_cycle;

            if (oil_age > 10 && oil_region.carbon_mass > 5.0f) { // Mature oil deposits
                float gas_conversion = oil_region.carbon_mass * 0.02f; // 2% oil→gas

                createNaturalGasFromOil(oil_region, gas_conversion);
            }
        }
    }
}

Petroleum Migration to Geological Traps:

struct GeologicalTrap {
    float center_x, center_y;
    float capacity;                    // Maximum hydrocarbon storage
    float current_fill;               // Current hydrocarbon content
    TrapType type;                    // ANTICLINE, FAULT, SALT_DOME, STRATIGRAPHIC
    float formation_cycle;            // When trap was formed
};

enum class TrapType {
    ANTICLINE,        // Upward fold in rock layers
    FAULT_TRAP,       // Hydrocarbon blocked by fault
    SALT_DOME,        // Hydrocarbon trapped around salt intrusion
    STRATIGRAPHIC     // Trapped between different rock types
};

void migratePetroleumToGeologicalTraps() {
    // 1. Identify geological traps from tectonic history
    std::vector<GeologicalTrap> traps = identifyGeologicalTraps();

    // 2. Migrate oil and gas to nearest suitable traps
    for (auto& hydrocarbon_region : carbon_regions) {
        if (hydrocarbon_region.type == CarbonType::OIL || hydrocarbon_region.type == CarbonType::NATURAL_GAS) {
            GeologicalTrap* nearest_trap = findNearestTrap(hydrocarbon_region, traps);

            if (nearest_trap && nearest_trap->current_fill < nearest_trap->capacity) {
                float migration_amount = calculateMigrationAmount(hydrocarbon_region, *nearest_trap);

                // Move hydrocarbons to trap
                nearest_trap->current_fill += migration_amount;
                hydrocarbon_region.carbon_mass -= migration_amount;

                // Concentrate hydrocarbons in trap location
                concentrateHydrocarbonsAtTrap(*nearest_trap, hydrocarbon_region.type, migration_amount);
            }
        }
    }
}

std::vector<GeologicalTrap> identifyGeologicalTraps() {
    std::vector<GeologicalTrap> traps;

    // Find anticlines (upward folds) from tectonic compression
    for (const auto& tectonic_region : tectonic_regions) {
        if (tectonic_region.type == RegionType::CONTINENTAL) {
            // Look for elevation highs within the region (potential anticlines)
            auto high_points = findElevationHighsInRegion(tectonic_region);

            for (const auto& point : high_points) {
                GeologicalTrap anticline = {
                    .center_x = point.x,
                    .center_y = point.y,
                    .capacity = calculateAnticlineCapacity(point),
                    .current_fill = 0.0f,
                    .type = TrapType::ANTICLINE,
                    .formation_cycle = tectonic_region.formation_cycle
                };
                traps.push_back(anticline);
            }
        }
    }

    // Find fault traps from tectonic collisions
    for (const auto& collision_zone : historical_tectonic_collisions) {
        GeologicalTrap fault_trap = {
            .center_x = collision_zone.center_x,
            .center_y = collision_zone.center_y,
            .capacity = collision_zone.intensity * 50.0f, // Capacity based on collision intensity
            .current_fill = 0.0f,
            .type = TrapType::FAULT_TRAP,
            .formation_cycle = collision_zone.formation_cycle
        };
        traps.push_back(fault_trap);
    }

    return traps;
}

Final Erosion and Surface Formation:

void applyFinalErosion() {
    const float FINAL_EROSION_RATE = 5.0f; // Reduced from active hydrological phase

    for (int y = 0; y < map_height; ++y) {
        for (int x = 0; x < map_width; ++x) {
            WorldTile tile = world_map->getTile(x, y);
            float elevation = tile.getElevation();

            if (elevation > current_sea_level) {
                // Expose coal seams through erosion
                exposeCoalSeams(tile, x, y);

                // Carve final river valleys
                carveRiverValleys(tile, x, y);

                // Form final coastal features
                if (elevation < current_sea_level + 100.0f) { // Coastal zone
                    formCoastalFeatures(tile, x, y);
                }
            }
        }
    }
}

void exposeCoalSeams(WorldTile& tile, int x, int y) {
    // Check if there are coal deposits below current elevation
    CarbonRegion* coal_region = findCoalRegionAt(x, y);

    if (coal_region && coal_region->type == CarbonType::COAL) {
        float erosion_depth = calculateErosionDepth(tile);

        if (erosion_depth > 50.0f) { // Significant erosion
            // Expose coal seam at surface
            tile.setFlag(TileFlags::SURFACE_COAL, true);

            // Add surface coal resource
            addResourceToTile(tile, "surface_coal", coal_region->carbon_mass * 0.1f);
        }
    }
}

Soil Development System:

void developSoilLayers() {
    for (int y = 0; y < map_height; ++y) {
        for (int x = 0; x < map_width; ++x) {
            WorldTile tile = world_map->getTile(x, y);

            if (tile.getElevation() > current_sea_level + 5.0f) { // Above sea level
                SoilType soil = determineSoilType(tile, x, y);
                float soil_depth = calculateSoilDepth(tile, x, y);

                // Apply soil influence
                RegionalInfluence soil_influence = {
                    .type = "soil_development",
                    .intensity = soil_depth / 10.0f, // Depth in meters / 10
                    .properties = json{
                        {"soil_type", getSoilTypeName(soil)},
                        {"fertility", calculateSoilFertility(soil, tile)},
                        {"drainage", calculateSoilDrainage(soil, tile)},
                        {"ph_level", calculateSoilPH(soil, tile)}
                    }
                };

                applyRegionalInfluence(soil_influence, x, y, 1.0f);
            }
        }
    }
}

enum class SoilType {
    CLAY,           // Heavy, nutrient-rich, poor drainage
    SAND,           // Light, good drainage, low nutrients
    LOAM,           // Balanced, ideal for agriculture
    PEAT,           // Organic-rich, acidic, wetland areas
    ROCKY,          // Thin soil, mountainous areas
    ALLUVIAL        // River deposits, very fertile
};

SoilType determineSoilType(const WorldTile& tile, int x, int y) {
    float elevation = tile.getElevation();
    float temperature = tile.getTemperature();
    float humidity = tile.getHumidity();

    // River deltas and floodplains
    if (tile.getFlag(TileFlags::HAS_RIVER) && elevation < current_sea_level + 50.0f) {
        return SoilType::ALLUVIAL;
    }

    // Wetland areas
    if (humidity > 0.8f && elevation < current_sea_level + 20.0f) {
        return SoilType::PEAT;
    }

    // Mountainous areas
    if (elevation > current_sea_level + 1000.0f) {
        return SoilType::ROCKY;
    }

    // Climate-based soil formation
    if (temperature > 20.0f && humidity < 0.3f) {
        return SoilType::SAND; // Arid regions
    } else if (temperature < 10.0f && humidity > 0.6f) {
        return SoilType::CLAY; // Cold, wet regions
    } else {
        return SoilType::LOAM; // Temperate regions
    }
}

Final Climate Stabilization:

void stabilizeClimate() {
    // Reduce climate variability and volcanic activity
    volcanic_activity_level *= 0.8f; // 20% reduction per cycle
    climate_stability += 0.1f;       // Increase stability

    // Establish final climate zones
    for (int y = 0; y < map_height; ++y) {
        for (int x = 0; x < map_width; ++x) {
            WorldTile tile = world_map->getTile(x, y);

            ClimateZone zone = determineClimateZone(tile, x, y);
            applyClimateZone(tile, zone);
        }
    }
}

enum class ClimateZone {
    ARCTIC,         // < -10°C, low precipitation
    SUBARCTIC,      // -10 to 0°C, moderate precipitation
    TEMPERATE,      // 0 to 20°C, variable precipitation
    SUBTROPICAL,    // 20 to 30°C, high precipitation
    TROPICAL,       // > 30°C, very high precipitation
    ARID,           // Any temperature, very low precipitation
    MEDITERRANEAN   // Warm, dry summers, mild wet winters
};

Final Geological State:

struct FinalGeologicalState {
    // TERRAIN
    ✅ stable_continents;              // Continental masses established
    ✅ ocean_basins;                   // Deep ocean basins
    ✅ mountain_ranges;                // Various ages, realistic erosion
    ✅ river_networks;                 // Mature drainage systems
    ✅ coastal_features;               // Beaches, cliffs, deltas, fjords

    // RESOURCES
    ✅ coal_deposits;                  // Continental basins, exposed seams
    ✅ oil_fields;                     // Sedimentary basins, geological traps
    ✅ natural_gas;                    // Associated with oil, separate fields
    ✅ metal_deposits;                 // From meteorite impacts and volcanism

    // CLIMATE
    ✅ climate_zones;                  // Stable temperature/precipitation patterns
    ✅ soil_types;                     // Mature soil development
    ✅ seasonal_patterns;              // Established weather cycles

    // READY FOR INDUSTRIAL GAMEPLAY
    ✅ resource_accessibility;         // Surface coal, shallow oil, metal ores
    ✅ transportation_routes;          // Rivers, coastal access, terrain variety
    ✅ strategic_locations;            // Resource clusters, defensive positions
    ✅ environmental_challenges;       // Climate zones, terrain obstacles
};

Expected Results:

• Complete geological maturity: All major processes stabilized
• Industrial-ready resources: Coal seams exposed, oil in accessible traps
• Realistic geography: Mountain ranges, river valleys, coastal plains
• Climate diversity: Multiple biomes and environmental conditions
• Strategic complexity: Resource distribution creates interesting gameplay choices
• Historical coherence: Geological features tell the story of planetary formation

Phase 6: Climate Simulation and Biome Generation

Overview

After geological stabilization, run advanced climate simulation using mobile WindRegions and Inter-Tropical Convergence Zones (ITCZ) to establish realistic temperature, humidity, and wind patterns. Creates emergent weather patterns through wind region interactions rather than predefined climate zones.

Key Innovation: Revolutionary climate simulation using mobile wind regions that spawn, evolve, and interact to create emergent weather patterns. Solves the "Sahara vs Congo" problem through Inter-Tropical Convergence Zones (ITCZ) and planetary rotation bands without complex 3D atmospheric physics.

Climate Simulation Architecture

Phase 6.1: Landmass Analysis and ITCZ Generation

Uses existing TectonicRegions from Phase 2:

// Analyze continental masses from existing tectonic data
std::vector<Continent> continents = groupTectonicRegions();

// Generate water masses by inverse analysis
std::vector<WaterMass> oceans = detectOceanBasins(continents);

// Place ITCZ zones on qualifying equatorial landmasses
for (auto& continent : continents) {
    if (continent.latitude >= 0.45 && continent.latitude <= 0.55 &&
        continent.area > MIN_LANDMASS_SIZE) {
        createITCZ(continent.center, sqrt(continent.area));
    }
}

ITCZ Requirements:

• Latitude Band: 45-55% of map height (equatorial)
• Minimum Landmass: 10,000+ tiles
• Ocean Proximity: Within 800km for moisture source
• Continental Heating: Large thermal mass for convection

Phase 6.2: WindRegion Spawning System

Mobile WindRegions spawn from ocean masses:

"wind_spawn_system": {
  "spawn_locations": "ocean_masses_only",
  "spawn_frequency": "water_mass_size / 1000",
  "initial_strength": "water_temperature * evaporation_factor",
  "spawn_distribution": "random_within_water_region_biased_toward_center"
}

WindRegion Mobile Entities:

"wind_region": {
  "position": [x, y],
  "wind_strength": 1.0,        // Base intensity (decays over time)
  "wetness": 0.0,              // Moisture content (gained over ocean)
  "velocity": [vx, vy],        // Movement vector
  "wind_tokens": 100,          // Distributed to tiles
  "decay_per_move": 0.02,      // -2% strength per movement
  "decay_per_tile": 0.01       // -1% strength per tile crossed
}

Phase 6.3: Movement and Planetary Circulation

Planetary rotation bands based on real atmospheric data:

"planetary_circulation": {
  "polar_jet_north": {
    "latitude_range": [0.10, 0.25],
    "direction": "west_to_east",
    "strength": 1.8
  },
  "equatorial_trades": {
    "latitude_range": [0.40, 0.60],
    "direction": "east_to_west",
    "strength": 1.5
  },
  "polar_jet_south": {
    "latitude_range": [0.75, 0.95],
    "direction": "west_to_east",
    "strength": 1.8
  }
}

Movement Calculation:

Vector2 movement =
  planetary_rotation_band * 0.6 +     // 60% planetary circulation
  equator_to_pole_bias * 0.2 +        // 20% thermal circulation
  terrain_deflection * 0.1 +          // 10% topographic influence
  random_variation * 0.1;             // 10% chaos

Phase 6.4: WindRegion Evolution and Interactions

Dynamic evolution with ITCZ gravitational effects:

void updateWindRegion(WindRegion& region) {
    // Gain moisture over water
    if (currentTile.isOcean()) {
        region.wetness += OCEAN_MOISTURE_GAIN;
    }

    // Repulsion from other regions = acceleration
    // NOTE: Not physical repulsion - proxy for spatial competition and turbulence
    // Prevents region stacking while creating realistic dispersion patterns
    // CRITIQUE POINT: May cause "force field" effect around ITCZ zones where regions
    // oscillate/scatter instead of converging due to attraction vs repulsion conflict.
    // Alternative approaches: density-based drift, no interaction, or collision division.
    // TODO: Implement as configurable algorithm options for empirical testing.
    float repulsion = calculateRepulsionForce(region, nearbyRegions);
    region.wind_strength += repulsion * ACCELERATION_FACTOR;

    // Movement decay
    region.wind_strength *= (1.0 - DECAY_PER_MOVE);

    // Die when too weak
    if (region.wind_strength < MINIMUM_THRESHOLD) {
        destroyRegion(region);
    }
}

ITCZ Gravitational Effects:

void applyITCZGravity(WindRegion& region) {
    for (auto& itcz : active_itcz_zones) {
        float distance = calculateDistance(region.position, itcz.center);

        if (distance < itcz.gravitational_range) {
            // Attraction force (inverse square law)
            // NOTE: "Gravitational" metaphor for influence strength, not literal physics
            // Like saying someone has "gravitas" - clear semantic meaning for developers
            float attraction = itcz.mass / (distance * distance);
            Vector2 pull_direction = normalize(itcz.center - region.position);

            // Apply attraction
            region.velocity += pull_direction * attraction;

            // Amplification effect as region approaches
            float proximity = (itcz.range - distance) / itcz.range;
            float amplification = 1.0 + (itcz.max_amplification * proximity);

            region.wind_strength *= amplification;
            region.wetness *= amplification;
        }
    }
}

Phase 6.5: Token Distribution and Climate Zone Formation

Climate zone classification using token accumulation:

void distributeTokens(WindRegion& region) {
    // Basic climate tokens for all regions
    int wind_tokens = static_cast<int>(region.wind_strength * 10);
    int rain_tokens = static_cast<int>(region.wetness * 10);

    WorldTile& tile = world_map.getTile(region.position);
    tile.addTokens("wind", wind_tokens);
    tile.addTokens("rain", rain_tokens);

    // Special climate zone tokens for extreme weather
    if (isHighWindZone(region)) {
        tile.addTokens("highWind", 1);  // Hostile to forests
    }
    if (isFloodZone(region)) {
        tile.addTokens("flood", 1);     // Forces wetlands/marshes
    }
    if (isHurricaneZone(region)) {
        tile.addTokens("hurricane", 1); // Specialized hurricane biome
    }
}

Climate Zone Effects on Biome Generation

Token-Based Biome Classification:

BiomeType classifyBiome(const WorldTile& tile) {
    int total_rain = tile.getAccumulatedTokens("rain");
    int total_wind = tile.getAccumulatedTokens("wind");
    int highWind_tokens = tile.getAccumulatedTokens("highWind");
    int flood_tokens = tile.getAccumulatedTokens("flood");
    int hurricane_tokens = tile.getAccumulatedTokens("hurricane");

    // Special climate zones override normal biome classification
    if (hurricane_tokens > 0) {
        return BiomeType::HURRICANE_ZONE;     // Specialized storm-resistant vegetation
    }
    if (flood_tokens > FLOOD_THRESHOLD) {
        return BiomeType::WETLANDS;           // Forced marshes/swamps
    }
    if (highWind_tokens > STORM_THRESHOLD) {
        // High wind prevents forest growth
        if (total_rain > 300) {
            return BiomeType::STORM_PRAIRIE;  // Grasslands that can handle wind
        } else {
            return BiomeType::BADLANDS;       // Sparse, wind-resistant vegetation
        }
    }

    // Normal biome classification using basic rain/wind tokens
    if (total_rain > 500) {
        return BiomeType::TROPICAL_RAINFOREST;
    } else if (total_rain < 50) {
        return BiomeType::HOT_DESERT;
    }
    // ... additional normal biome logic
}

Climate Zone Characteristics:

• Hurricane Zones → Storm-resistant palms, specialized coastal vegetation
• Flood Zones → Wetlands, marshes, swamp vegetation mandatory
• High Wind Zones → No forests allowed, prairie/badlands only
• Normal Zones → Standard biome classification by rain/temperature

Geographic Climate Patterns

Realistic Climate Formation Examples:

Congo Basin (Rainforest):

1. Large African landmass → Strong ITCZ at equator
2. Atlantic wind regions spawn → Move east via trade winds
3. ITCZ aspiration → Convergence at Congo → Amplification ×3
4. Super-humid storms → Massive rain token distribution
5. Result: Dense rainforest biome

Sahara Desert:

1. Sahara latitude (25-35°N) → Outside ITCZ band
2. No convergence zone → Wind regions pass through
3. Continental distance → Low initial moisture
4. Subtropical high pressure → Air descends (simulated via movement patterns)
5. Result: Minimal rain tokens → Desert biome

Configuration Integration

Climate Configuration (Hot-Reloadable):

{
  "climate_simulation": {
    "wind_spawn_system": {
      "base_spawn_rate": 0.1,
      "ocean_size_factor": 0.001,
      "max_concurrent_regions": 200
    },
    "planetary_circulation": {
      "trade_winds_strength": 1.5,
      "jet_stream_strength": 1.8,
      "calm_zone_chaos": 0.3
    },
    "itcz_system": {
      "latitude_band": [0.45, 0.55],
      "min_landmass_size": 10000,
      "max_ocean_distance": 800,
      "amplification_max": 3.0
    },
    "storm_thresholds": {
      "high_wind_min": 2.0,
      "flood_wetness_min": 1.5,
      "hurricane_wind_min": 2.5,
      "hurricane_rain_min": 2.0
    }
  }
}

Note: All parameters are hot-reloadable via the modular configuration system. Magic numbers are intentionally externalizable for real-time tuning during development - adjust values, save config, see immediate results without recompilation.

Performance Characteristics

Computational Complexity:

• Wind Regions: O(n) for n active regions (~50-200 simultaneously)
• ITCZ Calculations: O(m) for m convergence zones (~5-15 globally)
• Token Distribution: O(tiles_visited) per region movement
• Total per cycle: O(n × average_movement_distance)

Memory Usage:

• WindRegion: 32 bytes per region
• ITCZ Zone: 24 bytes per zone
• Token accumulation: Uses existing tile data structure
• Estimated total: <5MB for global weather simulation

Generation Time:

• Landmass analysis: 1-2 seconds (one-time setup)
• Per simulation cycle: 10-50ms for 100-200 wind regions
• Full climate stabilization: 100-500 cycles → 10-30 seconds total

Phase 7: Budget Assignment and Natural Features

Random Budget Assignment (Normal Distribution)

After climate and hydrology stabilization, assign budget scores to each tile using a bell curve distribution:

void assignBudgetScores(GMap& map) {
    std::random_device rd;
    std::mt19937 gen(rd());
    std::normal_distribution<float> budget_dist(0.0f, 3.0f); // Mean=0, σ=3

    for (int y = 0; y < map.getHeight(); y++) {
        for (int x = 0; x < map.getWidth(); x++) {
            float budget_value = budget_dist(gen);
            int8_t budget = static_cast<int8_t>(std::clamp(budget_value, -10.0f, 10.0f));

            WorldTile tile = map.getTile(x, y);
            tile.setTargetBudgetScore(budget);
        }
    }
}

Budget Distribution:

• 68% of tiles: budget score -3 to +3
• 95% of tiles: budget score -6 to +6
• Rare extremes: -10/-9 and +9/+10 scores for unique locations

Natural Features Placement (Uniform Random)

Place natural geological features randomly across the map with uniform distribution:

void placeNaturalFeatures(GMap& map, FeatureManager& feature_manager) {
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<float> prob_dist(0.0f, 1.0f);

    const float FEATURE_CHANCE = 0.05f; // 5% of tiles get features

    for (int y = 0; y < map.getHeight(); y++) {
        for (int x = 0; x < map.getWidth(); x++) {
            if (prob_dist(gen) < FEATURE_CHANCE) {
                // Natural features from gameData configuration
                std::vector<std::string> available_features = {
                    "cave", "hot_spring", "canyon", "plateau",
                    "marsh", "oasis", "geyser", "cliff", "gorge",
                    "natural_bridge", "sinkhole", "spring"
                };

                std::uniform_int_distribution<int> feature_dist(0, available_features.size() - 1);
                std::string feature = available_features[feature_dist(gen)];

                feature_manager.placeFeature(feature, x, y);
            }
        }
    }
}

Phase 8: Resource Region Conversion

Regional Influence to Resource Features

Convert abstract regional influences from geological simulation into concrete resource features on the map:

void convertRegionsToFeatures(GMap& map, FeatureManager& feature_manager) {
    std::random_device rd;
    std::mt19937 gen(rd());

    for (auto& region : map.getRegionalInfluences()) {
        // Semi-random feature count per region (1-8 features)
        std::uniform_int_distribution<int> feature_count_dist(1, 8);
        int feature_count = feature_count_dist(gen);

        for (int i = 0; i < feature_count; i++) {
            // Random position within region radius
            auto [x, y] = getRandomPositionInRegion(region);

            // Region influence strength determines feature quality
            float distance = getDistanceFromCenter(region, x, y);
            float influence_strength = region.getInfluenceAt(distance);

            // Place appropriate features based on region type
            placeRegionalFeature(region, x, y, influence_strength, feature_manager);
        }
    }
}

void placeRegionalFeature(const RegionalInfluence& region, int x, int y,
                         float strength, FeatureManager& feature_manager) {

    if (region.getType() == "tectonic_plate") {
        // Tectonic regions: metal deposits
        if (strength > 0.8f) {
            feature_manager.placeFeature("rich_iron_ore", x, y);
        } else if (strength > 0.6f) {
            feature_manager.placeFeature("copper_deposit", x, y);
        } else if (strength > 0.4f) {
            feature_manager.placeFeature("tin_deposit", x, y);
        } else {
            feature_manager.placeFeature("stone_quarry", x, y);
        }
    }

    else if (region.getType() == "carbon_region") {
        // Carbon regions: coal and oil deposits
        if (strength > 0.7f) {
            feature_manager.placeFeature("rich_coal_seam", x, y);
        } else if (strength > 0.5f) {
            feature_manager.placeFeature("oil_well", x, y);
        } else if (strength > 0.3f) {
            feature_manager.placeFeature("coal_outcrop", x, y);
        } else {
            feature_manager.placeFeature("peat_bog", x, y);
        }
    }

    else if (region.getType() == "volcanic_zone") {
        // Volcanic regions: geothermal and rare minerals
        if (strength > 0.8f) {
            feature_manager.placeFeature("geothermal_vent", x, y);
        } else if (strength > 0.6f) {
            feature_manager.placeFeature("sulfur_deposit", x, y);
        } else if (strength > 0.4f) {
            feature_manager.placeFeature("obsidian_field", x, y);
        } else {
            feature_manager.placeFeature("hot_spring", x, y);
        }
    }

    else if (region.getType() == "recent_meteorite_crater") {
        // Recent meteorite impacts: rare metals and exotic materials
        // NOTE: Different from Phase 1 planetary accretion meteorites
        if (strength > 0.9f) {
            feature_manager.placeFeature("platinum_deposit", x, y);
        } else if (strength > 0.7f) {
            feature_manager.placeFeature("rare_earth_deposit", x, y);
        } else if (strength > 0.5f) {
            feature_manager.placeFeature("gold_vein", x, y);
        } else {
            feature_manager.placeFeature("impact_glass", x, y);
        }
    }
}

std::pair<int, int> getRandomPositionInRegion(const RegionalInfluence& region) {
    std::random_device rd;
    std::mt19937 gen(rd());

    // Random angle and distance within region radius
    std::uniform_real_distribution<float> angle_dist(0.0f, 2.0f * M_PI);
    std::uniform_real_distribution<float> radius_dist(0.0f, region.getRadius());

    float angle = angle_dist(gen);
    float distance = radius_dist(gen);

    int x = static_cast<int>(region.getCenterX() + distance * cos(angle));
    int y = static_cast<int>(region.getCenterY() + distance * sin(angle));

    return {x, y};
}

Resource Feature Characteristics

Influence Strength Distribution:

• 0.8-1.0: Premium resources (rich deposits, rare materials)
• 0.6-0.8: High-quality resources (standard industrial deposits)
• 0.4-0.6: Medium-quality resources (adequate for basic industry)
• 0.2-0.4: Low-quality resources (marginal extraction)
• 0.0-0.2: Minimal resources (trace amounts only)

Regional Resource Mapping:

• Tectonic Regions: Iron, copper, tin, stone → Industrial base materials
• Carbon Regions: Coal, oil, peat → Energy resources
• Volcanic Zones: Geothermal, sulfur, obsidian → Specialized materials
• Recent Meteorite Craters: Platinum, rare earths, gold → Advanced technology materials

Feature Distribution:

• 1-8 features per region (semi-random count)
• Position: Random within region radius
• Quality: Determined by distance from region center (closer = better)
• Type: Fixed by region type, quality by influence strength

Planetary Mass Conservation System

Core-Based Mass Conservation

To maintain realistic mass conservation throughout geological simulation, the system uses a planetary core that absorbs eroded material and releases it through volcanic activity:

struct PlanetaryCore {
    float core_mass;                    // Accumulated eroded material
    float surface_mass;                 // Current surface terrain mass
    float max_core_capacity;            // Core saturation threshold
    float volcanic_overflow_rate;       // Rate of volcanic material expulsion (0.1-0.3)
    float total_planetary_mass;         // Constant after Phase 1 (meteorite bombardment)

    // Derived values
    float core_pressure_ratio;          // core_mass / max_core_capacity
    int pending_volcanic_events;        // Queued volcanic eruptions
};

void updateMassConservation(PlanetaryCore& core) {
    // Validation: Total mass conservation
    assert(core.core_mass + core.surface_mass == core.total_planetary_mass);

    // Core overflow triggers volcanic activity
    if (core.core_mass > core.max_core_capacity) {
        float overflow = core.core_mass - core.max_core_capacity;
        float volcanic_expulsion = overflow * core.volcanic_overflow_rate;

        // Transfer mass from core to pending volcanic events
        core.core_mass -= volcanic_expulsion;

        // Queue volcanic events proportional to overflow
        int new_volcanic_events = static_cast<int>(volcanic_expulsion / volcanic_event_mass_threshold);
        core.pending_volcanic_events += new_volcanic_events;

        // Store remaining material for future volcanism
        core.pending_volcanic_material += volcanic_expulsion;
    }
}

Erosion to Core Transfer

All erosion processes transfer material to the planetary core rather than redistributing on surface:

void erodeToCore(WorldTile& tile, float erosion_amount, PlanetaryCore& core) {
    // Remove material from surface
    float current_elevation = tile.getElevation();
    tile.setElevation(current_elevation - erosion_amount);

    // Transfer to planetary core
    core.core_mass += erosion_amount;
    core.surface_mass -= erosion_amount;

    // Track erosion for geological history
    tile.addFlag(TileFlags::ERODED_THIS_CYCLE);
}

// Apply to all erosion processes
void riverErosion(WorldTile& tile, float water_flow, PlanetaryCore& core) {
    if (water_flow > erosion_threshold) {
        float erosion_amount = water_flow * erosion_rate_factor;
        erodeToCore(tile, erosion_amount, core);
    }
}

void glacialErosion(WorldTile& tile, float ice_thickness, PlanetaryCore& core) {
    float erosion_amount = ice_thickness * glacial_erosion_factor;
    erodeToCore(tile, erosion_amount, core);
}

Volcanic Overflow System

Core overflow creates realistic volcanic activity that returns material to surface:

void processVolcanicOverflow(GMap& map, PlanetaryCore& core) {
    while (core.pending_volcanic_events > 0) {
        // Select volcanic location based on geological factors
        auto [x, y] = selectVolcanicLocation(map);

        // Calculate eruption magnitude
        float eruption_material = core.pending_volcanic_material / core.pending_volcanic_events;

        // Create volcanic eruption
        createVolcanicEruption(map, x, y, eruption_material);

        // Update core state
        core.surface_mass += eruption_material;
        core.pending_volcanic_material -= eruption_material;
        core.pending_volcanic_events--;
    }
}

std::pair<int, int> selectVolcanicLocation(const GMap& map) {
    // Prefer locations with:
    // - Active tectonic boundaries
    // - Existing volcanic history
    // - High core pressure influence

    std::vector<std::pair<int, int>> candidate_locations;

    for (auto& tectonic_region : map.getTectonicRegions()) {
        if (tectonic_region.getActivity() > volcanic_activity_threshold) {
            // Add boundary tiles as candidates
            auto boundary_tiles = tectonic_region.getBoundaryTiles();
            candidate_locations.insert(candidate_locations.end(),
                                     boundary_tiles.begin(), boundary_tiles.end());
        }
    }

    // Weight by distance from core pressure points
    return selectWeightedRandom(candidate_locations);
}

void createVolcanicEruption(GMap& map, int center_x, int center_y, float material_amount) {
    // Deposit material in volcanic pattern
    int eruption_radius = static_cast<int>(sqrt(material_amount / volcanic_density_factor));

    for (int dy = -eruption_radius; dy <= eruption_radius; dy++) {
        for (int dx = -eruption_radius; dx <= eruption_radius; dx++) {
            float distance = sqrt(dx*dx + dy*dy);
            if (distance <= eruption_radius) {
                int x = center_x + dx;
                int y = center_y + dy;

                if (map.isValidCoordinate(x, y)) {
                    // Volcanic deposition decreases with distance
                    float deposition_factor = 1.0f - (distance / eruption_radius);
                    float local_deposition = material_amount * deposition_factor / (M_PI * eruption_radius * eruption_radius);

                    WorldTile tile = map.getTile(x, y);
                    tile.setElevation(tile.getElevation() + local_deposition);

                    // Mark as volcanic terrain
                    tile.setFlag(TileFlags::VOLCANIC_DEPOSIT, true);
                }
            }
        }
    }
}

Mass Conservation Benefits

Perfect Conservation:

• Total planetary mass remains constant after Phase 1
• All eroded material is accounted for in core
• Volcanic activity returns material to surface
• No material created or destroyed

Realistic Geological Activity:

• Core pressure drives continuing volcanism
• Volcanic activity decreases as core pressure reduces
• Natural equilibrium between erosion and volcanic deposition
• Geological activity persists throughout simulation

Simplified Implementation:

• No complex sediment transport calculations
• No surface redistribution algorithms
• Single mass conservation equation
• Volcanic activity emerges naturally from core overflow

Gameplay Benefits:

• Ongoing geological activity creates dynamic world
• Volcanic regions provide unique resource opportunities
• Erosion/volcanism balance creates varied topography
• Long-term geological processes affect industrial planning

Technical Architecture

WorldTileData Structure (32 bytes)

struct WorldTileData {
    // Terrain (11 bytes) - Always accessed, geological simulation ready
    uint16_t terrain_type_id;       // 65k terrain types
    uint16_t biome_type_id;         // Biome classification
    uint16_t elevation;             // -11km to +32km range
    int16_t temperature;            // -3276°C to +3276°C (0.1°C precision)
    uint8_t humidity;               // 0-100% (0.4% precision)
    uint8_t wind_data;              // 4 bits direction + 4 bits intensity
    uint8_t water_level;            // Accumulated water for river formation

    // Metadata (21 bytes) - Generation and gameplay
    int8_t target_budget_score;     // -10 to +10
    uint32_t regional_influence_id; // → Regional influence data
    uint8_t influence_strength;     // 0-255
    uint32_t tile_flags;            // State flags
    uint32_t feature_set_id;        // → Feature collection
    uint8_t padding2[7];            // Future expansion space
};

RegionalInfluence System

• TectonicRegions: Mobile circular regions with physics
• CarbonRegions: Carbon deposit tracking with tectonic attachment
• VolcanicZones: Created dynamically at tectonic collisions
• SeaLevelInfluence: Global parameter affecting all processes

FeatureManager Integration

• Simple helper interface for placing geological features
• Handles feature set creation and tile updates automatically
• Used by generation algorithms, doesn't do generation itself

Performance Characteristics

Memory Usage

• 1M tiles: 32MB core tile data (increased from 24MB for hydrology)
• Feature sets: Sparse, shared between similar tiles
• Geological regions: ~10-50 regions vs millions of tiles
• Climate wind regions: ~50-200 mobile regions during simulation
• Total estimated: <125MB for complete planetary geology + advanced climate simulation

Computational Complexity

• Tectonic simulation: O(n²) for n regions (~25 regions = 625 operations)
• Meteorite impacts: O(k) for k impacts per wave
• Sea level effects: O(tiles) single pass
• Carbon processes: O(regions) sparse operations
• Climate simulation: O(n × movement_distance) for n wind regions (~50-200 regions)

Generation Time

• Geological phases: 10-20 seconds for complete 4.65 billion year simulation
• Climate simulation: 10-30 seconds for 100-500 climate cycles
• Total estimated: 20-50 seconds for complete world generation
• Progressive: Can display intermediate results during generation
• Deterministic: Same seed produces identical geology and climate

Implementation Priority

Phase 1: Core Architecture (1-2 weeks)

1. RegionalInfluence and RegionalInfluenceManager classes
2. TectonicRegion basic structure and physics
3. Simple meteorite impact system
4. Basic tile elevation/temperature updates

Phase 2: Full Geology (2-3 weeks)

1. Complete tectonic interaction system
2. Dynamic sea level integration
3. Carbon region formation and conversion
4. Volcanic zone creation

Phase 3: Polish (1-2 weeks)

1. Parameter tuning for realistic results
2. Performance optimization
3. Generation progress reporting
4. Validation and debugging tools

Scientific Accuracy vs Gameplay

Scientifically Inspired Elements

• ✅ Late Heavy Bombardment period
• ✅ Planetary differentiation (metals sink to core)
• ✅ Tectonic processes creating mountains/valleys
• ✅ Carboniferous period forest→coal formation
• ✅ Marine conditions for oil formation
• ✅ Sea level variations affecting geology
• ✅ ITCZ formation from continental heating
• ✅ Planetary circulation bands (trade winds, jet streams)
• ✅ Realistic climate differentiation (Congo vs Sahara)

Gameplay Simplifications

• ⚠️ Tectonic regions as circles (not realistic plate shapes)
• ⚠️ Fixed map size instead of planetary expansion
• ⚠️ Simplified core-mantle dynamics
• ⚠️ Compressed timescales for practical generation
• ⚠️ 2D wind simulation instead of 3D atmospheric layers
• ⚠️ Discrete token system instead of continuous climate fields
• ⚠️ Mobile wind regions as simplified weather systems

Result

Plausible geology and climate that creates interesting gameplay with emergent weather patterns without requiring PhD in atmospheric science to understand or debug.

Integration Points

WorldGenerationModule

• Orchestrates all geological phases
• Manages simulation state and progression
• Provides query interface for generated geology

FeatureManager

• Places geological features based on simulation results
• Handles feature set optimization and tile updates
• Simple interface for placement algorithms

RegionalInfluenceManager

• Core system managing all regional effects
• Handles tectonic regions, carbon regions, volcanic zones
• Provides influence application to tiles

Future Extensions

• Dynamic climate: Real-time weather systems during gameplay
• Mineral resource modeling: Detailed ore deposit formation
• Advanced erosion: River network evolution during gameplay
• Geological time compression: Speed up/slow down specific phases
• Seasonal climate: Monthly/yearly climate variations
• Climate change: Long-term climate evolution from industrial activity

---

Status: System designed and ready for implementation. All major components specified with clear interfaces and realistic performance targets. Scientific accuracy balanced with implementation complexity for practical game development.