use sciforge::hub::prelude::constants::EARTH_RADIUS;
fn legendrepnm(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.0f64;
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 pprev2 = pmm;
let mut pprev1 = 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 * pprev1 - b * pprev2);
pprev2 = pprev1;
pprev1 = result;
}
result
}
fn geoidundulation(latdeg: f64, londeg: f64) -> f64 {
let latrad = latdeg.to_radians();
let lonrad = londeg.to_radians();
let coscolat = latrad.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 = legendrepnm(n, m, coscolat);
let cosmlon = (m as f64 * lonrad).cos();
let sinmlon = (m as f64 * lonrad).sin();
undulation += pnm * (cnm * cosmlon + snm * sinmlon);
}
undulation * 1e-6 * EARTH_RADIUS
}
pub struct Heightmap {
pub resolution: u32,
pub data: Vec<f64>,
pub minelevationm: f64,
pub maxelevationm: f64,
}
fn shortestlondiff(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 smoothstep(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 landfraction(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 maxfrac = 0.0f64;
for &(clat, clon, slat, slon) in CONTINENTS {
let dlat = (lat - clat) / slat;
let dlon = shortestlondiff(lon, clon) / slon;
let d2 = dlat * dlat + dlon * dlon;
let frac = (-d2 * 0.5).exp();
maxfrac = maxfrac.max(frac);
}
smoothstep(maxfrac, 0.15, 0.35)
}
fn mountainelevation(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.0f64;
for &(clat, clon, slat, slon, peak) in MOUNTAINS {
let dlat = (lat - clat) / slat;
let dlon = shortestlondiff(lon, clon) / slon;
let d2 = dlat * dlat + dlon * dlon;
let contribution = peak * (-d2 * 0.5).exp();
elev = elev.max(contribution);
}
elev
}
fn trenchdepth(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.0f64;
for &(clat, clon, slat, slon, d) in TRENCHES {
let dlat = (lat - clat) / slat;
let dlon = shortestlondiff(lon, clon) / slon;
let d2 = dlat * dlat + dlon * dlon;
let contribution = d * (-d2 * 0.5).exp();
depth = depth.min(contribution);
}
depth
}
fn ridgeelevation(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.0f64;
for &(clat, clon, slat, slon, h) in RIDGES {
let dlat = (lat - clat) / slat;
let dlon = shortestlondiff(lon, clon) / slon;
let d2 = dlat * dlat + dlon * dlon;
let contribution = h * (-d2 * 0.5).exp();
elev = elev.max(contribution);
}
elev
}
pub fn earthelevation(lat: f64, lon: f64) -> f64 {
let geoid = geoidundulation(lat, lon);
let land = landfraction(lat, lon);
let baseocean = -3800.0;
let basecontinent = 200.0;
let base = baseocean * (1.0 - land) + basecontinent * land;
let mtn = mountainelevation(lat, lon);
let trench = trenchdepth(lat, lon);
let ridge = ridgeelevation(lat, lon) * (1.0 - land);
let elev = base + mtn + trench + ridge + geoid;
elev.clamp(-10994.0, 8848.86)
}
impl Heightmap {
pub fn new(resolution: u32) -> Self {
let size = (resolution * resolution) as usize;
Self {
resolution,
data: vec![0.0; size],
minelevationm: -10994.0,
maxelevationm: 8848.86,
}
}
pub fn sample(&self, latdeg: f64, londeg: f64) -> f64 {
let u = ((londeg + 180.0) / 360.0).clamp(0.0, 1.0);
let v = ((latdeg + 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 radiusat(&self, latdeg: f64, londeg: f64) -> f64 {
EARTH_RADIUS + self.sample(latdeg, londeg)
}
pub fn generateearthtopography(&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] = earthelevation(lat, lon);
}
}
}
pub fn generateprocedural(
&mut self,
seed: u64,
octaves: u32,
persistence: f64,
lacunarity: f64,
) {
let mut rngstate = 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 octave in 0..octaves {
rngstate = rngstate
.wrapping_mul(6364136223846793005)
.wrapping_add(frequency as u64 + octave as u64);
let noise = ((rngstate >> 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.minelevationm.abs()
} else {
elevation * self.maxelevationm
};
let idx = (y * self.resolution + x) as usize;
self.data[idx] = mapped;
}
}
}
}