use sciforge::hub::domain::geology::seismology::{
    epicenter_distance, gutenberg_richter, moment_magnitude, omori_aftershock, p_wave_velocity,
    peak_ground_acceleration, s_wave_velocity, seismic_energy, seismic_moment, snell_seismic,
    travel_time,
};

pub struct Icequake {
    pub latitude: f64,
    pub longitude: f64,
    pub depthkm: f64,
    pub seismicmomentnm: f64,
}

impl Icequake {
    pub fn frommoment(lat: f64, lon: f64, depthkm: f64, momentnm: f64) -> Self {
        Self {
            latitude: lat,
            longitude: lon,
            depthkm,
            seismicmomentnm: momentnm,
        }
    }

    pub fn frommagnitude(lat: f64, lon: f64, depthkm: f64, mw: f64) -> Self {
        let moment = seismic_moment(mw);
        Self::frommoment(lat, lon, depthkm, moment)
    }

    pub fn momentmagnitude(&self) -> f64 {
        moment_magnitude(self.seismicmomentnm)
    }

    pub fn energyjoules(&self) -> f64 {
        seismic_energy(self.momentmagnitude())
    }

    pub fn pwavearrivaltime(&self, distancekm: f64, vp: f64) -> f64 {
        travel_time(distancekm * 1000.0, vp)
    }

    pub fn swavearrivaltime(&self, distancekm: f64, vs: f64) -> f64 {
        travel_time(distancekm * 1000.0, vs)
    }

    pub fn pgaatdistance(&self, distancekm: f64) -> f64 {
        peak_ground_acceleration(1.0, 1.0, self.momentmagnitude(), distancekm)
    }
}

pub fn aftershockrate(mainshockmagnitude: f64, timedays: f64) -> f64 {
    let a = mainshockmagnitude - 1.2;
    let c = 0.05;
    let p = 1.1;
    omori_aftershock(a.powf(10.0), c, p, timedays)
}

pub fn annualfrequency(magnitude: f64) -> f64 {
    let a = 4.0;
    let b = 0.8;
    gutenberg_richter(a, b, magnitude)
}

pub fn icepwavevelocity(density: f64, bulkmodulus: f64, shearmodulus: f64) -> f64 {
    p_wave_velocity(bulkmodulus, shearmodulus, density)
}

pub fn iceswavevelocity(density: f64, shearmodulus: f64) -> f64 {
    s_wave_velocity(shearmodulus, density)
}

pub fn distancefromsptime(tp: f64, ts: f64, vp: f64, vs: f64) -> f64 {
    epicenter_distance(vp, vs, ts - tp)
}

pub fn refractionangle(v1: f64, v2: f64, angle1rad: f64) -> f64 {
    snell_seismic(v1, angle1rad, v2)
}
