use sciforge::hub::domain::common::constants::SIGMA_SB;

pub struct ThermalCycle {
    pub subsolar_temp_k: f64,
    pub nightside_temp_k: f64,
    pub thermal_inertia: f64,
    pub latitude_deg: f64,
}

impl ThermalCycle {
    pub fn equatorial() -> Self {
        Self {
            subsolar_temp_k: crate::SUBSOLAR_TEMP_K,
            nightside_temp_k: crate::NIGHTSIDE_TEMP_K,
            thermal_inertia: crate::THERMAL_INERTIA,
            latitude_deg: 0.0,
        }
    }

    pub fn polar() -> Self {
        Self {
            subsolar_temp_k: 400.0,
            nightside_temp_k: crate::POLAR_COLD_TRAP_TEMP_K,
            thermal_inertia: crate::THERMAL_INERTIA,
            latitude_deg: 85.0,
        }
    }

    pub fn temperature_swing_k(&self) -> f64 {
        self.subsolar_temp_k - self.nightside_temp_k
    }

    pub fn temperature_at_phase(&self, phase: f64) -> f64 {
        let mean = (self.subsolar_temp_k + self.nightside_temp_k) / 2.0;
        let amplitude = (self.subsolar_temp_k - self.nightside_temp_k) / 2.0;
        let lag = 0.02 * self.thermal_inertia / 75.0;
        let angle = 2.0 * std::f64::consts::PI * (phase - lag);
        mean + amplitude * (-angle).cos()
    }

    pub fn cooling_rate_at_phase(&self, phase: f64) -> f64 {
        let dt = 60.0;
        let t1 = self.temperature_at_phase(phase);
        let t2 = self.temperature_at_phase(phase + dt / crate::SOLAR_DAY_S);
        (t2 - t1) / dt
    }

    pub fn thermal_stress_mpa(&self) -> f64 {
        let youngs_modulus = 7e10;
        let thermal_expansion = 8e-6;
        youngs_modulus * thermal_expansion * self.temperature_swing_k() / 1e6
    }

    pub fn peak_emission(&self) -> f64 {
        crate::SURFACE_EMISSIVITY * SIGMA_SB * self.subsolar_temp_k.powi(4)
    }

    pub fn thermal_skin_depth(&self) -> f64 {
        let period_s = crate::SOLAR_DAY_S;
        let kappa = self.thermal_inertia;
        let rho = 1500.0;
        let cp = 800.0;
        let diffusivity = (kappa / (rho * cp)).powi(2);
        (diffusivity * period_s / std::f64::consts::PI).sqrt()
    }
}

pub fn hot_pole_longitudes_deg() -> (f64, f64) {
    (0.0, 180.0)
}

pub fn warm_pole_longitudes_deg() -> (f64, f64) {
    (90.0, 270.0)
}

pub fn hot_pole_max_temp_k() -> f64 {
    crate::SUBSOLAR_TEMP_K
}

pub fn warm_pole_max_temp_k() -> f64 {
    let flux_ratio = (0.307_f64 / 0.467).powi(2);
    (crate::SUBSOLAR_TEMP_K.powi(4) * flux_ratio).powf(0.25)
}
