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

pub fn icedensity() -> f64 {
    *crate::WATERICEDENSITY
}

pub const ICECREEPEXPONENT: f64 = 3.0;
pub const GLENFLOWLAWA: f64 = 2.4e-24;

pub fn watermolarmassfromelements() -> f64 {
    crate::mh2o()
}

pub struct IceLayer {
    pub name: &'static str,
    pub thicknessm: f64,
    pub surfaceslopedeg: f64,
    pub accumulationmyr: f64,
    pub ablationmyr: f64,
}

impl IceLayer {
    pub fn massbalancemyr(&self) -> f64 {
        self.accumulationmyr - self.ablationmyr
    }

    pub fn basalshearstresspa(&self) -> f64 {
        let g = *crate::SURFACEGRAVITY;
        let sloperad = self.surfaceslopedeg.to_radians();
        icedensity() * g * self.thicknessm * sloperad.sin()
    }

    pub fn deformationvelocitymyr(&self) -> f64 {
        let tau = self.basalshearstresspa();
        let n = ICECREEPEXPONENT;
        let velocityms = 2.0 * GLENFLOWLAWA * tau.powf(n) * self.thicknessm / (n + 1.0);
        velocityms * crate::SECONDSPERYEAR
    }

    pub fn volumem3(&self, areakm2: f64) -> f64 {
        areakm2 * 1e6 * self.thicknessm * 0.7
    }
}

pub fn watericemantle() -> IceLayer {
    IceLayer {
        name: "Water Ice Mantle",
        thicknessm: 7000000.0,
        surfaceslopedeg: 0.001,
        accumulationmyr: 0.0,
        ablationmyr: 0.0,
    }
}

pub fn methaneicecap() -> IceLayer {
    IceLayer {
        name: "Methane Ice Cap",
        thicknessm: 1000.0,
        surfaceslopedeg: 0.01,
        accumulationmyr: 0.001,
        ablationmyr: 0.0008,
    }
}

pub fn icerheologyviscosity(temperaturek: f64, stresspa: f64) -> f64 {
    let q = 60000.0;
    let a = GLENFLOWLAWA * (-q / (R_GAS * temperaturek)).exp();
    1.0 / (2.0 * a * stresspa.powf(ICECREEPEXPONENT - 1.0))
}

pub fn glenstrainrate(temperaturek: f64, stresspa: f64) -> f64 {
    let q = 60000.0;
    let atemp = GLENFLOWLAWA * (-q / (R_GAS * temperaturek)).exp();
    atemp * stresspa.powf(ICECREEPEXPONENT)
}

pub fn shallowicevelocity(thicknessm: f64, sloperad: f64, temperaturek: f64) -> f64 {
    let q = 60000.0;
    let atemp = GLENFLOWLAWA * (-q / (R_GAS * temperaturek)).exp();
    let g = *crate::SURFACEGRAVITY;
    let n = ICECREEPEXPONENT;
    let drivingstress = icedensity() * g * sloperad.sin().abs();
    2.0 * atemp / (n + 1.0) * drivingstress.powf(n) * thicknessm.powf(n + 1.0)
}
