use serde::{Deserialize, Serialize};
#[derive(Debug, Clone, Serialize, Deserialize, PartialEq)]
pub struct GeologicUnit {
pub name: String,
pub rock_type: String,
pub age_ma: f64,
}
#[derive(Debug, Clone, Copy, Serialize, Deserialize, PartialEq)]
pub struct StrikeDip {
pub strike_deg: f64,
pub dip_deg: f64,
}
impl StrikeDip {
#[must_use]
pub fn dip_direction(&self) -> f64 {
(self.strike_deg + 90.0) % 360.0
}
}
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct GeologicGrid {
pub nx: usize,
pub ny: usize,
pub cell_size_m: f64,
cells: Vec<Option<GeologicUnit>>,
}
impl GeologicGrid {
#[must_use]
pub fn new(nx: usize, ny: usize, cell_size_m: f64) -> Self {
Self {
nx,
ny,
cell_size_m,
cells: vec![None; nx * ny],
}
}
pub fn set(&mut self, x: usize, y: usize, unit: GeologicUnit) {
assert!(x < self.nx, "x index {x} out of bounds (nx={})", self.nx);
assert!(y < self.ny, "y index {y} out of bounds (ny={})", self.ny);
self.cells[y * self.nx + x] = Some(unit);
}
#[must_use]
pub fn get(&self, x: usize, y: usize) -> Option<&GeologicUnit> {
assert!(x < self.nx, "x index {x} out of bounds (nx={})", self.nx);
assert!(y < self.ny, "y index {y} out of bounds (ny={})", self.ny);
self.cells[y * self.nx + x].as_ref()
}
pub fn fill(&mut self, unit: GeologicUnit) {
for cell in &mut self.cells {
*cell = Some(unit.clone());
}
}
}
#[derive(Debug, Clone, Serialize, Deserialize, PartialEq)]
pub struct GeologicLayer {
pub name: String,
pub thickness_m: f64,
pub rock_type: String,
pub depth_top_m: f64,
}
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct StratigraphicColumn {
layers: Vec<GeologicLayer>,
}
impl StratigraphicColumn {
#[must_use]
pub fn new() -> Self {
Self { layers: Vec::new() }
}
pub fn add_layer(&mut self, name: String, thickness_m: f64, rock_type: String) {
let depth_top_m = self.total_thickness();
self.layers.push(GeologicLayer {
name,
thickness_m,
rock_type,
depth_top_m,
});
}
#[must_use]
pub fn total_thickness(&self) -> f64 {
self.layers.iter().map(|l| l.thickness_m).sum()
}
#[must_use]
pub fn layer_at_depth(&self, depth_m: f64) -> Option<&GeologicLayer> {
self.layers
.iter()
.find(|l| depth_m >= l.depth_top_m && depth_m < l.depth_top_m + l.thickness_m)
}
}
impl Default for StratigraphicColumn {
fn default() -> Self {
Self::new()
}
}
#[cfg(test)]
mod tests {
use super::*;
fn sample_unit(name: &str) -> GeologicUnit {
GeologicUnit {
name: name.to_string(),
rock_type: "granite".to_string(),
age_ma: 250.0,
}
}
#[test]
fn grid_creation() {
let grid = GeologicGrid::new(4, 3, 10.0);
assert_eq!(grid.nx, 4);
assert_eq!(grid.ny, 3);
assert_eq!(grid.cell_size_m, 10.0);
assert!(grid.get(0, 0).is_none());
}
#[test]
fn grid_set_and_get() {
let mut grid = GeologicGrid::new(5, 5, 1.0);
let unit = sample_unit("Basalt Flow");
grid.set(2, 3, unit.clone());
let retrieved = grid.get(2, 3).unwrap();
assert_eq!(retrieved, &unit);
}
#[test]
fn grid_get_empty_cell() {
let grid = GeologicGrid::new(3, 3, 5.0);
assert!(grid.get(1, 1).is_none());
}
#[test]
fn grid_fill() {
let mut grid = GeologicGrid::new(2, 2, 1.0);
let unit = sample_unit("Sandstone");
grid.fill(unit.clone());
for y in 0..2 {
for x in 0..2 {
assert_eq!(grid.get(x, y).unwrap(), &unit);
}
}
}
#[test]
#[should_panic(expected = "out of bounds")]
fn grid_set_out_of_bounds() {
let mut grid = GeologicGrid::new(3, 3, 1.0);
grid.set(5, 0, sample_unit("oob"));
}
#[test]
#[should_panic(expected = "out of bounds")]
fn grid_get_out_of_bounds() {
let grid = GeologicGrid::new(3, 3, 1.0);
let _ = grid.get(0, 10);
}
#[test]
fn strike_dip_direction() {
let sd = StrikeDip {
strike_deg: 45.0,
dip_deg: 30.0,
};
assert!((sd.dip_direction() - 135.0).abs() < f64::EPSILON);
}
#[test]
fn strike_dip_wraps_around() {
let sd = StrikeDip {
strike_deg: 300.0,
dip_deg: 60.0,
};
assert!((sd.dip_direction() - 30.0).abs() < f64::EPSILON);
}
#[test]
fn column_construction_and_thickness() {
let mut col = StratigraphicColumn::new();
col.add_layer("Topsoil".into(), 2.0, "soil".into());
col.add_layer("Limestone".into(), 15.0, "carbonate".into());
col.add_layer("Shale".into(), 8.0, "mudrock".into());
assert!((col.total_thickness() - 25.0).abs() < f64::EPSILON);
assert_eq!(col.layers.len(), 3);
}
#[test]
fn column_layer_at_depth() {
let mut col = StratigraphicColumn::new();
col.add_layer("Topsoil".into(), 2.0, "soil".into());
col.add_layer("Limestone".into(), 15.0, "carbonate".into());
col.add_layer("Shale".into(), 8.0, "mudrock".into());
let layer = col.layer_at_depth(0.5).unwrap();
assert_eq!(layer.name, "Topsoil");
let layer = col.layer_at_depth(10.0).unwrap();
assert_eq!(layer.name, "Limestone");
let layer = col.layer_at_depth(20.0).unwrap();
assert_eq!(layer.name, "Shale");
assert!(col.layer_at_depth(25.0).is_none());
}
#[test]
fn column_default() {
let col = StratigraphicColumn::default();
assert!((col.total_thickness() - 0.0).abs() < f64::EPSILON);
}
#[test]
fn geologic_unit_serde_roundtrip() {
let unit = sample_unit("Gneiss");
let json = serde_json::to_string(&unit).unwrap();
let deserialized: GeologicUnit = serde_json::from_str(&json).unwrap();
assert_eq!(unit, deserialized);
}
}