use sciforge::hub::prelude::constants::EARTH_RADIUS;
fn legendre_pnm(n: u32, m: u32, x: f64) -> f64 {
if m > n {
return 0.0;
}
let u = (1.0 - x * x).sqrt().max(0.0);
let mut pmm = 1.0_f64;
for i in 1..=m {
let fi = i as f64;
pmm *= u * ((2.0 * fi + 1.0) / (2.0 * fi)).sqrt();
}
if n == m {
return pmm;
}
let pm1m = x * ((2 * m + 3) as f64).sqrt() * pmm;
if n == m + 1 {
return pm1m;
}
let mut p_prev2 = pmm;
let mut p_prev1 = pm1m;
let mut result = 0.0;
for l in (m + 2)..=n {
let fl = l as f64;
let fm = m as f64;
let a = ((4.0 * fl * fl - 1.0) / (fl * fl - fm * fm)).sqrt();
let b =
(((fl - 1.0) * (fl - 1.0) - fm * fm) / (4.0 * (fl - 1.0) * (fl - 1.0) - 1.0)).sqrt();
result = a * (x * p_prev1 - b * p_prev2);
p_prev2 = p_prev1;
p_prev1 = result;
}
result
}
fn geoid_undulation(lat_deg: f64, lon_deg: f64) -> f64 {
let lat_rad = lat_deg.to_radians();
let lon_rad = lon_deg.to_radians();
let cos_colat = lat_rad.sin();
const COEFFS: &[(u32, u32, f64, f64)] = &[
(2, 0, -484.165, 0.0),
(2, 1, -0.000, 0.001),
(2, 2, 2.439, -1.400),
(3, 0, 0.957, 0.0),
(3, 1, 2.030, 0.249),
(3, 2, 0.905, -0.619),
(3, 3, 0.721, 1.414),
(4, 0, 0.540, 0.0),
(4, 1, -0.536, -0.474),
(4, 2, 0.351, 0.663),
(4, 3, 0.991, -0.201),
(4, 4, -0.189, 0.309),
(5, 0, 0.069, 0.0),
(5, 1, -0.062, -0.094),
(5, 2, 0.652, -0.323),
(5, 3, -0.452, -0.215),
(5, 4, -0.295, 0.050),
(5, 5, 0.175, -0.669),
(6, 0, -0.150, 0.0),
(6, 1, -0.076, 0.027),
(6, 2, 0.049, -0.374),
(6, 3, 0.057, 0.009),
(6, 4, -0.087, -0.471),
(6, 5, -0.267, -0.536),
(6, 6, 0.010, -0.237),
(7, 0, 0.091, 0.0),
(7, 1, 0.028, 0.095),
(7, 2, 0.330, 0.093),
(7, 3, 0.251, -0.217),
(7, 4, -0.275, -0.124),
(7, 5, 0.002, 0.016),
(7, 6, -0.359, 0.152),
(7, 7, 0.011, 0.024),
(8, 0, 0.050, 0.0),
(8, 1, 0.023, 0.059),
(8, 2, 0.080, 0.065),
(8, 3, -0.002, -0.086),
(8, 4, -0.250, 0.070),
(9, 0, 0.028, 0.0),
(9, 1, 0.142, 0.022),
(9, 2, 0.029, -0.002),
(9, 3, -0.153, 0.052),
(9, 4, 0.003, -0.059),
(9, 5, -0.013, -0.005),
(10, 0, 0.053, 0.0),
(10, 1, 0.083, -0.013),
(10, 2, -0.078, -0.088),
(10, 3, -0.048, 0.017),
(10, 4, -0.012, -0.050),
(10, 5, -0.027, 0.054),
(11, 0, -0.051, 0.0),
(11, 1, 0.006, -0.033),
(11, 2, -0.009, -0.027),
(11, 3, -0.024, -0.034),
(11, 4, -0.007, 0.016),
(12, 0, -0.036, 0.0),
(12, 1, -0.029, -0.005),
(12, 2, 0.038, 0.025),
(12, 3, -0.012, -0.004),
(12, 4, 0.017, -0.028),
(12, 5, 0.010, -0.001),
(12, 6, -0.043, -0.058),
];
let mut undulation = 0.0;
for &(n, m, cnm, snm) in COEFFS {
let pnm = legendre_pnm(n, m, cos_colat);
let cos_mlon = (m as f64 * lon_rad).cos();
let sin_mlon = (m as f64 * lon_rad).sin();
undulation += pnm * (cnm * cos_mlon + snm * sin_mlon);
}
undulation * 1e-6 * EARTH_RADIUS
}
pub struct Heightmap {
pub resolution: u32,
pub data: Vec<f64>,
pub min_elevation_m: f64,
pub max_elevation_m: f64,
}
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 smooth_step(x: f64, edge0: f64, edge1: f64) -> f64 {
let t = ((x - edge0) / (edge1 - edge0)).clamp(0.0, 1.0);
t * t * (3.0 - 2.0 * t)
}
fn land_fraction(lat: f64, lon: f64) -> f64 {
const CONTINENTS: &[(f64, f64, f64, f64)] = &[
(5.0, 22.0, 20.0, 15.0),
(8.0, 45.0, 8.0, 5.0),
(50.0, 15.0, 12.0, 20.0),
(63.0, 15.0, 7.0, 8.0),
(45.0, 80.0, 20.0, 40.0),
(20.0, 78.0, 10.0, 7.0),
(15.0, 105.0, 10.0, 10.0),
(35.0, 115.0, 12.0, 15.0),
(62.0, 100.0, 10.0, 35.0),
(45.0, -100.0, 15.0, 25.0),
(64.0, -153.0, 5.0, 12.0),
(15.0, -88.0, 5.0, 5.0),
(-5.0, -60.0, 12.0, 15.0),
(-35.0, -65.0, 15.0, 8.0),
(-25.0, 134.0, 12.0, 18.0),
(-82.0, 0.0, 8.0, 180.0),
(72.0, -42.0, 8.0, 10.0),
(23.0, 45.0, 8.0, 8.0),
(-2.0, 115.0, 4.0, 15.0),
(54.0, -2.0, 4.0, 3.0),
(-19.0, 47.0, 6.0, 2.5),
(-42.0, 172.0, 5.0, 2.0),
(36.0, 138.0, 5.0, 3.0),
(12.0, 122.0, 6.0, 3.0),
];
let mut max_frac = 0.0_f64;
for &(clat, clon, slat, slon) in CONTINENTS {
let dlat = (lat - clat) / slat;
let dlon = shortest_lon_diff(lon, clon) / slon;
let d2 = dlat * dlat + dlon * dlon;
let frac = (-d2 * 0.5).exp();
max_frac = max_frac.max(frac);
}
smooth_step(max_frac, 0.15, 0.35)
}
fn mountain_elevation(lat: f64, lon: f64) -> f64 {
const MOUNTAINS: &[(f64, f64, f64, f64, f64)] = &[
(28.0, 85.0, 3.0, 12.0, 5500.0),
(33.0, 88.0, 8.0, 15.0, 4500.0),
(36.0, 76.0, 2.0, 5.0, 5000.0),
(-15.0, -70.0, 30.0, 3.0, 4000.0),
(45.0, -110.0, 15.0, 5.0, 3000.0),
(46.5, 10.0, 2.0, 6.0, 2500.0),
(58.0, 60.0, 12.0, 2.0, 1000.0),
(-2.0, 37.0, 8.0, 4.0, 2000.0),
(10.0, 38.0, 5.0, 5.0, 2000.0),
(37.0, -80.0, 8.0, 3.0, 1000.0),
(63.0, 14.0, 8.0, 3.0, 1200.0),
(32.0, -1.0, 3.0, 8.0, 2000.0),
(-30.0, 149.0, 12.0, 3.0, 1000.0),
(42.5, 44.0, 2.0, 5.0, 3000.0),
(64.0, -153.0, 3.0, 5.0, 3500.0),
(-3.5, 138.0, 3.0, 2.0, 3500.0),
(72.0, -42.0, 5.0, 6.0, 2000.0),
(-82.0, 0.0, 6.0, 90.0, 2500.0),
];
let mut elev = 0.0_f64;
for &(clat, clon, slat, slon, peak) in MOUNTAINS {
let dlat = (lat - clat) / slat;
let dlon = shortest_lon_diff(lon, clon) / slon;
let d2 = dlat * dlat + dlon * dlon;
let contribution = peak * (-d2 * 0.5).exp();
elev = elev.max(contribution);
}
elev
}
fn trench_depth(lat: f64, lon: f64) -> f64 {
const TRENCHES: &[(f64, f64, f64, f64, f64)] = &[
(11.35, 142.2, 2.0, 1.0, -7000.0),
(-22.0, -175.0, 3.0, 1.0, -7000.0),
(10.0, 127.0, 3.0, 1.0, -6500.0),
(19.7, -66.0, 1.0, 4.0, -4500.0),
(-55.0, -26.0, 3.0, 1.0, -4400.0),
(-10.0, 110.0, 1.0, 10.0, -3900.0),
];
let mut depth = 0.0_f64;
for &(clat, clon, slat, slon, d) in TRENCHES {
let dlat = (lat - clat) / slat;
let dlon = shortest_lon_diff(lon, clon) / slon;
let d2 = dlat * dlat + dlon * dlon;
let contribution = d * (-d2 * 0.5).exp();
depth = depth.min(contribution);
}
depth
}
fn ridge_elevation(lat: f64, lon: f64) -> f64 {
const RIDGES: &[(f64, f64, f64, f64, f64)] = &[
(30.0, -40.0, 30.0, 3.0, 1500.0),
(-20.0, -14.0, 30.0, 3.0, 1500.0),
(-50.0, 140.0, 10.0, 40.0, 1200.0),
(10.0, 60.0, 15.0, 3.0, 1000.0),
(-30.0, -110.0, 40.0, 3.0, 1500.0),
];
let mut elev = 0.0_f64;
for &(clat, clon, slat, slon, h) in RIDGES {
let dlat = (lat - clat) / slat;
let dlon = shortest_lon_diff(lon, clon) / slon;
let d2 = dlat * dlat + dlon * dlon;
let contribution = h * (-d2 * 0.5).exp();
elev = elev.max(contribution);
}
elev
}
pub fn earth_elevation(lat: f64, lon: f64) -> f64 {
let geoid = geoid_undulation(lat, lon);
let land = land_fraction(lat, lon);
let base_ocean = -3800.0;
let base_continent = 200.0;
let base = base_ocean * (1.0 - land) + base_continent * land;
let mtn = mountain_elevation(lat, lon);
let trench = trench_depth(lat, lon);
let ridge = ridge_elevation(lat, lon) * (1.0 - land);
let elev = base + mtn + trench + ridge + geoid;
elev.clamp(-10_994.0, 8_848.86)
}
impl Heightmap {
pub fn new(resolution: u32) -> Self {
let size = (resolution * resolution) as usize;
Self {
resolution,
data: vec![0.0; size],
min_elevation_m: -10_994.0,
max_elevation_m: 8_848.86,
}
}
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 {
EARTH_RADIUS + self.sample(lat_deg, lon_deg)
}
pub fn generate_earth_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] = earth_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;
}
}
}
}