earths 0.0.1

pub struct River {
    pub name: &'static str,
    pub length_km: f64,
    pub discharge_m3_s: f64,
    pub drainage_area_km2: f64,
    pub slope: f64,
}
impl River {
    pub fn manning_velocity(&self, hydraulic_radius_m: f64, roughness_n: f64) -> f64 {
        if roughness_n.abs() < f64::EPSILON {
            return 0.0;
        }
        (1.0 / roughness_n) * hydraulic_radius_m.powf(2.0 / 3.0) * self.slope.sqrt()
    }
    pub fn froude_number(&self, velocity_m_s: f64, depth_m: f64) -> f64 {
        if depth_m.abs() < 1e-10 {
            return 0.0;
        }
        velocity_m_s / (*crate::SURFACE_GRAVITY * depth_m.abs()).sqrt()
    }
    pub fn reynolds_number(&self, velocity_m_s: f64, depth_m: f64) -> f64 {
        let kinematic_viscosity = 1.004e-6;
        if kinematic_viscosity < 1e-15 {
            return 0.0;
        }
        velocity_m_s * depth_m / kinematic_viscosity
    }
    pub fn sediment_transport_capacity(
        &self,
        velocity_m_s: f64,
        depth_m: f64,
        grain_size_m: f64,
    ) -> f64 {
        if grain_size_m.abs() < 1e-10 {
            return 0.0;
        }
        let shields_param = SEAWATER_DENSITY_FRESH * velocity_m_s * velocity_m_s
            / ((*crate::QUARTZ_DENSITY - SEAWATER_DENSITY_FRESH)
                * *crate::SURFACE_GRAVITY
                * grain_size_m);
        shields_param * velocity_m_s * depth_m
    }
    pub fn stream_power_w_per_m(&self) -> f64 {
        SEAWATER_DENSITY_FRESH * *crate::SURFACE_GRAVITY * self.discharge_m3_s * self.slope
    }
    pub fn specific_discharge_m_s(&self) -> f64 {
        self.discharge_m3_s / (self.drainage_area_km2 * 1e6)
    }
}
const SEAWATER_DENSITY_FRESH: f64 = crate::FRESHWATER_DENSITY;
pub fn amazon() -> River {
    River {
        name: "Amazon",
        length_km: 6992.0,
        discharge_m3_s: 209_000.0,
        drainage_area_km2: 7_050_000.0,
        slope: 1.5e-5,
    }
}
pub fn nile() -> River {
    River {
        name: "Nile",
        length_km: 6650.0,
        discharge_m3_s: 2830.0,
        drainage_area_km2: 3_349_000.0,
        slope: 8.5e-5,
    }
}
pub fn chezy_velocity(chezy_c: f64, hydraulic_radius_m: f64, slope: f64) -> f64 {
    chezy_c * (hydraulic_radius_m * slope).sqrt()
}
pub fn hjulstrom_erosion_threshold(grain_diameter_mm: f64) -> f64 {
    if grain_diameter_mm < 0.5 {
        0.3 * grain_diameter_mm.powf(-0.6)
    } else {
        0.22 * grain_diameter_mm.powf(0.5)
    }
}
pub struct RiverSegment {
    pub river: River,
    pub downstream_idx: Option<usize>,
}
pub struct RiverNetwork {
    pub segments: Vec<RiverSegment>,
}
impl Default for RiverNetwork {
    fn default() -> Self {
        Self::new()
    }
}
impl RiverNetwork {
    pub fn new() -> Self {
        Self {
            segments: Vec::new(),
        }
    }
    pub fn add_headwater(&mut self, river: River) -> usize {
        let idx = self.segments.len();
        self.segments.push(RiverSegment {
            river,
            downstream_idx: None,
        });
        idx
    }
    pub fn connect(&mut self, tributary_idx: usize, main_idx: usize) {
        self.segments[tributary_idx].downstream_idx = Some(main_idx);
    }
    pub fn route_discharge(&self) -> Vec<f64> {
        let n = self.segments.len();
        let mut discharge = vec![0.0_f64; n];
        for (i, seg) in self.segments.iter().enumerate() {
            discharge[i] = seg.river.discharge_m3_s;
        }
        let mut in_degree = vec![0u32; n];
        for seg in &self.segments {
            if let Some(ds) = seg.downstream_idx {
                in_degree[ds] += 1;
            }
        }
        let mut queue: Vec<usize> = (0..n).filter(|&i| in_degree[i] == 0).collect();
        let mut order = Vec::with_capacity(n);
        let mut remaining_in = in_degree.clone();
        while let Some(idx) = queue.pop() {
            order.push(idx);
            if let Some(ds) = self.segments[idx].downstream_idx {
                remaining_in[ds] -= 1;
                if remaining_in[ds] == 0 {
                    queue.push(ds);
                }
            }
        }
        for &idx in &order {
            if let Some(ds) = self.segments[idx].downstream_idx {
                discharge[ds] += discharge[idx];
            }
        }
        discharge
    }
    pub fn outlet_discharge(&self) -> f64 {
        let discharge = self.route_discharge();
        self.segments
            .iter()
            .enumerate()
            .filter(|(_, seg)| seg.downstream_idx.is_none())
            .map(|(i, _)| discharge[i])
            .fold(0.0, f64::max)
    }
}
pub fn amazon_basin() -> RiverNetwork {
    let mut net = RiverNetwork::new();
    let solimoes = net.add_headwater(River {
        name: "Solimões",
        length_km: 1600.0,
        discharge_m3_s: 103_000.0,
        drainage_area_km2: 2_200_000.0,
        slope: 2.0e-5,
    });
    let negro = net.add_headwater(River {
        name: "Rio Negro",
        length_km: 2250.0,
        discharge_m3_s: 28_400.0,
        drainage_area_km2: 700_000.0,
        slope: 3.0e-5,
    });
    let madeira = net.add_headwater(River {
        name: "Madeira",
        length_km: 3380.0,
        discharge_m3_s: 31_200.0,
        drainage_area_km2: 1_380_000.0,
        slope: 4.0e-5,
    });
    let amazon_main = net.add_headwater(River {
        name: "Amazon (main)",
        length_km: 1400.0,
        discharge_m3_s: 20_000.0,
        drainage_area_km2: 500_000.0,
        slope: 1.5e-5,
    });
    net.connect(solimoes, amazon_main);
    net.connect(negro, amazon_main);
    net.connect(madeira, amazon_main);
    net
}