use crate::physics::orbit::MERCURY_RADIUS;
fn shortest_lon_diff(a: f64, b: f64) -> f64 {
let mut d = a - b;
if d > 180.0 {
d -= 360.0;
}
if d < -180.0 {
d += 360.0;
}
d
}
fn basin_depth(lat: f64, lon: f64) -> f64 {
const BASINS: &[(f64, f64, f64, f64)] = &[
(30.5, 170.0, 12.0, -2000.0),
(-33.0, -67.0, 5.5, -1500.0),
(27.0, -57.0, 2.5, -1200.0),
(-20.0, -164.0, 5.0, -1500.0),
(-16.0, -165.0, 4.0, -1000.0),
(-4.0, -151.0, 3.5, -900.0),
];
let mut depth = 0.0_f64;
for &(clat, clon, radius_deg, max_d) in BASINS {
let dlat = (lat - clat) / radius_deg;
let dlon = shortest_lon_diff(lon, clon) / radius_deg;
let d2 = dlat * dlat + dlon * dlon;
if d2 < 4.0 {
let basin_profile = max_d * (-d2 * 0.8).exp();
let rim = -max_d * 0.3 * (-((d2.sqrt() - 1.0).powi(2)) * 3.0).exp();
depth = depth.min(basin_profile + rim);
}
}
depth
}
fn crater_rim_height(lat: f64, lon: f64) -> f64 {
const CRATERS: &[(f64, f64, f64, f64)] = &[
(58.4, -30.0, 1.5, 1200.0),
(-51.0, 180.0, 1.0, 800.0),
(12.0, 137.0, 2.0, 900.0),
(44.0, -134.0, 1.5, 700.0),
(-3.0, -14.0, 0.8, 600.0),
(7.0, -167.0, 1.0, 1000.0),
];
let mut elev = 0.0_f64;
for &(clat, clon, radius_deg, rim_h) in CRATERS {
let dlat = (lat - clat) / radius_deg;
let dlon = shortest_lon_diff(lon, clon) / radius_deg;
let d2 = dlat * dlat + dlon * dlon;
let rim = rim_h * (-((d2.sqrt() - 1.0).powi(2)) * 4.0).exp();
elev = elev.max(rim);
}
elev
}
fn scarp_height(lat: f64, lon: f64) -> f64 {
const SCARPS: &[(f64, f64, f64, f64, f64, f64)] = &[
(-53.0, -38.0, 5.0, 0.5, 0.3, 2000.0),
(-3.0, -100.0, 4.5, 0.5, 1.2, 1000.0),
(26.0, -65.0, 6.0, 0.6, 0.8, 3000.0),
(50.0, -32.0, 3.0, 0.4, 0.5, 1500.0),
];
let mut h = 0.0_f64;
for &(clat, clon, len, width, az, height) in SCARPS {
let dlat = lat - clat;
let dlon = shortest_lon_diff(lon, clon);
let along = dlat * az.cos() + dlon * az.sin();
let perp = -dlat * az.sin() + dlon * az.cos();
if along.abs() < len / 2.0 {
let scarp = height * (-(perp / width).powi(2)).exp();
h = h.max(scarp);
}
}
h
}
fn smooth_plains_offset(lat: f64, lon: f64) -> f64 {
const PLAINS: &[(f64, f64, f64, f64)] = &[
(65.0, 30.0, 20.0, -500.0),
(70.0, -30.0, 15.0, -400.0),
(30.5, 170.0, 8.0, -300.0),
];
let mut offset = 0.0_f64;
for &(clat, clon, radius_deg, depth) in PLAINS {
let dlat = (lat - clat) / radius_deg;
let dlon = shortest_lon_diff(lon, clon) / radius_deg;
let d2 = dlat * dlat + dlon * dlon;
let contribution = depth * (-d2 * 0.5).exp();
offset = offset.min(contribution);
}
offset
}
pub fn mercury_elevation(lat: f64, lon: f64) -> f64 {
let basin = basin_depth(lat, lon);
let crater = crater_rim_height(lat, lon);
let scarp = scarp_height(lat, lon);
let plains = smooth_plains_offset(lat, lon);
let base = 0.0;
let elev = base + basin + crater + scarp + plains;
elev.clamp(crate::LOWEST_POINT_M, crate::HIGHEST_POINT_M)
}
pub struct Heightmap {
pub resolution: u32,
pub data: Vec<f64>,
pub min_elevation_m: f64,
pub max_elevation_m: f64,
}
impl Heightmap {
pub fn new(resolution: u32) -> Self {
let size = (resolution * resolution) as usize;
Self {
resolution,
data: vec![0.0; size],
min_elevation_m: crate::LOWEST_POINT_M,
max_elevation_m: crate::HIGHEST_POINT_M,
}
}
pub fn sample(&self, lat_deg: f64, lon_deg: f64) -> f64 {
let u = ((lon_deg + 180.0) / 360.0).clamp(0.0, 1.0);
let v = ((lat_deg + 90.0) / 180.0).clamp(0.0, 1.0);
let x = (u * (self.resolution - 1) as f64) as usize;
let y = (v * (self.resolution - 1) as f64) as usize;
let idx = y * self.resolution as usize + x;
self.data.get(idx).copied().unwrap_or(0.0)
}
pub fn radius_at(&self, lat_deg: f64, lon_deg: f64) -> f64 {
MERCURY_RADIUS + self.sample(lat_deg, lon_deg)
}
pub fn generate_mercury_topography(&mut self) {
for y in 0..self.resolution {
for x in 0..self.resolution {
let lon = (x as f64 / (self.resolution - 1) as f64) * 360.0 - 180.0;
let lat = (y as f64 / (self.resolution - 1) as f64) * 180.0 - 90.0;
let idx = (y * self.resolution + x) as usize;
self.data[idx] = mercury_elevation(lat, lon);
}
}
}
pub fn generate_procedural(
&mut self,
seed: u64,
octaves: u32,
persistence: f64,
lacunarity: f64,
) {
let mut rng_state = seed;
for y in 0..self.resolution {
for x in 0..self.resolution {
let mut amplitude = 1.0;
let mut frequency = 1.0;
let mut value = 0.0;
for _ in 0..octaves {
rng_state = rng_state
.wrapping_mul(6364136223846793005)
.wrapping_add(frequency as u64);
let noise = ((rng_state >> 33) as f64 / u32::MAX as f64) * 2.0 - 1.0;
value += noise * amplitude;
amplitude *= persistence;
frequency *= lacunarity;
}
let elevation = value.clamp(-1.0, 1.0);
let mapped = if elevation < 0.0 {
elevation * self.min_elevation_m.abs()
} else {
elevation * self.max_elevation_m
};
let idx = (y * self.resolution + x) as usize;
self.data[idx] = mapped;
}
}
}
}