using System;

public static class PolarIceCapModel
{
    // Current masses of the ice caps, updated dynamically
    public static double massAntarctica_kg { get; private set; } = EarthConstantsAndCalculations.ANTARCTICA_ICE_MASS_KG;
    public static double massGreenland_kg { get; private set; } = EarthConstantsAndCalculations.GREENLAND_ICE_MASS_KG;

    /// <summary>
    /// Updates the mass of each ice cap based on a 6-month seasonal cycle.
    /// Antarctica melts during Southern summer, Greenland melts during Northern summer.
    /// </summary>
    /// <param name="simulationTime_s">Total elapsed simulation time in seconds.</param>
    public static void UpdateSeasonalMass(double simulationTime_s)
    {
        const double SECONDS_PER_YEAR = 365.25 * 24 * 3600;
        double phase_rad = (2.0 * Math.PI * simulationTime_s) / SECONDS_PER_YEAR;

        // Southern Hemisphere summer (Antarctic melt) peaks around day 0/365.
        // We use -cos(phase) so that at t=0, the value is -1 (max melt).
        double antarcticMassChange = -EarthConstantsAndCalculations.SEASONAL_ICE_MASS_SWING_KG * Math.Cos(phase_rad);
        
        // Northern Hemisphere summer (Greenland melt) is 6 months out of phase.
        // We use +cos(phase) so that at t=0, the value is +1 (max accumulation), and at 6 months, -1 (max melt).
        // To offset by 6 months (half a year), we can add Math.PI to the phase.
        double greenlandMassChange = EarthConstantsAndCalculations.SEASONAL_ICE_MASS_SWING_KG * Math.Cos(phase_rad + Math.PI);

        massAntarctica_kg = EarthConstantsAndCalculations.ANTARCTICA_ICE_MASS_KG + antarcticMassChange;
        massGreenland_kg = EarthConstantsAndCalculations.GREENLAND_ICE_MASS_KG + greenlandMassChange;
    }

    /// <summary>
    /// Calculates the total moment of inertia from both polar ice caps.
    /// Treats each cap as a point mass at its average latitude.
    /// MOI for a point mass = m * r^2, where r is the perpendicular distance to the axis of rotation.
    /// r = R_Earth * cos(latitude).
    /// </summary>
    /// <returns>Total moment of inertia from ice in kg·m^2.</returns>
    public static double GetTotalIceMomentOfInertia()
    {
        // MOI for Antarctica
        double latAntarctica_rad = EarthConstantsAndCalculations.ANTARCTICA_AVG_LATITUDE_DEG * Math.PI / 180.0;
        double r_antarctica = EarthConstantsAndCalculations.WGS84_EQUATORIAL_RADIUS_M * Math.Cos(latAntarctica_rad);
        double moi_antarctica = massAntarctica_kg * r_antarctica * r_antarctica;

        // MOI for Greenland
        double latGreenland_rad = EarthConstantsAndCalculations.GREENLAND_AVG_LATITUDE_DEG * Math.PI / 180.0;
        double r_greenland = EarthConstantsAndCalculations.WGS84_EQUATORIAL_RADIUS_M * Math.Cos(latGreenland_rad);
        double moi_greenland = massGreenland_kg * r_greenland * r_greenland;

        return moi_antarctica + moi_greenland;
    }
}

Planetary Instability Model PIM - Copyright (C) 2025 James Pacha