pub struct Aeroplankton {
    pub speciesname: &'static str,
    pub count: f64,
    pub carryingcapacity: f64,
    pub intrinsicgrowthrate: f64,
    pub effectivemasskg: f64,
}

impl Aeroplankton {
    pub fn growthrate(&self) -> f64 {
        self.intrinsicgrowthrate
            * self.count
            * (1.0 - self.count / self.carryingcapacity.max(1e-30))
    }

    pub fn projectforward(&self, dtyears: f64) -> f64 {
        let r = self.intrinsicgrowthrate;
        let k = self.carryingcapacity;
        let n = self.count;
        if n.abs() < 1e-30 || k.abs() < 1e-30 {
            return n;
        }
        k / (1.0 + ((k - n) / n) * (-r * dtyears).exp())
    }

    pub fn metabolicratew(&self) -> f64 {
        3.5 * self.effectivemasskg.powf(0.75)
    }

    pub fn doublingtimeyr(&self) -> f64 {
        if self.intrinsicgrowthrate <= 0.0 {
            return f64::INFINITY;
        }
        (2.0f64).ln() / self.intrinsicgrowthrate
    }

    pub fn dragcoefficient(reynoldsnumber: f64) -> f64 {
        if reynoldsnumber < 1.0 {
            24.0 / reynoldsnumber.max(1e-30)
        } else if reynoldsnumber < 1000.0 {
            24.0 / reynoldsnumber * (1.0 + 0.15 * reynoldsnumber.powf(0.687))
        } else {
            0.44
        }
    }

    pub fn terminalvelocityms(&self, cellradiusm: f64, gasdensitykgm3: f64) -> f64 {
        let cd = 0.44;
        let area = std::f64::consts::PI * cellradiusm * cellradiusm;
        let denom = 0.5 * cd * gasdensitykgm3 * area;
        if denom < 1e-30 {
            return 0.0;
        }
        (self.effectivemasskg * 24.79 / denom).sqrt()
    }
}

pub struct CloudLayerInteraction {
    pub producers: Aeroplankton,
    pub consumers: Aeroplankton,
    pub interactionrate: f64,
    pub conversionefficiency: f64,
    pub consumerlossrate: f64,
}

impl CloudLayerInteraction {
    pub fn producergrowthrate(&self) -> f64 {
        let logistic = self.producers.growthrate();
        let predation = self.interactionrate * self.producers.count * self.consumers.count;
        logistic - predation
    }

    pub fn consumergrowthrate(&self) -> f64 {
        let gains = self.conversionefficiency
            * self.interactionrate
            * self.producers.count
            * self.consumers.count;
        let losses = self.consumerlossrate * self.consumers.count;
        gains - losses
    }

    pub fn step(&mut self, dt: f64) {
        let (n, p) = (self.producers.count, self.consumers.count);
        let deriv = |prod: f64, cons: f64| -> (f64, f64) {
            let prodlogistic = self.producers.intrinsicgrowthrate
                * prod
                * (1.0 - prod / self.producers.carryingcapacity.max(1e-30));
            let predation = self.interactionrate * prod * cons;
            let dn = prodlogistic - predation;
            let dp = self.conversionefficiency * self.interactionrate * prod * cons
                - self.consumerlossrate * cons;
            (dn, dp)
        };
        let (k1n, k1p) = deriv(n, p);
        let (k2n, k2p) = deriv(n + 0.5 * dt * k1n, p + 0.5 * dt * k1p);
        let (k3n, k3p) = deriv(n + 0.5 * dt * k2n, p + 0.5 * dt * k2p);
        let (k4n, k4p) = deriv(n + dt * k3n, p + dt * k3p);
        self.producers.count = (n + dt / 6.0 * (k1n + 2.0 * k2n + 2.0 * k3n + k4n)).max(0.0);
        self.consumers.count = (p + dt / 6.0 * (k1p + 2.0 * k2p + 2.0 * k3p + k4p)).max(0.0);
    }
}

pub fn hypotheticalh2floater() -> Aeroplankton {
    Aeroplankton {
        speciesname: "Hypothetical H2 Floater",
        count: 0.0,
        carryingcapacity: 0.0,
        intrinsicgrowthrate: 0.0,
        effectivemasskg: 1.0,
    }
}

pub fn hypotheticalchemoautotroph() -> Aeroplankton {
    Aeroplankton {
        speciesname: "Hypothetical Chemoautotroph",
        count: 0.0,
        carryingcapacity: 0.0,
        intrinsicgrowthrate: 0.0,
        effectivemasskg: 1e-12,
    }
}
