pub struct River {
pub name: &'static str,
pub lengthkm: f64,
pub dischargem3s: f64,
pub drainageareakm2: f64,
pub slope: f64,
}
impl River {
pub fn manningvelocity(&self, hydraulicradiusm: f64, roughnessn: f64) -> f64 {
if roughnessn.abs() < f64::EPSILON {
return 0.0;
}
(1.0 / roughnessn) * hydraulicradiusm.powf(2.0 / 3.0) * self.slope.sqrt()
}
pub fn froudenumber(&self, velocityms: f64, depthm: f64) -> f64 {
if depthm.abs() < 1e-10 {
return 0.0;
}
velocityms / (*crate::SURFACEGRAVITY * depthm.abs()).sqrt()
}
pub fn reynoldsnumber(&self, velocityms: f64, depthm: f64) -> f64 {
let kinematicviscosity = 1.004e-6;
if kinematicviscosity < 1e-15 {
return 0.0;
}
velocityms * depthm / kinematicviscosity
}
pub fn sedimenttransportcapacity(&self, velocityms: f64, depthm: f64, grainsizem: f64) -> f64 {
if grainsizem.abs() < 1e-10 {
return 0.0;
}
let shieldsparam = SEAWATERDENSITYFRESH * velocityms * velocityms
/ ((*crate::QUARTZDENSITY - SEAWATERDENSITYFRESH)
* *crate::SURFACEGRAVITY
* grainsizem);
shieldsparam * velocityms * depthm
}
pub fn streampowerwperm(&self) -> f64 {
SEAWATERDENSITYFRESH * *crate::SURFACEGRAVITY * self.dischargem3s * self.slope
}
pub fn specificdischargems(&self) -> f64 {
self.dischargem3s / (self.drainageareakm2 * 1e6)
}
}
const SEAWATERDENSITYFRESH: f64 = crate::FRESHWATERDENSITY;
pub fn amazon() -> River {
River {
name: "Amazon",
lengthkm: 6992.0,
dischargem3s: 209000.0,
drainageareakm2: 7050000.0,
slope: 1.5e-5,
}
}
pub fn nile() -> River {
River {
name: "Nile",
lengthkm: 6650.0,
dischargem3s: 2830.0,
drainageareakm2: 3349000.0,
slope: 8.5e-5,
}
}
pub fn chezyvelocity(chezyc: f64, hydraulicradiusm: f64, slope: f64) -> f64 {
chezyc * (hydraulicradiusm * slope).sqrt()
}
pub fn hjulstromerosionthreshold(graindiametermm: f64) -> f64 {
if graindiametermm < 0.5 {
0.3 * graindiametermm.powf(-0.6)
} else {
0.22 * graindiametermm.powf(0.5)
}
}
pub struct RiverSegment {
pub river: River,
pub downstreamidx: 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 addheadwater(&mut self, river: River) -> usize {
let idx = self.segments.len();
self.segments.push(RiverSegment {
river,
downstreamidx: None,
});
idx
}
pub fn connect(&mut self, tributaryidx: usize, mainidx: usize) {
self.segments[tributaryidx].downstreamidx = Some(mainidx);
}
pub fn routedischarge(&self) -> Vec<f64> {
let n = self.segments.len();
let mut discharge = vec![0.0f64; n];
for (i, seg) in self.segments.iter().enumerate() {
discharge[i] = seg.river.dischargem3s;
}
let mut indegree = vec![0u32; n];
for seg in &self.segments {
if let Some(ds) = seg.downstreamidx {
indegree[ds] += 1;
}
}
let mut queue: Vec<usize> = (0..n).filter(|&i| indegree[i] == 0).collect();
let mut order = Vec::with_capacity(n);
let mut remainingin = indegree.clone();
while let Some(idx) = queue.pop() {
order.push(idx);
if let Some(ds) = self.segments[idx].downstreamidx {
remainingin[ds] -= 1;
if remainingin[ds] == 0 {
queue.push(ds);
}
}
}
for &idx in &order {
if let Some(ds) = self.segments[idx].downstreamidx {
discharge[ds] += discharge[idx];
}
}
discharge
}
pub fn outletdischarge(&self) -> f64 {
let discharge = self.routedischarge();
self.segments
.iter()
.enumerate()
.filter(|entry| entry.1.downstreamidx.is_none())
.map(|entry| discharge[entry.0])
.fold(0.0, f64::max)
}
}
pub fn amazonbasin() -> RiverNetwork {
let mut net = RiverNetwork::new();
let solimoes = net.addheadwater(River {
name: "Solimões",
lengthkm: 1600.0,
dischargem3s: 103000.0,
drainageareakm2: 2200000.0,
slope: 2.0e-5,
});
let negro = net.addheadwater(River {
name: "Rio Negro",
lengthkm: 2250.0,
dischargem3s: 28400.0,
drainageareakm2: 700000.0,
slope: 3.0e-5,
});
let madeira = net.addheadwater(River {
name: "Madeira",
lengthkm: 3380.0,
dischargem3s: 31200.0,
drainageareakm2: 1380000.0,
slope: 4.0e-5,
});
let amazonmain = net.addheadwater(River {
name: "Amazon (main)",
lengthkm: 1400.0,
dischargem3s: 20000.0,
drainageareakm2: 500000.0,
slope: 1.5e-5,
});
net.connect(solimoes, amazonmain);
net.connect(negro, amazonmain);
net.connect(madeira, amazonmain);
net
}