From e91578ffdb838f93480aeed48c71bfc339b65296 Mon Sep 17 00:00:00 2001 From: Joshua Yanovski Date: Sat, 11 Jan 2020 21:38:30 +0100 Subject: [PATCH] Cargo fmt most things (except erosion.rs). --- client/src/lib.rs | 7 +- common/src/msg/server.rs | 2 +- server/src/lib.rs | 24 +- server/src/settings.rs | 3 +- voxygen/src/hud/map.rs | 2 +- voxygen/src/hud/minimap.rs | 2 +- voxygen/src/hud/mod.rs | 5 +- voxygen/src/scene/mod.rs | 2 +- voxygen/src/scene/terrain.rs | 4 +- world/examples/view.rs | 11 +- world/examples/water.rs | 104 +++--- world/src/block/mod.rs | 3 +- world/src/column/mod.rs | 37 +- world/src/sim/diffusion.rs | 550 ++++++++++++++-------------- world/src/sim/mod.rs | 683 ++++++++++++++++++----------------- world/src/sim/util.rs | 24 +- 16 files changed, 768 insertions(+), 695 deletions(-) diff --git a/client/src/lib.rs b/client/src/lib.rs index 4264dceb80..7a0b633dc9 100644 --- a/client/src/lib.rs +++ b/client/src/lib.rs @@ -113,11 +113,8 @@ impl Client { log::debug!("Preparing image..."); let world_map = Arc::new(image::DynamicImage::ImageRgba8({ // Should not fail if the dimensions are correct. - let world_map = image::ImageBuffer::from_raw( - map_size.x, - map_size.y, - world_map_raw, - ); + let world_map = + image::ImageBuffer::from_raw(map_size.x, map_size.y, world_map_raw); world_map.ok_or(Error::Other("Server sent a bad world map image".into()))? })); log::debug!("Done preparing image..."); diff --git a/common/src/msg/server.rs b/common/src/msg/server.rs index 8636ffa25e..6e3f9ad661 100644 --- a/common/src/msg/server.rs +++ b/common/src/msg/server.rs @@ -37,7 +37,7 @@ pub enum ServerMsg { entity_package: sync::EntityPackage, server_info: ServerInfo, time_of_day: state::TimeOfDay, - world_map: (Vec2, Vec) + world_map: (Vec2, Vec), }, PlayerListUpdate(PlayerListUpdate), StateAnswer(Result), diff --git a/server/src/lib.rs b/server/src/lib.rs index 2c42102011..2b04f88b62 100644 --- a/server/src/lib.rs +++ b/server/src/lib.rs @@ -46,7 +46,10 @@ use std::{ }; use uvth::{ThreadPool, ThreadPoolBuilder}; use vek::*; -use world::{sim::{FileOpts, WORLD_SIZE, WorldOpts}, World}; +use world::{ + sim::{FileOpts, WorldOpts, WORLD_SIZE}, + World, +}; const CLIENT_TIMEOUT: f64 = 20.0; // Seconds pub enum Event { @@ -104,15 +107,18 @@ impl Server { state.ecs_mut().register::(); state.ecs_mut().register::(); - let world = World::generate(settings.world_seed, WorldOpts { - seed_elements: true, - world_file: if let Some(ref file) = settings.map_file { - FileOpts::Load(file.clone()) - } else { - FileOpts::Generate + let world = World::generate( + settings.world_seed, + WorldOpts { + seed_elements: true, + world_file: if let Some(ref file) = settings.map_file { + FileOpts::Load(file.clone()) + } else { + FileOpts::Generate + }, + ..WorldOpts::default() }, - ..WorldOpts::default() - }); + ); let spawn_point = { // NOTE: all of these `.map(|e| e as [type])` calls should compile into no-ops, diff --git a/server/src/settings.rs b/server/src/settings.rs index 134fc10e50..92d6f7409b 100644 --- a/server/src/settings.rs +++ b/server/src/settings.rs @@ -104,8 +104,7 @@ impl ServerSettings { max_players: 100, start_time: 9.0 * 3600.0, admins: vec!["singleplayer".to_string()], // TODO: Let the player choose if they want to use admin commands or not - .. - Self::load() // Fill in remaining fields from settings.ron. + ..Self::load() // Fill in remaining fields from settings.ron. } } diff --git a/voxygen/src/hud/map.rs b/voxygen/src/hud/map.rs index fdc11fdf42..e1f4c613e7 100644 --- a/voxygen/src/hud/map.rs +++ b/voxygen/src/hud/map.rs @@ -162,7 +162,7 @@ impl<'a> Widget for Map<'a> { let (world_map, worldsize) = self.world_map; let worldsize = worldsize.map2(TerrainChunkSize::RECT_SIZE, |e, f| e as f64 * f as f64); - Image::new(world_map/*self.imgs.map_placeholder*/) + Image::new(world_map /*self.imgs.map_placeholder*/) .middle_of(state.ids.map_bg) .color(Some(Color::Rgba(1.0, 1.0, 1.0, fade - 0.1))) .w_h(700.0, 700.0) diff --git a/voxygen/src/hud/minimap.rs b/voxygen/src/hud/minimap.rs index 2bbf2702b2..a5597661a9 100644 --- a/voxygen/src/hud/minimap.rs +++ b/voxygen/src/hud/minimap.rs @@ -137,7 +137,7 @@ impl<'a> Widget for MiniMap<'a> { // Map Image let (world_map, worldsize) = self.world_map; let worldsize = worldsize.map2(TerrainChunkSize::RECT_SIZE, |e, f| e as f64 * f as f64); - Image::new(world_map/*self.imgs.map_placeholder*/) + Image::new(world_map /*self.imgs.map_placeholder*/) .middle_of(state.ids.mmap_frame_bg) .w_h(92.0 * 4.0 * zoom, 82.0 * 4.0 * zoom) .parent(state.ids.mmap_frame_bg) diff --git a/voxygen/src/hud/mod.rs b/voxygen/src/hud/mod.rs index c6b04c2bb0..2bf5e07b25 100644 --- a/voxygen/src/hud/mod.rs +++ b/voxygen/src/hud/mod.rs @@ -454,7 +454,10 @@ impl Hud { // Generate ids. let ids = Ids::new(ui.id_generator()); // Load world map - let world_map = (ui.add_graphic(Graphic::Image(client.world_map.0.clone())), client.world_map.1); + let world_map = ( + ui.add_graphic(Graphic::Image(client.world_map.0.clone())), + client.world_map.1, + ); // Load images. let imgs = Imgs::load(&mut ui).expect("Failed to load images!"); // Load rotation images. diff --git a/voxygen/src/scene/mod.rs b/voxygen/src/scene/mod.rs index f0691d18af..4faa0f7bf7 100644 --- a/voxygen/src/scene/mod.rs +++ b/voxygen/src/scene/mod.rs @@ -20,7 +20,7 @@ use client::Client; use common::{ comp, terrain::{BlockKind, TerrainChunk, TerrainChunkSize}, - vol::{RectVolSize, ReadVol}, + vol::{ReadVol, RectVolSize}, }; use specs::{Join, WorldExt}; use vek::*; diff --git a/voxygen/src/scene/terrain.rs b/voxygen/src/scene/terrain.rs index 5726f2bc28..536714f650 100644 --- a/voxygen/src/scene/terrain.rs +++ b/voxygen/src/scene/terrain.rs @@ -175,7 +175,9 @@ fn mesh_worker + RectRasterableVol + ReadVol + Debug>( let kind = volume.get(wpos).unwrap_or(&Block::empty()).kind(); if let Some(cfg) = sprite_config_for(kind) { - let seed = wpos.x as u64 * 3 + wpos.y as u64 * 7 + wpos.x as u64 * wpos.y as u64; // Awful PRNG + let seed = wpos.x as u64 * 3 + + wpos.y as u64 * 7 + + wpos.x as u64 * wpos.y as u64; // Awful PRNG let instance = SpriteInstance::new( Mat4::identity() diff --git a/world/examples/view.rs b/world/examples/view.rs index eab8ef1398..6ff94ee124 100644 --- a/world/examples/view.rs +++ b/world/examples/view.rs @@ -6,10 +6,13 @@ const W: usize = 640; const H: usize = 480; fn main() { - let world = World::generate(0, WorldOpts { - seed_elements: true, - ..WorldOpts::default() - }); + let world = World::generate( + 0, + WorldOpts { + seed_elements: true, + ..WorldOpts::default() + }, + ); let sampler = world.sample_columns(); diff --git a/world/examples/water.rs b/world/examples/water.rs index 7b3c977ebe..68a7d48aa3 100644 --- a/world/examples/water.rs +++ b/world/examples/water.rs @@ -36,12 +36,15 @@ fn main() { let mut _map_file = PathBuf::from("./maps"); _map_file.push(map_file); - let world = World::generate(1337, WorldOpts { - seed_elements: false, - // world_file: sim::FileOpts::Load(_map_file), - world_file: sim::FileOpts::Save, - ..WorldOpts::default() - }); + let world = World::generate( + 1337, + WorldOpts { + seed_elements: false, + // world_file: sim::FileOpts::Load(_map_file), + world_file: sim::FileOpts::Save, + ..WorldOpts::default() + }, + ); let sampler = world.sim(); @@ -87,25 +90,35 @@ fn main() { for i in 0..W { for j in 0..H { - let pos = (focus_rect + Vec2::new(i as f64, j as f64) * scale).map(|e: f64| e as i32); + let pos = + (focus_rect + Vec2::new(i as f64, j as f64) * scale).map(|e: f64| e as i32); /* let top_left = pos; let top_right = focus + Vec2::new(i as i32 + light_res, j as i32) * scale; let bottom_left = focus + Vec2::new(i as i32, j as i32 + light_res) * scale; */ - let (alt, basement, water_alt, humidity, temperature, downhill, river_kind) = sampler - .get(pos) - .map(|sample| { - ( - sample.alt, - sample.basement, - sample.water_alt, - sample.humidity, - sample.temp, - sample.downhill, - sample.river.river_kind, - ) - }) - .unwrap_or((CONFIG.sea_level, CONFIG.sea_level, CONFIG.sea_level, 0.0, 0.0, None, None)); + let (alt, basement, water_alt, humidity, temperature, downhill, river_kind) = + sampler + .get(pos) + .map(|sample| { + ( + sample.alt, + sample.basement, + sample.water_alt, + sample.humidity, + sample.temp, + sample.downhill, + sample.river.river_kind, + ) + }) + .unwrap_or(( + CONFIG.sea_level, + CONFIG.sea_level, + CONFIG.sea_level, + 0.0, + 0.0, + None, + None, + )); let humidity = humidity.min(1.0).max(0.0); let temperature = temperature.min(1.0).max(-1.0) * 0.5 + 0.5; let pos = pos * TerrainChunkSize::RECT_SIZE.map(|e| e as i32); @@ -177,15 +190,9 @@ fn main() { let true_water_alt = (alt.max(water_alt) as f64 - focus.z) / gain as f64; let true_alt = (alt as f64 - focus.z) / gain as f64; - let water_depth = (true_water_alt - true_alt) - .min(1.0) - .max(0.0); - let water_alt = true_water_alt - .min(1.0) - .max(0.0); - let alt = true_alt - .min(1.0) - .max(0.0); + let water_depth = (true_water_alt - true_alt).min(1.0).max(0.0); + let water_alt = true_water_alt.min(1.0).max(0.0); + let alt = true_alt.min(1.0).max(0.0); let quad = |x: f32| ((x as f64 * QUADRANTS as f64).floor() as usize).min(QUADRANTS - 1); if river_kind.is_none() || humidity != 0.0 { @@ -205,17 +212,29 @@ fn main() { } buf[j * W + i] = match (river_kind, (is_water, true_alt >= true_sea_level)) { - (_, (false, _)) | ( None, (_, true)) => { + (_, (false, _)) | (None, (_, true)) => { let (r, g, b) = ( - (if is_shaded { alt } else { alt } * if is_temperature { temperature as f64 } else if is_shaded { alt } else { 0.0 }).sqrt(), + (if is_shaded { alt } else { alt } + * if is_temperature { + temperature as f64 + } else if is_shaded { + alt + } else { + 0.0 + }) + .sqrt(), if is_shaded { 0.2 + (alt * 0.8) } else { alt }, - (if is_shaded { alt } else { alt } * if is_humidity { humidity as f64 } else if is_shaded { alt } else { 0.0 }).sqrt(), + (if is_shaded { alt } else { alt } + * if is_humidity { + humidity as f64 + } else if is_shaded { + alt + } else { + 0.0 + }) + .sqrt(), ); - let light = if is_shaded { - light - } else { - 1.0 - }; + let light = if is_shaded { light } else { 1.0 }; u32::from_le_bytes([ (b * light * 255.0) as u8, (g * light * 255.0) as u8, @@ -228,7 +247,7 @@ fn main() { (/*alt*//*alt * *//*(1.0 - humidity)*/(alt * temperature).sqrt() * 255.0) as u8, 255, ]) */ - }, + } (Some(RiverKind::Ocean), _) => u32::from_le_bytes([ ((64.0 - water_depth * 64.0) * 1.0) as u8, ((32.0 - water_depth * 32.0) * 1.0) as u8, @@ -242,8 +261,8 @@ fn main() { 255, ]), (None, _) | (Some(RiverKind::Lake { .. }), _) => u32::from_le_bytes([ - (((64.0 + water_alt * 191.0) + (- water_depth * 64.0)) * 1.0) as u8, - (((32.0 + water_alt * 95.0) + (- water_depth * 32.0)) * 1.0) as u8, + (((64.0 + water_alt * 191.0) + (-water_depth * 64.0)) * 1.0) as u8, + (((32.0 + water_alt * 95.0) + (-water_depth * 32.0)) * 1.0) as u8, 0, 255, ]), @@ -280,7 +299,8 @@ fn main() { } if win.get_mouse_down(minifb::MouseButton::Left) { if let Some((mx, my)) = win.get_mouse_pos(minifb::MouseMode::Clamp) { - let pos = (focus_rect + (Vec2::new(mx as f64, my as f64) * scale)).map(|e| e as i32); + let pos = + (focus_rect + (Vec2::new(mx as f64, my as f64) * scale)).map(|e| e as i32); println!( "Chunk position: {:?}", pos.map2(TerrainChunkSize::RECT_SIZE, |e, f| e * f as i32) diff --git a/world/src/block/mod.rs b/world/src/block/mod.rs index 0aefcd7b72..83d6553ff7 100644 --- a/world/src/block/mod.rs +++ b/world/src/block/mod.rs @@ -265,7 +265,8 @@ impl<'a> BlockGen<'a> { sub_surface_color, stone_col.map(|e| e as f32 / 255.0), (height - grass_depth - wposf.z as f32) * 0.15, - ).map(|e| (e * 255.0) as u8); + ) + .map(|e| (e * 255.0) as u8); // Underground if (wposf.z as f32) > alt - 32.0 * chaos { diff --git a/world/src/column/mod.rs b/world/src/column/mod.rs index bb006f148e..0a2f859572 100644 --- a/world/src/column/mod.rs +++ b/world/src/column/mod.rs @@ -437,13 +437,14 @@ impl<'a> Sampler<'a> for ColumnGen<'a> { downhill_alt.sub(CONFIG.sea_level) >= CONFIG.mountain_scale * 0.25)*/ is_rocky { sim.get_interpolated_monotone(wpos, |chunk| chunk.alt)? - // sim.get_interpolated_bilinear(wpos, |chunk| chunk.alt)? - // sim.get_interpolated(wpos, |chunk| chunk.alt)? + // sim.get_interpolated_bilinear(wpos, |chunk| chunk.alt)? + // sim.get_interpolated(wpos, |chunk| chunk.alt)? } else { sim.get_interpolated_monotone(wpos, |chunk| chunk.alt)? // sim.get_interpolated(wpos, |chunk| chunk.alt)? }; - let basement = alt + sim./*get_interpolated*/get_interpolated_monotone(wpos, |chunk| chunk.basement.sub(chunk.alt))?; + let basement = alt + + sim./*get_interpolated*/get_interpolated_monotone(wpos, |chunk| chunk.basement.sub(chunk.alt))?; // Find the average distance to each neighboring body of water. let mut river_count = 0.0f64; @@ -830,7 +831,6 @@ impl<'a> Sampler<'a> for ColumnGen<'a> { .mul(1.0 - humidity) */ /* .mul(32.0) */; - let riverless_alt_delta = Lerp::lerp(0.0, riverless_alt_delta, warp_factor); let alt = alt_ + riverless_alt_delta; let basement = basement.min(alt); @@ -895,7 +895,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> { ); let tundra = Lerp::lerp(snow, Rgb::new(0.01, 0.3, 0.0), 0.4 + marble * 0.6); let dead_tundra = Lerp::lerp(warm_stone, Rgb::new(0.3, 0.12, 0.2), marble); - let cliff = Rgb::lerp(cold_stone, /*warm_stone*/hot_stone, marble); + let cliff = Rgb::lerp(cold_stone, /*warm_stone*/ hot_stone, marble); let grass = Rgb::lerp( cold_grass, @@ -953,11 +953,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> { .mul(1.0), ); - let sub_surface_color = Lerp::lerp( - cliff, - ground, - alt.sub(basement).mul(0.25) - ); + let sub_surface_color = Lerp::lerp(cliff, ground, alt.sub(basement).mul(0.25)); /* let ground = Rgb::lerp( dead_tundra, @@ -1087,15 +1083,18 @@ impl<'a> Sampler<'a> for ColumnGen<'a> { ); */ // Snow covering - let snow_cover = - temp.sub(CONFIG.snow_temp) - .max(-humidity.sub(CONFIG.desert_hum)) - .mul(16.0) - .add((marble_small - 0.5) * 0.5); + let snow_cover = temp + .sub(CONFIG.snow_temp) + .max(-humidity.sub(CONFIG.desert_hum)) + .mul(16.0) + .add((marble_small - 0.5) * 0.5); let (alt, ground, sub_surface_color) = if snow_cover /*< 0.1*/<= 0.5 && alt > water_level { // Allow snow cover. - (alt + 1.0 - snow_cover.max(0.0), Rgb::lerp(snow, ground, snow_cover), - Lerp::lerp(sub_surface_color, ground, alt.sub(basement).mul(0.15))) + ( + alt + 1.0 - snow_cover.max(0.0), + Rgb::lerp(snow, ground, snow_cover), + Lerp::lerp(sub_surface_color, ground, alt.sub(basement).mul(0.15)), + ) } else { (alt, ground, sub_surface_color) }; @@ -1165,7 +1164,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> { // Beach ((ocean_level - 1.0) / 2.0).max(0.0), ), - sub_surface_color,// /*warm_grass*/Lerp::lerp(cliff, dirt, alt.sub(basement).mul(0.25)), + sub_surface_color, // /*warm_grass*/Lerp::lerp(cliff, dirt, alt.sub(basement).mul(0.25)), // No growing directly on bedrock. tree_density: Lerp::lerp(0.0, tree_density, alt.sub(2.0).sub(basement).mul(0.5)), forest_kind: sim_chunk.forest_kind, @@ -1183,7 +1182,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> { humidity, spawn_rate, location: sim_chunk.location.as_ref(), - stone_col, + stone_col, chunk: sim_chunk, spawn_rules: sim_chunk diff --git a/world/src/sim/diffusion.rs b/world/src/sim/diffusion.rs index eebeb5e229..5e7cc82575 100644 --- a/world/src/sim/diffusion.rs +++ b/world/src/sim/diffusion.rs @@ -1,5 +1,5 @@ -use rayon::prelude::*; use super::Alt; +use rayon::prelude::*; /// From https://github.com/fastscape-lem/fastscapelib-fortran/blob/master/src/Diffusion.f90 /// @@ -36,68 +36,76 @@ use super::Alt; implicit none */ -pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (), - h: &mut [Alt], b: &mut [Alt], - kd: impl Fn(usize) -> f64, kdsed: f64, - ) { +pub fn diffusion( + nx: usize, + ny: usize, + xl: f64, + yl: f64, + dt: f64, + _ibc: (), + h: &mut [Alt], + b: &mut [Alt], + kd: impl Fn(usize) -> f64, + kdsed: f64, +) { let aij = |i: usize, j: usize| j * nx + i; -/* - double precision, dimension(:), allocatable :: f,diag,sup,inf,res - double precision, dimension(:,:), allocatable :: zint,kdint,zintp - integer i,j,ij - double precision factxp,factxm,factyp,factym,dx,dy -*/ - let mut f : Vec; - let mut diag : Vec; - let mut sup : Vec; - let mut inf : Vec; - let mut res : Vec; - let mut zint : Vec; - let mut kdint : Vec; - let mut zintp : Vec; - let mut i : usize; - let mut j : usize; - let mut ij : usize; - let mut factxp : f64; - let mut factxm : f64; - let mut factyp : f64; - let mut factym : f64; - let mut dx : f64; - let mut dy : f64; -/* - character cbc*4 + /* + double precision, dimension(:), allocatable :: f,diag,sup,inf,res + double precision, dimension(:,:), allocatable :: zint,kdint,zintp + integer i,j,ij + double precision factxp,factxm,factyp,factym,dx,dy + */ + let mut f: Vec; + let mut diag: Vec; + let mut sup: Vec; + let mut inf: Vec; + let mut res: Vec; + let mut zint: Vec; + let mut kdint: Vec; + let mut zintp: Vec; + let mut i: usize; + let mut j: usize; + let mut ij: usize; + let mut factxp: f64; + let mut factxm: f64; + let mut factyp: f64; + let mut factym: f64; + let mut dx: f64; + let mut dy: f64; + /* + character cbc*4 - !print*,'Diffusion' + !print*,'Diffusion' - write (cbc,'(i4)') ibc + write (cbc,'(i4)') ibc - dx=xl/(nx-1) - dy=yl/(ny-1) -*/ + dx=xl/(nx-1) + dy=yl/(ny-1) + */ // 2048*32/2048/2048 // 1 / 64 m dx = xl / /*(nx - 1)*/nx as f64; dy = yl / /*(ny - 1)*/ny as f64; -/* - ! creates 2D internal arrays to store topo and kd + /* + ! creates 2D internal arrays to store topo and kd - allocate (zint(nx,ny),kdint(nx,ny),zintp(nx,ny)) -*/ + allocate (zint(nx,ny),kdint(nx,ny),zintp(nx,ny)) + */ zint = vec![Default::default(); nx * ny]; kdint = vec![Default::default(); nx * ny]; zintp = vec![Default::default(); nx * ny]; -/* - do j=1,ny - do i=1,nx - ij=(j-1)*nx+i - zint(i,j)=h(ij) - kdint(i,j)=kd(ij) - if (kdsed.gt.0.d0 .and. (h(ij)-b(ij)).gt.1.d-6) kdint(i,j)=kdsed - enddo - enddo + /* + do j=1,ny + do i=1,nx + ij=(j-1)*nx+i + zint(i,j)=h(ij) + kdint(i,j)=kd(ij) + if (kdsed.gt.0.d0 .and. (h(ij)-b(ij)).gt.1.d-6) kdint(i,j)=kdsed + enddo + enddo - zintp = zint -*/ + zintp = zint + */ for j in 0..ny { for i in 0..nx { // ij = vec2_as_uniform_idx(i, j); @@ -111,63 +119,62 @@ pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (), } zintp = zint.clone(); -/* - ! first pass along the x-axis + /* + ! first pass along the x-axis - allocate (f(nx),diag(nx),sup(nx),inf(nx),res(nx)) - f=0.d0 - diag=0.d0 - sup=0.d0 - inf=0.d0 - res=0.d0 - do j=2,ny-1 -*/ + allocate (f(nx),diag(nx),sup(nx),inf(nx),res(nx)) + f=0.d0 + diag=0.d0 + sup=0.d0 + inf=0.d0 + res=0.d0 + do j=2,ny-1 + */ f = vec![0.0; nx]; diag = vec![0.0; nx]; sup = vec![0.0; nx]; inf = vec![0.0; nx]; res = vec![0.0; nx]; - for j in 1..ny-1 { -/* - do i=2,nx-1 - factxp=(kdint(i+1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2 - factxm=(kdint(i-1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2 - factyp=(kdint(i,j+1)+kdint(i,j))/2.d0*(dt/2.)/dy**2 - factym=(kdint(i,j-1)+kdint(i,j))/2.d0*(dt/2.)/dy**2 - diag(i)=1.d0+factxp+factxm - sup(i)=-factxp - inf(i)=-factxm - f(i)=zintp(i,j)+factyp*zintp(i,j+1)-(factyp+factym)*zintp(i,j)+factym*zintp(i,j-1) - enddo -*/ - for i in 1..nx-1 { - factxp = (kdint[aij(i+1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx); - factxm = (kdint[aij(i-1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx); - factyp = (kdint[aij(i, j+1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy); - factym = (kdint[aij(i, j-1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy); + for j in 1..ny - 1 { + /* + do i=2,nx-1 + factxp=(kdint(i+1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2 + factxm=(kdint(i-1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2 + factyp=(kdint(i,j+1)+kdint(i,j))/2.d0*(dt/2.)/dy**2 + factym=(kdint(i,j-1)+kdint(i,j))/2.d0*(dt/2.)/dy**2 + diag(i)=1.d0+factxp+factxm + sup(i)=-factxp + inf(i)=-factxm + f(i)=zintp(i,j)+factyp*zintp(i,j+1)-(factyp+factym)*zintp(i,j)+factym*zintp(i,j-1) + enddo + */ + for i in 1..nx - 1 { + factxp = (kdint[aij(i + 1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx); + factxm = (kdint[aij(i - 1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx); + factyp = (kdint[aij(i, j + 1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy); + factym = (kdint[aij(i, j - 1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy); diag[i] = 1.0 + factxp + factxm; sup[i] = -factxp; inf[i] = -factxm; - f[i] = - zintp[aij(i, j)] + factyp * zintp[aij(i, j+1)] - - (factyp + factym) * - zintp[aij(i, j)] + factym * zintp[aij(i, j-1)]; + f[i] = zintp[aij(i, j)] + factyp * zintp[aij(i, j + 1)] + - (factyp + factym) * zintp[aij(i, j)] + + factym * zintp[aij(i, j - 1)]; } -/* -! left bc - if (cbc(4:4).eq.'1') then - diag(1)=1. - sup(1)=0. - f(1)=zintp(1,j) - else - factxp=(kdint(2,j)+kdint(1,j))/2.d0*(dt/2.)/dx**2 - factyp=(kdint(1,j+1)+kdint(1,j))/2.d0*(dt/2.)/dy**2 - factym=(kdint(1,j-1)+kdint(1,j))/2.d0*(dt/2.)/dy**2 - diag(1)=1.d0+factxp - sup(1)=-factxp - f(1)=zintp(1,j)+factyp*zintp(1,j+1)-(factyp+factym)*zintp(1,j)+factym*zintp(1,j-1) - endif -*/ + /* + ! left bc + if (cbc(4:4).eq.'1') then + diag(1)=1. + sup(1)=0. + f(1)=zintp(1,j) + else + factxp=(kdint(2,j)+kdint(1,j))/2.d0*(dt/2.)/dx**2 + factyp=(kdint(1,j+1)+kdint(1,j))/2.d0*(dt/2.)/dy**2 + factym=(kdint(1,j-1)+kdint(1,j))/2.d0*(dt/2.)/dy**2 + diag(1)=1.d0+factxp + sup(1)=-factxp + f(1)=zintp(1,j)+factyp*zintp(1,j+1)-(factyp+factym)*zintp(1,j)+factym*zintp(1,j-1) + endif + */ if true { diag[0] = 1.0; sup[0] = 0.0; @@ -175,119 +182,118 @@ pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (), } else { // reflective boundary factxp = (kdint[aij(1, j)] + kdint[aij(0, j)]) / 2.0 * (dt / 2.0) / (dx * dx); - factyp = (kdint[aij(0, j+1)] + kdint[aij(0, j)]) / 2.0 * (dt / 2.0) / (dy * dy); - factym = (kdint[aij(0, j-1)] + kdint[aij(0, j)]) / 2.0 * (dt / 2.0) / (dy * dy); + factyp = (kdint[aij(0, j + 1)] + kdint[aij(0, j)]) / 2.0 * (dt / 2.0) / (dy * dy); + factym = (kdint[aij(0, j - 1)] + kdint[aij(0, j)]) / 2.0 * (dt / 2.0) / (dy * dy); diag[0] = 1.0 + factxp; sup[0] = -factxp; - f[0] = - zintp[aij(0, j)] + factyp * zintp[aij(0, j+1)] - - (factyp + factym) * - zintp[aij(0, j)] + factym * zintp[aij(0, j-1)]; + f[0] = zintp[aij(0, j)] + factyp * zintp[aij(0, j + 1)] + - (factyp + factym) * zintp[aij(0, j)] + + factym * zintp[aij(0, j - 1)]; } -/* -! right bc - if (cbc(2:2).eq.'1') then - diag(nx)=1. - inf(nx)=0. - f(nx)=zintp(nx,j) - else - factxm=(kdint(nx-1,j)+kdint(nx,j))/2.d0*(dt/2.)/dx**2 - factyp=(kdint(nx,j+1)+kdint(nx,j))/2.d0*(dt/2.)/dy**2 - factym=(kdint(nx,j-1)+kdint(nx,j))/2.d0*(dt/2.)/dy**2 - diag(nx)=1.d0+factxm - inf(nx)=-factxm - f(nx)=zintp(nx,j)+factyp*zintp(nx,j+1)-(factyp+factym)*zintp(nx,j)+factym*zintp(nx,j-1) - endif -*/ + /* + ! right bc + if (cbc(2:2).eq.'1') then + diag(nx)=1. + inf(nx)=0. + f(nx)=zintp(nx,j) + else + factxm=(kdint(nx-1,j)+kdint(nx,j))/2.d0*(dt/2.)/dx**2 + factyp=(kdint(nx,j+1)+kdint(nx,j))/2.d0*(dt/2.)/dy**2 + factym=(kdint(nx,j-1)+kdint(nx,j))/2.d0*(dt/2.)/dy**2 + diag(nx)=1.d0+factxm + inf(nx)=-factxm + f(nx)=zintp(nx,j)+factyp*zintp(nx,j+1)-(factyp+factym)*zintp(nx,j)+factym*zintp(nx,j-1) + endif + */ if true { diag[nx - 1] = 1.0; inf[nx - 1] = 0.0; f[nx - 1] = zintp[aij(nx - 1, j)]; } else { // reflective boundary - factxm = (kdint[aij(nx-2, j)] + kdint[aij(nx-1, j)]) / 2.0 * (dt / 2.0) / (dx * dx); - factyp = (kdint[aij(nx-1, j+1)] + kdint[aij(nx-1, j)]) / 2.0 * (dt / 2.0) / (dy * dy); - factym = (kdint[aij(nx-1, j-1)] + kdint[aij(nx-1, j)]) / 2.0 * (dt / 2.0) / (dy * dy); + factxm = (kdint[aij(nx - 2, j)] + kdint[aij(nx - 1, j)]) / 2.0 * (dt / 2.0) / (dx * dx); + factyp = + (kdint[aij(nx - 1, j + 1)] + kdint[aij(nx - 1, j)]) / 2.0 * (dt / 2.0) / (dy * dy); + factym = + (kdint[aij(nx - 1, j - 1)] + kdint[aij(nx - 1, j)]) / 2.0 * (dt / 2.0) / (dy * dy); diag[nx - 1] = 1.0 + factxm; inf[nx - 1] = -factxm; - f[nx - 1] = - zintp[aij(nx-1, j)] + factyp * zintp[aij(nx-1, j+1)] - - (factyp + factym) * - zintp[aij(nx-1, j)] + factym * zintp[aij(nx-1, j-1)]; + f[nx - 1] = zintp[aij(nx - 1, j)] + factyp * zintp[aij(nx - 1, j + 1)] + - (factyp + factym) * zintp[aij(nx - 1, j)] + + factym * zintp[aij(nx - 1, j - 1)]; } -/* - call tridag (inf,diag,sup,f,res,nx) - do i=1,nx - zint(i,j)=res(i) - enddo -*/ + /* + call tridag (inf,diag,sup,f,res,nx) + do i=1,nx + zint(i,j)=res(i) + enddo + */ tridag(&inf, &diag, &sup, &f, &mut res, nx); for i in 0..nx { zint[aij(i, j)] = res[i]; } -/* - enddo - deallocate (f,diag,sup,inf,res) -*/ + /* + enddo + deallocate (f,diag,sup,inf,res) + */ } -/* - ! second pass along y-axis + /* + ! second pass along y-axis - allocate (f(ny),diag(ny),sup(ny),inf(ny),res(ny)) - f=0.d0 - diag=0.d0 - sup=0.d0 - inf=0.d0 - res=0.d0 - do i=2,nx-1 -*/ + allocate (f(ny),diag(ny),sup(ny),inf(ny),res(ny)) + f=0.d0 + diag=0.d0 + sup=0.d0 + inf=0.d0 + res=0.d0 + do i=2,nx-1 + */ f = vec![0.0; ny]; diag = vec![0.0; ny]; sup = vec![0.0; ny]; inf = vec![0.0; ny]; res = vec![0.0; ny]; - for i in 1..nx-1 { -/* - do j=2,ny-1 - factxp=(kdint(i+1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2 - factxm=(kdint(i-1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2 - factyp=(kdint(i,j+1)+kdint(i,j))/2.d0*(dt/2.)/dy**2 - factym=(kdint(i,j-1)+kdint(i,j))/2.d0*(dt/2.)/dy**2 - diag(j)=1.d0+factyp+factym - sup(j)=-factyp - inf(j)=-factym - f(j)=zint(i,j)+factxp*zint(i+1,j)-(factxp+factxm)*zint(i,j)+factxm*zint(i-1,j) - enddo -*/ - for j in 1..ny-1 { - factxp = (kdint[aij(i+1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx); - factxm = (kdint[aij(i-1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx); - factyp = (kdint[aij(i, j+1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy); - factym = (kdint[aij(i, j-1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy); + for i in 1..nx - 1 { + /* + do j=2,ny-1 + factxp=(kdint(i+1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2 + factxm=(kdint(i-1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2 + factyp=(kdint(i,j+1)+kdint(i,j))/2.d0*(dt/2.)/dy**2 + factym=(kdint(i,j-1)+kdint(i,j))/2.d0*(dt/2.)/dy**2 + diag(j)=1.d0+factyp+factym + sup(j)=-factyp + inf(j)=-factym + f(j)=zint(i,j)+factxp*zint(i+1,j)-(factxp+factxm)*zint(i,j)+factxm*zint(i-1,j) + enddo + */ + for j in 1..ny - 1 { + factxp = (kdint[aij(i + 1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx); + factxm = (kdint[aij(i - 1, j)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx); + factyp = (kdint[aij(i, j + 1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy); + factym = (kdint[aij(i, j - 1)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dy * dy); diag[j] = 1.0 + factyp + factym; sup[j] = -factyp; inf[j] = -factym; - f[j] = - zint[aij(i, j)] + factxp * zint[aij(i+1, j)] - - (factxp + factxm) * - zint[aij(i, j)] + factxm * zint[aij(i-1, j)]; + f[j] = zint[aij(i, j)] + factxp * zint[aij(i + 1, j)] + - (factxp + factxm) * zint[aij(i, j)] + + factxm * zint[aij(i - 1, j)]; } -/* -! bottom bc - if (cbc(1:1).eq.'1') then - diag(1)=1. - sup(1)=0. - f(1)=zint(i,1) - else - factxp=(kdint(i+1,1)+kdint(i,j))/2.d0*(dt/2.)/dx**2 - factxm=(kdint(i-1,1)+kdint(i,1))/2.d0*(dt/2.)/dx**2 - factyp=(kdint(i,2)+kdint(i,1))/2.d0*(dt/2.)/dy**2 - diag(1)=1.d0+factyp - sup(1)=-factyp - f(1)=zint(i,1)+factxp*zint(i+1,1)-(factxp+factxm)*zint(i,1)+factxm*zint(i-1,1) - endif -*/ + /* + ! bottom bc + if (cbc(1:1).eq.'1') then + diag(1)=1. + sup(1)=0. + f(1)=zint(i,1) + else + factxp=(kdint(i+1,1)+kdint(i,j))/2.d0*(dt/2.)/dx**2 + factxm=(kdint(i-1,1)+kdint(i,1))/2.d0*(dt/2.)/dx**2 + factyp=(kdint(i,2)+kdint(i,1))/2.d0*(dt/2.)/dy**2 + diag(1)=1.d0+factyp + sup(1)=-factyp + f(1)=zint(i,1)+factxp*zint(i+1,1)-(factxp+factxm)*zint(i,1)+factxm*zint(i-1,1) + endif + */ if true { diag[0] = 1.0; sup[0] = 0.0; @@ -297,76 +303,76 @@ pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (), // TODO: Check whether this j was actually supposed to be a 0 in the original // (probably). // factxp = (kdint[aij(i+1, 0)] + kdint[aij(i, j)]) / 2.0 * (dt / 2.0) / (dx * dx); - factxp = (kdint[aij(i+1, 0)] + kdint[aij(i, 0)]) / 2.0 * (dt / 2.0) / (dx * dx); - factxm = (kdint[aij(i-1, 0)] + kdint[aij(i, 0)]) / 2.0 * (dt / 2.0) / (dx * dx); + factxp = (kdint[aij(i + 1, 0)] + kdint[aij(i, 0)]) / 2.0 * (dt / 2.0) / (dx * dx); + factxm = (kdint[aij(i - 1, 0)] + kdint[aij(i, 0)]) / 2.0 * (dt / 2.0) / (dx * dx); factyp = (kdint[aij(i, 1)] + kdint[aij(i, 0)]) / 2.0 * (dt / 2.0) / (dy * dy); diag[0] = 1.0 + factyp; sup[0] = -factyp; - f[0] = - zint[aij(i, 0)] + factxp * zint[aij(i+1, 0)] - - (factxp + factxm) * - zint[aij(i, 0)] + factxm * zint[aij(i-1, 0)]; + f[0] = zint[aij(i, 0)] + factxp * zint[aij(i + 1, 0)] + - (factxp + factxm) * zint[aij(i, 0)] + + factxm * zint[aij(i - 1, 0)]; } -/* -! top bc - if (cbc(3:3).eq.'1') then - diag(ny)=1. - inf(ny)=0. - f(ny)=zint(i,ny) - else - factxp=(kdint(i+1,ny)+kdint(i,ny))/2.d0*(dt/2.)/dx**2 - factxm=(kdint(i-1,ny)+kdint(i,ny))/2.d0*(dt/2.)/dx**2 - factym=(kdint(i,ny-1)+kdint(i,ny))/2.d0*(dt/2.)/dy**2 - diag(ny)=1.d0+factym - inf(ny)=-factym - f(ny)=zint(i,ny)+factxp*zint(i+1,ny)-(factxp+factxm)*zint(i,ny)+factxm*zint(i-1,ny) - endif -*/ + /* + ! top bc + if (cbc(3:3).eq.'1') then + diag(ny)=1. + inf(ny)=0. + f(ny)=zint(i,ny) + else + factxp=(kdint(i+1,ny)+kdint(i,ny))/2.d0*(dt/2.)/dx**2 + factxm=(kdint(i-1,ny)+kdint(i,ny))/2.d0*(dt/2.)/dx**2 + factym=(kdint(i,ny-1)+kdint(i,ny))/2.d0*(dt/2.)/dy**2 + diag(ny)=1.d0+factym + inf(ny)=-factym + f(ny)=zint(i,ny)+factxp*zint(i+1,ny)-(factxp+factxm)*zint(i,ny)+factxm*zint(i-1,ny) + endif + */ if true { diag[ny - 1] = 1.0; inf[ny - 1] = 0.0; f[ny - 1] = zint[aij(i, ny - 1)]; } else { // reflective boundary - factxp = (kdint[aij(i+1, ny-1)] + kdint[aij(i, ny-1)]) / 2.0 * (dt / 2.0) / (dx * dx); - factxm = (kdint[aij(i-1, ny-1)] + kdint[aij(i, ny-1)]) / 2.0 * (dt / 2.0) / (dx * dx); - factym = (kdint[aij(i, ny-2)] + kdint[aij(i, ny-1)]) / 2.0 * (dt / 2.0) / (dy * dy); + factxp = + (kdint[aij(i + 1, ny - 1)] + kdint[aij(i, ny - 1)]) / 2.0 * (dt / 2.0) / (dx * dx); + factxm = + (kdint[aij(i - 1, ny - 1)] + kdint[aij(i, ny - 1)]) / 2.0 * (dt / 2.0) / (dx * dx); + factym = (kdint[aij(i, ny - 2)] + kdint[aij(i, ny - 1)]) / 2.0 * (dt / 2.0) / (dy * dy); diag[ny - 1] = 1.0 + factym; inf[ny - 1] = -factym; - f[ny - 1] = - zint[aij(i, ny-1)] + factxp * zint[aij(i+1, ny-1)] - - (factxp + factxm) * - zint[aij(i, ny-1)] + factxm * zint[aij(i-1, ny-1)]; + f[ny - 1] = zint[aij(i, ny - 1)] + factxp * zint[aij(i + 1, ny - 1)] + - (factxp + factxm) * zint[aij(i, ny - 1)] + + factxm * zint[aij(i - 1, ny - 1)]; } -/* - call tridag (inf,diag,sup,f,res,ny) - do j=1,ny - zintp(i,j)=res(j) - enddo -*/ + /* + call tridag (inf,diag,sup,f,res,ny) + do j=1,ny + zintp(i,j)=res(j) + enddo + */ tridag(&inf, &diag, &sup, &f, &mut res, ny); for j in 0..ny { zintp[aij(i, j)] = res[j]; } -/* - enddo - deallocate (f,diag,sup,inf,res) -*/ + /* + enddo + deallocate (f,diag,sup,inf,res) + */ } -/* - ! stores result in 1D array + /* + ! stores result in 1D array - do j=1,ny - do i=1,nx - ij=(j-1)*nx+i - etot(ij)=etot(ij)+h(ij)-zintp(i,j) - erate(ij)=erate(ij)+(h(ij)-zintp(i,j))/dt - h(ij)=zintp(i,j) - enddo - enddo + do j=1,ny + do i=1,nx + ij=(j-1)*nx+i + etot(ij)=etot(ij)+h(ij)-zintp(i,j) + erate(ij)=erate(ij)+(h(ij)-zintp(i,j))/dt + h(ij)=zintp(i,j) + enddo + enddo - b=min(h,b) -*/ + b=min(h,b) + */ for j in 0..ny { for i in 0..nx { ij = aij(i, j); @@ -378,13 +384,13 @@ pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (), b.par_iter_mut().zip(h).for_each(|(mut b, h)| { *b = h.min(*b); }); -/* - deallocate (zint,kdint,zintp) + /* + deallocate (zint,kdint,zintp) - return + return -end subroutine Diffusion -*/ + end subroutine Diffusion + */ } /* @@ -400,39 +406,39 @@ end subroutine Diffusion double precision a(n),b(n),c(n),r(n),u(n) */ pub fn tridag(a: &[f64], b: &[f64], c: &[f64], r: &[f64], u: &mut [f64], n: usize) { -/* - INTEGER j - double precision bet - double precision,dimension(:),allocatable::gam + /* + INTEGER j + double precision bet + double precision,dimension(:),allocatable::gam - allocate (gam(n)) + allocate (gam(n)) - if(b(1).eq.0.d0) stop 'in tridag' -*/ - let mut j : usize; - let mut bet : f64; - let mut precision : f64; + if(b(1).eq.0.d0) stop 'in tridag' + */ + let mut j: usize; + let mut bet: f64; + let mut precision: f64; let mut gam: Vec; gam = vec![Default::default(); n]; assert!(b[0] != 0.0); -/* + /* -! first pass + ! first pass -bet=b(1) -u(1)=r(1)/bet -do 11 j=2,n - gam(j)=c(j-1)/bet - bet=b(j)-a(j)*gam(j) - if(bet.eq.0.) then - print*,'tridag failed' - stop - endif - u(j)=(r(j)-a(j)*u(j-1))/bet - 11 continue -*/ + bet=b(1) + u(1)=r(1)/bet + do 11 j=2,n + gam(j)=c(j-1)/bet + bet=b(j)-a(j)*gam(j) + if(bet.eq.0.) then + print*,'tridag failed' + stop + endif + u(j)=(r(j)-a(j)*u(j-1))/bet + 11 continue + */ bet = b[0]; u[0] = r[0] / bet; for j in 1..n { @@ -448,21 +454,21 @@ do 11 j=2,n // = r'[j] / b'[j] u[j] = (r[j] - a[j] * u[j - 1]) / bet; } -/* - ! second pass + /* + ! second pass - do 12 j=n-1,1,-1 - u(j)=u(j)-gam(j+1)*u(j+1) - 12 continue -*/ + do 12 j=n-1,1,-1 + u(j)=u(j)-gam(j+1)*u(j+1) + 12 continue + */ for j in (0..n - 1).rev() { u[j] = u[j] - gam[j + 1] * u[j + 1]; } -/* - deallocate (gam) + /* + deallocate (gam) - return + return - END -*/ + END + */ } diff --git a/world/src/sim/mod.rs b/world/src/sim/mod.rs index 9fddd9bd25..b8fb96acff 100644 --- a/world/src/sim/mod.rs +++ b/world/src/sim/mod.rs @@ -12,9 +12,9 @@ pub use self::erosion::{ pub use self::location::Location; pub use self::settlement::Settlement; pub use self::util::{ - cdf_irwin_hall, downhill, get_oceans, HybridMulti as HybridMulti_, local_cells, map_edge_factor, neighbors, - NEIGHBOR_DELTA, ScaleBias, - uniform_idx_as_vec2, uniform_noise, uphill, vec2_as_uniform_idx, InverseCdf, + cdf_irwin_hall, downhill, get_oceans, local_cells, map_edge_factor, neighbors, + uniform_idx_as_vec2, uniform_noise, uphill, vec2_as_uniform_idx, HybridMulti as HybridMulti_, + InverseCdf, ScaleBias, NEIGHBOR_DELTA, }; use crate::{ @@ -29,8 +29,8 @@ use common::{ vol::RectVolSize, }; use noise::{ - BasicMulti, Billow, Fbm, HybridMulti, MultiFractal, NoiseFn, RangeFunction, - RidgedMulti, Seedable, SuperSimplex, Worley, + BasicMulti, Billow, Fbm, HybridMulti, MultiFractal, NoiseFn, RangeFunction, RidgedMulti, + Seedable, SuperSimplex, Worley, }; use num::{Float, Signed}; use rand::{Rng, SeedableRng}; @@ -39,9 +39,9 @@ use rayon::prelude::*; use serde_derive::{Deserialize, Serialize}; use std::{ collections::HashMap, - io::{BufReader, BufWriter}, f32, f64, fs::File, + io::{BufReader, BufWriter}, ops::{Add, Div, Mul, Neg, Sub}, path::PathBuf, sync::Arc, @@ -54,10 +54,7 @@ use vek::*; // cleanly representable in f32 (that stops around 1024 * 4 * 1024 * 4, for signed floats anyway) // but I think that is probably less important since I don't think we actually cast a chunk id to // float, just coordinates... could be wrong though! -pub const WORLD_SIZE: Vec2 = Vec2 { - x: 1024, - y: 1024, -}; +pub const WORLD_SIZE: Vec2 = Vec2 { x: 1024, y: 1024 }; /// A structure that holds cached noise values and cumulative distribution functions for the input /// that led to those values. See the definition of InverseCdf for a description of how to @@ -107,7 +104,7 @@ pub(crate) struct GenCtx { pub town_gen: StructureGen2d, pub river_seed: RandomField, - pub rock_strength_nz: /*HybridMulti_*/Fbm, + pub rock_strength_nz: Fbm, pub uplift_nz: Worley, } @@ -146,7 +143,7 @@ impl Default for WorldOpts { /// A way to store certain components between runs of map generation. Only intended for /// development purposes--no attempt is made to detect map invalidation or make sure that the map /// is synchronized with updates to noise-rs, changes to other parameters, etc. -#[derive(Serialize,Deserialize)] +#[derive(Serialize, Deserialize)] pub struct WorldFile { /// Saved altitude height map. pub alt: Box<[Alt]>, @@ -166,7 +163,9 @@ pub struct WorldSim { impl WorldSim { pub fn generate(seed: u32, opts: WorldOpts) -> Self { let mut rng = ChaChaRng::from_seed(seed_expan::rng_state(seed)); - let continent_scale = 5_000.0f64/*32768.0*/.div(32.0).mul(TerrainChunkSize::RECT_SIZE.x as f64); + let continent_scale = 5_000.0f64 /*32768.0*/ + .div(32.0) + .mul(TerrainChunkSize::RECT_SIZE.x as f64); let rock_lacunarity = /*0.5*/2.0/*HybridMulti::DEFAULT_LACUNARITY*/; let uplift_scale = /*512.0*//*256.0*/128.0; let uplift_turb_scale = uplift_scale / 4.0/*32.0*//*64.0*/; @@ -290,9 +289,9 @@ impl WorldSim { let erosion_pow_low = /*0.25*//*1.5*//*2.0*//*0.5*//*4.0*//*0.25*//*1.0*//*2.0*//*1.5*//*1.5*//*0.35*//*0.43*//*0.5*//*0.45*//*0.37*/1.002; let erosion_pow_high = /*1.5*//*1.0*//*0.55*//*0.51*//*2.0*/1.002; let erosion_center = /*0.45*//*0.75*//*0.75*//*0.5*//*0.75*/0.5; - let n_steps = /*200*//*10_000*//*1000*//*50*//*100*/100;//100; // /*100*//*50*//*100*//*100*//*50*//*25*/25/*100*//*37*/;//150;//37/*100*/;//50;//50;//37;//50;//37; // /*37*//*29*//*40*//*150*/37; //150;//200; - let n_small_steps = 0;//25;//8;//50;//50;//8;//8;//8;//8;//8; // 8 - let n_post_load_steps = 0;//25;//8 + let n_steps = /*200*//*10_000*//*1000*//*50*//*100*/100; //100; // /*100*//*50*//*100*//*100*//*50*//*25*/25/*100*//*37*/;//150;//37/*100*/;//50;//50;//37;//50;//37; // /*37*//*29*//*40*//*150*/37; //150;//200; + let n_small_steps = 0; //25;//8;//50;//50;//8;//8;//8;//8;//8; // 8 + let n_post_load_steps = 0; //25;//8 // Logistic regression. Make sure x ∈ (0, 1). let logit = |x: f64| x.ln() - (-x).ln_1p(); @@ -304,7 +303,8 @@ impl WorldSim { let exp_inverse_cdf = |x: f64/*, pow: f64*/| -(-x).ln_1p()/* / ln(pow)*/; // 2 / pi * ln(tan(pi/2 * p)) - let hypsec_inverse_cdf = |x: f64| f64::consts::FRAC_2_PI * ((x * f64::consts::FRAC_PI_2).tan().ln()); + let hypsec_inverse_cdf = + |x: f64| f64::consts::FRAC_2_PI * ((x * f64::consts::FRAC_PI_2).tan().ln()); let min_epsilon = 1.0 / (WORLD_SIZE.x as f64 * WORLD_SIZE.y as f64).max(f64::EPSILON as f64 * 0.5); @@ -364,10 +364,9 @@ impl WorldSim { .alt_nz .get((wposf.div(10_000.0)).into_array()) .min(1.0) - .max(-1.0) - /* .mul(0.25) - .add(0.125) */) - // .add(0.5) + .max(-1.0)/* .mul(0.25) + .add(0.125) */) + // .add(0.5) .sub(0.05) // .add(0.05) // .add(0.075) @@ -525,102 +524,110 @@ impl WorldSim { }); // Calculate oceans. - let old_height = |posi: usize| alt_old[posi].1 * CONFIG.mountain_scale * height_scale as f32; + let old_height = + |posi: usize| alt_old[posi].1 * CONFIG.mountain_scale * height_scale as f32; /* let is_ocean = (0..WORLD_SIZE.x * WORLD_SIZE.y) - .into_par_iter() - .map(|i| map_edge_factor(i) == 0.0) - .collect::>(); */ + .into_par_iter() + .map(|i| map_edge_factor(i) == 0.0) + .collect::>(); */ let is_ocean = get_oceans(old_height); let is_ocean_fn = |posi: usize| is_ocean[posi]; let turb_wposf_div = 8.0/*64.0*/; - let uplift_nz_dist = gen_ctx.uplift_nz - .clone() - .enable_range(true); + let uplift_nz_dist = gen_ctx.uplift_nz.clone().enable_range(true); // Recalculate altitudes without oceans. // NaNs in these uniform vectors wherever is_ocean_fn returns true. - let (alt_old_no_ocean, alt_old_inverse) = - uniform_noise(|posi, _| { - if is_ocean_fn(posi) { - None - } else { - Some(old_height(posi) /*.abs()*/) - } - }); - let (uplift_uniform, _) = - uniform_noise(|posi, wposf| { - if is_ocean_fn(posi) { - None - } else { - let turb_wposf = - wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(turb_wposf_div); - let turb = Vec2::new( - gen_ctx.turb_x_nz.get(turb_wposf.into_array()), - gen_ctx.turb_y_nz.get(turb_wposf.into_array()), - ) * uplift_turb_scale * TerrainChunkSize::RECT_SIZE.map(|e| e as f64); - // let turb = Vec2::zero(); - let turb_wposf = wposf + turb; - let turb_wposi = turb_wposf - .map2(TerrainChunkSize::RECT_SIZE, |e, f| e / f as f64) - .map2(WORLD_SIZE, |e, f| (e as i32).max(f as i32 - 1).min(0)); - let turb_posi = vec2_as_uniform_idx(turb_wposi); - let udist = uplift_nz_dist.get(turb_wposf.into_array()) - .min(1.0) - .max(-1.0) - .mul(0.5) - .add(0.5); - let uheight = gen_ctx.uplift_nz.get(turb_wposf.into_array()) - /* .min(0.5) - .max(-0.5)*/ - .min(1.0) - .max(-1.0) - .mul(0.5) - .add(0.5); - let uchaos = /* gen_ctx.chaos_nz.get((wposf.div(3_000.0)).into_array()) + let (alt_old_no_ocean, alt_old_inverse) = uniform_noise(|posi, _| { + if is_ocean_fn(posi) { + None + } else { + Some(old_height(posi) /*.abs()*/) + } + }); + let (uplift_uniform, _) = uniform_noise(|posi, wposf| { + if is_ocean_fn(posi) { + None + } else { + let turb_wposf = wposf + .div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)) + .div(turb_wposf_div); + let turb = Vec2::new( + gen_ctx.turb_x_nz.get(turb_wposf.into_array()), + gen_ctx.turb_y_nz.get(turb_wposf.into_array()), + ) * uplift_turb_scale + * TerrainChunkSize::RECT_SIZE.map(|e| e as f64); + // let turb = Vec2::zero(); + let turb_wposf = wposf + turb; + let turb_wposi = turb_wposf + .map2(TerrainChunkSize::RECT_SIZE, |e, f| e / f as f64) + .map2(WORLD_SIZE, |e, f| (e as i32).max(f as i32 - 1).min(0)); + let turb_posi = vec2_as_uniform_idx(turb_wposi); + let udist = uplift_nz_dist + .get(turb_wposf.into_array()) + .min(1.0) + .max(-1.0) + .mul(0.5) + .add(0.5); + let uheight = gen_ctx + .uplift_nz + .get(turb_wposf.into_array()) + /* .min(0.5) + .max(-0.5)*/ + .min(1.0) + .max(-1.0) + .mul(0.5) + .add(0.5); + let uchaos = /* gen_ctx.chaos_nz.get((wposf.div(3_000.0)).into_array()) .min(1.0) .max(-1.0) .mul(0.5) .add(0.5); */ chaos[posi].1; - let uchaos_1 = (uchaos as f64) / 1.32; + let uchaos_1 = (uchaos as f64) / 1.32; - let oheight = /*alt_old*//*alt_base*/alt_old_no_ocean[/*(turb_posi / 64) * 64*/posi].0 as f64 - 0.5; - assert!(udist >= 0.0); - assert!(udist <= 1.0); - let uheight_1 = uheight;//.powf(2.0); - let udist_1 = (0.5 - udist).mul(2.0).max(0.0); - let udist_2 = udist.mul(2.0).min(1.0); - let udist_3 = (1.0 - udist).max(0.0); - let udist_4 = udist.min(1.0); - let variation = 1.0.min(64.0 * 64.0 / (WORLD_SIZE.x as f64 * WORLD_SIZE.y as f64 * (TerrainChunkSize::RECT_SIZE.x as f64 * TerrainChunkSize::RECT_SIZE.y as f64 / 128.0 / 128.0))); - let variation_1 = (uheight * /*udist_2*/udist_4).min(variation); - let height = - (oheight + 0.5).powf(2.0); - // 1.0 - variation + variation * uchaos_1; - // uheight * /*udist_2*/udist_4 - variation_1 + variation_1 * uchaos_1; - // uheight * (0.5 + 0.5 * ((uchaos as f64) / 1.32)) - 0.125; - // 0.2; - // 1.0; - // uheight_1; - // uheight_1 * (0.8 + 0.2 * oheight.signum() * oheight.abs().powf(0.25)); - // uheight_1 * (/*udist_2*/udist.powf(2.0) * (f64::consts::PI * 2.0 * (1.0 / (1.0 - udist).max(f64::EPSILON)).min(2.5)/*udist * 5.0*/ * 2.0).cos().mul(0.5).add(0.5)); - // uheight * udist_ * (udist_ * 4.0 * 2 * f64::consts::PI).sin() - // uheight; - // (0.8 * uheight + oheight.powf(2.0) * 0.2).max(0.0).min(1.0); - // ((0.8 - 0.2) * uheight + 0.2 + oheight.signum() * oheight.abs().powf(/*0.5*/2.0) * udist_2.powf(2.0)).max(0.0).min(1.0); - // ((0.8 - 0.2) * uheight + 0.2 + oheight.signum() * oheight.abs().powf(/*0.5*/2.0) * 0.2).max(0.0).min(1.0); - // (0.8 * uheight * udist_1 + 0.8 * udist_2 + oheight.powf(2.0) * 0.2).max(0.0).min(1.0); - /* uheight * 0.8 * udist_1.powf(2.0) + - /*exp_inverse_cdf*/(oheight/*.max(0.0).min(max_epsilon).abs()*/).powf(2.0) * 0.2 * udist_2.powf(2.0); */ - // (uheight + oheight.powf(2.0) * 0.05).max(0.0).min(1.0); - // (uheight + oheight.powf(2.0) * 0.2).max(0.0).min(1.0); - // * (1.0 - udist);// uheight * (1.0 - udist)/*oheight*//* * udist*/ + oheight * udist;/*uheight * (1.0 - udist);*/ - // let height = uheight * (0.5 - udist) * 0.8 + (oheight.signum() * oheight.max(0.0).abs().powf(2.0)) * 0.2;// * (1.0 - udist);// uheight * (1.0 - udist)/*oheight*//* * udist*/ + oheight * udist;/*uheight * (1.0 - udist);*/ - - Some(height) - } - }); + let oheight = /*alt_old*//*alt_base*/alt_old_no_ocean[/*(turb_posi / 64) * 64*/posi].0 as f64 - 0.5; + assert!(udist >= 0.0); + assert!(udist <= 1.0); + let uheight_1 = uheight; //.powf(2.0); + let udist_1 = (0.5 - udist).mul(2.0).max(0.0); + let udist_2 = udist.mul(2.0).min(1.0); + let udist_3 = (1.0 - udist).max(0.0); + let udist_4 = udist.min(1.0); + let variation = 1.0.min( + 64.0 * 64.0 + / (WORLD_SIZE.x as f64 + * WORLD_SIZE.y as f64 + * (TerrainChunkSize::RECT_SIZE.x as f64 + * TerrainChunkSize::RECT_SIZE.y as f64 + / 128.0 + / 128.0)), + ); + let variation_1 = (uheight * /*udist_2*/udist_4).min(variation); + let height = (oheight + 0.5).powf(2.0); + // 1.0 - variation + variation * uchaos_1; + // uheight * /*udist_2*/udist_4 - variation_1 + variation_1 * uchaos_1; + // uheight * (0.5 + 0.5 * ((uchaos as f64) / 1.32)) - 0.125; + // 0.2; + // 1.0; + // uheight_1; + // uheight_1 * (0.8 + 0.2 * oheight.signum() * oheight.abs().powf(0.25)); + // uheight_1 * (/*udist_2*/udist.powf(2.0) * (f64::consts::PI * 2.0 * (1.0 / (1.0 - udist).max(f64::EPSILON)).min(2.5)/*udist * 5.0*/ * 2.0).cos().mul(0.5).add(0.5)); + // uheight * udist_ * (udist_ * 4.0 * 2 * f64::consts::PI).sin() + // uheight; + // (0.8 * uheight + oheight.powf(2.0) * 0.2).max(0.0).min(1.0); + // ((0.8 - 0.2) * uheight + 0.2 + oheight.signum() * oheight.abs().powf(/*0.5*/2.0) * udist_2.powf(2.0)).max(0.0).min(1.0); + // ((0.8 - 0.2) * uheight + 0.2 + oheight.signum() * oheight.abs().powf(/*0.5*/2.0) * 0.2).max(0.0).min(1.0); + // (0.8 * uheight * udist_1 + 0.8 * udist_2 + oheight.powf(2.0) * 0.2).max(0.0).min(1.0); + /* uheight * 0.8 * udist_1.powf(2.0) + + /*exp_inverse_cdf*/(oheight/*.max(0.0).min(max_epsilon).abs()*/).powf(2.0) * 0.2 * udist_2.powf(2.0); */ + // (uheight + oheight.powf(2.0) * 0.05).max(0.0).min(1.0); + // (uheight + oheight.powf(2.0) * 0.2).max(0.0).min(1.0); + // * (1.0 - udist);// uheight * (1.0 - udist)/*oheight*//* * udist*/ + oheight * udist;/*uheight * (1.0 - udist);*/ + // let height = uheight * (0.5 - udist) * 0.8 + (oheight.signum() * oheight.max(0.0).abs().powf(2.0)) * 0.2;// * (1.0 - udist);// uheight * (1.0 - udist)/*oheight*//* * udist*/ + oheight * udist;/*uheight * (1.0 - udist);*/ + Some(height) + } + }); let old_height_uniform = |posi: usize| alt_old_no_ocean[posi].0; let alt_old_min_uniform = 0.0; @@ -684,22 +691,25 @@ impl WorldSim { if is_ocean_fn(posi) { return 1.0; } - let wposf = (uniform_idx_as_vec2(posi) - * TerrainChunkSize::RECT_SIZE.map(|e| e as i32)) - .map(|e| e as f64); - let turb_wposf = - wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(turb_wposf_div); + let wposf = (uniform_idx_as_vec2(posi) * TerrainChunkSize::RECT_SIZE.map(|e| e as i32)) + .map(|e| e as f64); + let turb_wposf = wposf + .div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)) + .div(turb_wposf_div); let turb = Vec2::new( gen_ctx.turb_x_nz.get(turb_wposf.into_array()), gen_ctx.turb_y_nz.get(turb_wposf.into_array()), - ) * uplift_turb_scale * TerrainChunkSize::RECT_SIZE.map(|e| e as f64); + ) * uplift_turb_scale + * TerrainChunkSize::RECT_SIZE.map(|e| e as f64); // let turb = Vec2::zero(); let turb_wposf = wposf + turb; let turb_wposi = turb_wposf .map2(TerrainChunkSize::RECT_SIZE, |e, f| e / f as f64) .map2(WORLD_SIZE, |e, f| (e as i32).max(f as i32 - 1).min(0)); let turb_posi = vec2_as_uniform_idx(turb_wposi); - let uheight = gen_ctx.uplift_nz.get(turb_wposf.into_array()) + let uheight = gen_ctx + .uplift_nz + .get(turb_wposf.into_array()) /* .min(0.5) .max(-0.5)*/ .min(1.0) @@ -716,9 +726,7 @@ impl WorldSim { // ((3.5 - 1.5) * (1.0 - uheight) + 1.5) as f32 1.0 }; - let theta_func = |posi| { - 0.4 - }; + let theta_func = |posi| 0.4; let kf_func = { |posi| { if is_ocean_fn(posi) { @@ -731,19 +739,23 @@ impl WorldSim { let wposf = (uniform_idx_as_vec2(posi) * TerrainChunkSize::RECT_SIZE.map(|e| e as i32)) .map(|e| e as f64); - let turb_wposf = - wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(turb_wposf_div); + let turb_wposf = wposf + .div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)) + .div(turb_wposf_div); let turb = Vec2::new( gen_ctx.turb_x_nz.get(turb_wposf.into_array()), gen_ctx.turb_y_nz.get(turb_wposf.into_array()), - ) * uplift_turb_scale * TerrainChunkSize::RECT_SIZE.map(|e| e as f64); + ) * uplift_turb_scale + * TerrainChunkSize::RECT_SIZE.map(|e| e as f64); // let turb = Vec2::zero(); let turb_wposf = wposf + turb; let turb_wposi = turb_wposf .map2(TerrainChunkSize::RECT_SIZE, |e, f| e / f as f64) .map2(WORLD_SIZE, |e, f| (e as i32).max(f as i32 - 1).min(0)); let turb_posi = vec2_as_uniform_idx(turb_wposi); - let uheight = gen_ctx.uplift_nz.get(turb_wposf.into_array()) + let uheight = gen_ctx + .uplift_nz + .get(turb_wposf.into_array()) /* .min(0.5) .max(-0.5)*/ .min(1.0) @@ -784,7 +796,7 @@ impl WorldSim { // ((1.0 - uheight) * (0.5 - 0.5 * /*((1.32 - uchaos as f64) / 1.32)*/oheight_2) * (1.5e-4 - 2.0e-6) + 2.0e-6) // ((1.0 - uheight) * (0.5 - 0.5 * /*((1.32 - uchaos as f64) / 1.32)*/oheight) * (1.5e-4 - 2.0e-6) + 2.0e-6) // 2e-5 - 2.5e-6 * 16.0.powf(0.4)/* / 4.0 * 0.25 *//* * 4.0*/ + 2.5e-6 * 16.0.powf(0.4) /* / 4.0 * 0.25 *//* * 4.0*/ // 2.9e-10 // ((1.0 - uheight) * (5e-5 - 2.9e-10) + 2.9e-10) // ((1.0 - uheight) * (5e-5 - 2.9e-14) + 2.9e-14) @@ -798,19 +810,23 @@ impl WorldSim { let wposf = (uniform_idx_as_vec2(posi) * TerrainChunkSize::RECT_SIZE.map(|e| e as i32)) .map(|e| e as f64); - let turb_wposf = - wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(turb_wposf_div); + let turb_wposf = wposf + .div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)) + .div(turb_wposf_div); let turb = Vec2::new( gen_ctx.turb_x_nz.get(turb_wposf.into_array()), gen_ctx.turb_y_nz.get(turb_wposf.into_array()), - ) * uplift_turb_scale * TerrainChunkSize::RECT_SIZE.map(|e| e as f64); + ) * uplift_turb_scale + * TerrainChunkSize::RECT_SIZE.map(|e| e as f64); // let turb = Vec2::zero(); let turb_wposf = wposf + turb; let turb_wposi = turb_wposf .map2(TerrainChunkSize::RECT_SIZE, |e, f| e / f as f64) .map2(WORLD_SIZE, |e, f| (e as i32).max(f as i32 - 1).min(0)); let turb_posi = vec2_as_uniform_idx(turb_wposi); - let uheight = gen_ctx.uplift_nz.get(turb_wposf.into_array()) + let uheight = gen_ctx + .uplift_nz + .get(turb_wposf.into_array()) /* .min(0.5) .max(-0.5)*/ .min(1.0) @@ -821,209 +837,215 @@ impl WorldSim { // kd = 1.5e-2: normal-high (plateau [fan sediment]) // kd = 1e-2: normal (plateau) 1.0e-2 * (1.0 / 16.0) // m_old^2 / y * (1 m_new / 4 m_old)^2 - // (uheight * (1.0e-1 - 1.0e-2) + 1.0e-2) + // (uheight * (1.0e-1 - 1.0e-2) + 1.0e-2) } }; - let g_func = - |posi| { - if /*is_ocean_fn(posi)*/map_edge_factor(posi) == 0.0 { - return 0.0; - // return 5.0; - } - let wposf = (uniform_idx_as_vec2(posi) - * TerrainChunkSize::RECT_SIZE.map(|e| e as i32)) + let g_func = |posi| { + if + /*is_ocean_fn(posi)*/ + map_edge_factor(posi) == 0.0 { + return 0.0; + // return 5.0; + } + let wposf = (uniform_idx_as_vec2(posi) * TerrainChunkSize::RECT_SIZE.map(|e| e as i32)) .map(|e| e as f64); - let turb_wposf = - wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(turb_wposf_div); - let turb = Vec2::new( - gen_ctx.turb_x_nz.get(turb_wposf.into_array()), - gen_ctx.turb_y_nz.get(turb_wposf.into_array()), - ) * uplift_turb_scale * TerrainChunkSize::RECT_SIZE.map(|e| e as f64); - // let turb = Vec2::zero(); - let turb_wposf = wposf + turb; - let turb_wposi = turb_wposf - .map2(TerrainChunkSize::RECT_SIZE, |e, f| e / f as f64) - .map2(WORLD_SIZE, |e, f| (e as i32).max(f as i32 - 1).min(0)); - let turb_posi = vec2_as_uniform_idx(turb_wposi); - let uheight = gen_ctx.uplift_nz.get(turb_wposf.into_array()) - /* .min(0.5) - .max(-0.5)*/ - .min(1.0) - .max(-1.0) - .mul(0.5) - .add(0.5); + let turb_wposf = wposf + .div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)) + .div(turb_wposf_div); + let turb = Vec2::new( + gen_ctx.turb_x_nz.get(turb_wposf.into_array()), + gen_ctx.turb_y_nz.get(turb_wposf.into_array()), + ) * uplift_turb_scale + * TerrainChunkSize::RECT_SIZE.map(|e| e as f64); + // let turb = Vec2::zero(); + let turb_wposf = wposf + turb; + let turb_wposi = turb_wposf + .map2(TerrainChunkSize::RECT_SIZE, |e, f| e / f as f64) + .map2(WORLD_SIZE, |e, f| (e as i32).max(f as i32 - 1).min(0)); + let turb_posi = vec2_as_uniform_idx(turb_wposi); + let uheight = gen_ctx + .uplift_nz + .get(turb_wposf.into_array()) + /* .min(0.5) + .max(-0.5)*/ + .min(1.0) + .max(-1.0) + .mul(0.5) + .add(0.5); - let uchaos = /* gen_ctx.chaos_nz.get((wposf.div(3_000.0)).into_array()) + let uchaos = /* gen_ctx.chaos_nz.get((wposf.div(3_000.0)).into_array()) .min(1.0) .max(-1.0) .mul(0.5) .add(0.5); */ chaos[posi].1; - assert!(uchaos <= 1.32); + assert!(uchaos <= 1.32); - // G = d* v_s / p_0, where - // v_s is the settling velocity of sediment grains - // p_0 is the mean precipitation rate - // d* is the sediment concentration ratio (between concentration near riverbed - // interface, and average concentration over the water column). - // d* varies with Rouse number which defines relative contribution of bed, suspended, - // and washed loads. - // - // G is typically on the order of 1 or greater. However, we are only guaranteed to - // converge for G ≤ 1, so we keep it in the chaos range of [0.12, 1.32]. - // (((1.32 - uchaos) / 1.32).powf(0.75) * 1.32).min(/*1.1*/1.0) - // ((1.32 - 0.12) * (1.0 - uheight) + 0.12) as f32 - // 1.1 * (1.0 - uheight) as f32 - // 1.0 * (1.0 - uheight) as f32 - // 1.0 - // 5.0 - // 10.0 - // 2.0 - // 0.0 - 1.0 - // 1.5 - }; - let uplift_fn = - |posi| { - if is_ocean_fn(posi) { - /* return 1e-2 - .mul(max_erosion_per_delta_t) as f32; */ - return 0.0; - } - let wposf = (uniform_idx_as_vec2(posi) - * TerrainChunkSize::RECT_SIZE.map(|e| e as i32)) + // G = d* v_s / p_0, where + // v_s is the settling velocity of sediment grains + // p_0 is the mean precipitation rate + // d* is the sediment concentration ratio (between concentration near riverbed + // interface, and average concentration over the water column). + // d* varies with Rouse number which defines relative contribution of bed, suspended, + // and washed loads. + // + // G is typically on the order of 1 or greater. However, we are only guaranteed to + // converge for G ≤ 1, so we keep it in the chaos range of [0.12, 1.32]. + // (((1.32 - uchaos) / 1.32).powf(0.75) * 1.32).min(/*1.1*/1.0) + // ((1.32 - 0.12) * (1.0 - uheight) + 0.12) as f32 + // 1.1 * (1.0 - uheight) as f32 + // 1.0 * (1.0 - uheight) as f32 + // 1.0 + // 5.0 + // 10.0 + // 2.0 + // 0.0 + 1.0 + // 1.5 + }; + let uplift_fn = |posi| { + if is_ocean_fn(posi) { + /* return 1e-2 + .mul(max_erosion_per_delta_t) as f32; */ + return 0.0; + } + let wposf = (uniform_idx_as_vec2(posi) * TerrainChunkSize::RECT_SIZE.map(|e| e as i32)) .map(|e| e as f64); - let alt_main = { - // Extension upwards from the base. A positive number from 0 to 1 curved to be - // maximal at 0. Also to be multiplied by CONFIG.mountain_scale. - let alt_main = (gen_ctx - .alt_nz - .get((wposf.div(2_000.0)).into_array()) + let alt_main = { + // Extension upwards from the base. A positive number from 0 to 1 curved to be + // maximal at 0. Also to be multiplied by CONFIG.mountain_scale. + let alt_main = (gen_ctx + .alt_nz + .get((wposf.div(2_000.0)).into_array()) + .min(1.0) + .max(-1.0)) + .abs() + .powf(1.35); + + fn spring(x: f64, pow: f64) -> f64 { + x.abs().powf(pow) * x.signum() + } + + (0.0 + alt_main + + (gen_ctx + .small_nz + .get((wposf.div(300.0)).into_array()) .min(1.0) .max(-1.0)) - .abs() - .powf(1.35); - - fn spring(x: f64, pow: f64) -> f64 { - x.abs().powf(pow) * x.signum() - } - - (0.0 + alt_main - + (gen_ctx - .small_nz - .get((wposf.div(300.0)).into_array()) - .min(1.0) - .max(-1.0)) - .mul(alt_main.powf(0.8).max(/*0.25*/ 0.15)) - .mul(0.3) - .add(1.0) - .mul(0.4) - /* + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0) - .mul(0.045)*/) - }; - let height = + .mul(alt_main.powf(0.8).max(/*0.25*/ 0.15)) + .mul(0.3) + .add(1.0) + .mul(0.4)/* + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0) + .mul(0.045)*/) + }; + let height = ((/*old_height_uniform*/uplift_uniform[posi]./*0*/1 - alt_old_min_uniform) as f64 / (alt_old_max_uniform - alt_old_min_uniform) as f64) /*((old_height(posi) - alt_old_min) as f64 / (alt_old_max - alt_old_min) as f64)*/ ; - let height = height.mul(max_epsilon - min_epsilon).add(min_epsilon); - /*.max(1e-7 / CONFIG.mountain_scale as f64) - .min(1.0f64 - 1e-7);*/ - /* let alt_main = { - // Extension upwards from the base. A positive number from 0 to 1 curved to be - // maximal at 0. Also to be multiplied by CONFIG.mountain_scale. - let alt_main = (gen_ctx - .alt_nz - .get((wposf.div(2_000.0)).into_array()) + let height = height.mul(max_epsilon - min_epsilon).add(min_epsilon); + /*.max(1e-7 / CONFIG.mountain_scale as f64) + .min(1.0f64 - 1e-7);*/ + /* let alt_main = { + // Extension upwards from the base. A positive number from 0 to 1 curved to be + // maximal at 0. Also to be multiplied by CONFIG.mountain_scale. + let alt_main = (gen_ctx + .alt_nz + .get((wposf.div(2_000.0)).into_array()) + .min(1.0) + .max(-1.0)) + .abs() + .powf(1.35); + + fn spring(x: f64, pow: f64) -> f64 { + x.abs().powf(pow) * x.signum() + } + + (0.0 + alt_main + + (gen_ctx + .small_nz + .get((wposf.div(300.0)).into_array()) .min(1.0) .max(-1.0)) - .abs() - .powf(1.35); - - fn spring(x: f64, pow: f64) -> f64 { - x.abs().powf(pow) * x.signum() - } - - (0.0 + alt_main - + (gen_ctx - .small_nz - .get((wposf.div(300.0)).into_array()) - .min(1.0) - .max(-1.0)) - .mul(alt_main.powf(0.8).max(/*0.25*/ 0.15)) - .mul(0.3) - .add(1.0) - .mul(0.4) - + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0).mul(0.045)) - }; */ - // let height = height + (alt_main./*to_le_bytes()[7]*/to_bits() & 1) as f64 * ((1.0 / CONFIG.mountain_scale as f64).powf(1.0 / erosion_pow_low)); - let height = erosion_factor(height); - assert!(height >= 0.0); - assert!(height <= 1.0); - // assert!(alt_main >= 0.0); - let (bump_factor, bump_max) = if - /*height < f32::EPSILON as f64 * 0.5*//*false*/ - /*true*/false { - ( - /*(alt_main./*to_le_bytes()[7]*/to_bits() & 1) as f64*/ - (alt_main / CONFIG.mountain_scale as f64 * 128.0).mul(0.1).powf(1.2) * /*(1.0 / CONFIG.mountain_scale as f64)*/(f32::EPSILON * 0.5) as f64, - (f32::EPSILON * 0.5) as f64, - ) - } else { - (0.0, 0.0) - }; - // tan(6/360*2*pi)*32 ~ 3.4 - // 3.4/32*512 ~ 54 - // 18/32*512 ~ 288 - // tan(pi/6)*32 ~ 18 - // tan(54/360*2*pi)*32 - // let height = 1.0f64; - let turb_wposf = - wposf.div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(turb_wposf_div); - let turb = Vec2::new( - gen_ctx.turb_x_nz.get(turb_wposf.into_array()), - gen_ctx.turb_y_nz.get(turb_wposf.into_array()), - ) * uplift_turb_scale * TerrainChunkSize::RECT_SIZE.map(|e| e as f64); - let turb_wposf = wposf + turb; - let uheight = gen_ctx.uplift_nz.get(turb_wposf.into_array()) - /* .min(0.5) - .max(-0.5)*/ - .min(1.0) - .max(-1.0) - .mul(0.5) - .add(0.5); - // u = 1e-3: normal-high (dike, mountain) - // u = 5e-4: normal (mid example in Yuan, average mountain uplift) - // u = 2e-4: low (low example in Yuan; known that lagoons etc. may have u ~ 0.05). - // u = 0: low (plateau [fan, altitude = 0.0]) - // let height = uheight; - // let height = 1.0f64; - - // let height = 1.0 / 7.0f64; - // let height = 0.0 / 31.0f64; - let bfrac = /*erosion_factor(0.5);*/0.0; - let height = (height - bfrac).abs().div(1.0 - bfrac); - let height = height - /* .mul(31.0 / 32.0) - .add(1.0 / 32.0) */ - /* .mul(15.0 / 16.0) - .add(1.0 / 16.0) */ - /* .mul(5.0 / 8.0) - .add(3.0 / 8.0) */ - /* .mul(6.0 / 8.0) - .add(2.0 / 8.0) */ - /* .mul(7.0 / 8.0) - .add(1.0 / 8.0) */ - .mul(max_erosion_per_delta_t) - .sub(/*1.0 / CONFIG.mountain_scale as f64*/ bump_max) - .add(bump_factor); - /* .sub(/*1.0 / CONFIG.mountain_scale as f64*/(f32::EPSILON * 0.5) as f64) - .add(bump_factor); */ - height as f32 + .mul(alt_main.powf(0.8).max(/*0.25*/ 0.15)) + .mul(0.3) + .add(1.0) + .mul(0.4) + + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0).mul(0.045)) + }; */ + // let height = height + (alt_main./*to_le_bytes()[7]*/to_bits() & 1) as f64 * ((1.0 / CONFIG.mountain_scale as f64).powf(1.0 / erosion_pow_low)); + let height = erosion_factor(height); + assert!(height >= 0.0); + assert!(height <= 1.0); + // assert!(alt_main >= 0.0); + let (bump_factor, bump_max) = if + /*height < f32::EPSILON as f64 * 0.5*//*false*/ + /*true*/ + false { + ( + /*(alt_main./*to_le_bytes()[7]*/to_bits() & 1) as f64*/ + (alt_main / CONFIG.mountain_scale as f64 * 128.0).mul(0.1).powf(1.2) * /*(1.0 / CONFIG.mountain_scale as f64)*/(f32::EPSILON * 0.5) as f64, + (f32::EPSILON * 0.5) as f64, + ) + } else { + (0.0, 0.0) }; + // tan(6/360*2*pi)*32 ~ 3.4 + // 3.4/32*512 ~ 54 + // 18/32*512 ~ 288 + // tan(pi/6)*32 ~ 18 + // tan(54/360*2*pi)*32 + // let height = 1.0f64; + let turb_wposf = wposf + .div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)) + .div(turb_wposf_div); + let turb = Vec2::new( + gen_ctx.turb_x_nz.get(turb_wposf.into_array()), + gen_ctx.turb_y_nz.get(turb_wposf.into_array()), + ) * uplift_turb_scale + * TerrainChunkSize::RECT_SIZE.map(|e| e as f64); + let turb_wposf = wposf + turb; + let uheight = gen_ctx + .uplift_nz + .get(turb_wposf.into_array()) + /* .min(0.5) + .max(-0.5)*/ + .min(1.0) + .max(-1.0) + .mul(0.5) + .add(0.5); + // u = 1e-3: normal-high (dike, mountain) + // u = 5e-4: normal (mid example in Yuan, average mountain uplift) + // u = 2e-4: low (low example in Yuan; known that lagoons etc. may have u ~ 0.05). + // u = 0: low (plateau [fan, altitude = 0.0]) + // let height = uheight; + // let height = 1.0f64; + + // let height = 1.0 / 7.0f64; + // let height = 0.0 / 31.0f64; + let bfrac = /*erosion_factor(0.5);*/0.0; + let height = (height - bfrac).abs().div(1.0 - bfrac); + let height = height + /* .mul(31.0 / 32.0) + .add(1.0 / 32.0) */ + /* .mul(15.0 / 16.0) + .add(1.0 / 16.0) */ + /* .mul(5.0 / 8.0) + .add(3.0 / 8.0) */ + /* .mul(6.0 / 8.0) + .add(2.0 / 8.0) */ + /* .mul(7.0 / 8.0) + .add(1.0 / 8.0) */ + .mul(max_erosion_per_delta_t) + .sub(/*1.0 / CONFIG.mountain_scale as f64*/ bump_max) + .add(bump_factor); + /* .sub(/*1.0 / CONFIG.mountain_scale as f64*/(f32::EPSILON * 0.5) as f64) + .add(bump_factor); */ + height as f32 + }; let alt_func = |posi| { if is_ocean_fn(posi) { // -max_erosion_per_delta_t as f32 @@ -1062,9 +1084,8 @@ impl WorldSim { .mul(alt_main.powf(0.8).max(/*0.25*/ 0.15)) .mul(0.3) .add(1.0) - .mul(0.4) - /* + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0) - .mul(0.045)*/) + .mul(0.4)/* + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0) + .mul(0.045)*/) }; // (kf_func(posi) / 1.5e-4 * CONFIG.mountain_scale as f64) as f32 @@ -1112,7 +1133,7 @@ impl WorldSim { }; let reader = BufReader::new(file); - let map : WorldFile = match bincode::deserialize_from(reader) { + let map: WorldFile = match bincode::deserialize_from(reader) { Ok(map) => map, Err(err) => { log::warn!("Couldn't parse map: {:?})", err); @@ -1120,7 +1141,9 @@ impl WorldSim { } }; - if map.alt.len() != map.basement.len() || map.alt.len() != WORLD_SIZE.x as usize * WORLD_SIZE.y as usize { + if map.alt.len() != map.basement.len() + || map.alt.len() != WORLD_SIZE.x as usize * WORLD_SIZE.y as usize + { log::warn!("World size of map is invalid."); return None; } @@ -1149,8 +1172,16 @@ impl WorldSim { n_steps, &river_seed, &rock_strength_nz, - |posi| alt_func(posi),// + if is_ocean_fn(posi) { 0.0 } else { 128.0 }, - |posi| alt_func(posi) - if is_ocean_fn(posi) { 0.0 } else { /*1400.0*//*CONFIG.mountain_scale * 0.75*/0.0 },// if is_ocean_fn(posi) { old_height(posi) } else { 0.0 }, + |posi| alt_func(posi), // + if is_ocean_fn(posi) { 0.0 } else { 128.0 }, + |posi| { + alt_func(posi) + - if is_ocean_fn(posi) { + 0.0 + } else { + /*1400.0*//*CONFIG.mountain_scale * 0.75*/ + 0.0 + } + }, // if is_ocean_fn(posi) { old_height(posi) } else { 0.0 }, is_ocean_fn, uplift_fn, |posi| n_func(posi), @@ -1180,13 +1211,10 @@ impl WorldSim { }; // Save map, if necessary. - let map = WorldFile { - alt, - basement, - }; + let map = WorldFile { alt, basement }; (|| { if let FileOpts::Save = opts.world_file { - use std::{time::SystemTime}; + use std::time::SystemTime; // Check if folder exists and create it if it does not let mut path = PathBuf::from("./maps"); if !path.exists() { @@ -1243,7 +1271,10 @@ impl WorldSim { let is_ocean = get_oceans(|posi| alt[posi]); let is_ocean_fn = |posi: usize| is_ocean[posi]; - let mut dh = downhill(|posi| alt[posi] as f32/*&alt*/, /*old_height*/ is_ocean_fn); + let mut dh = downhill( + |posi| alt[posi] as f32, /*&alt*/ + /*old_height*/ is_ocean_fn, + ); let (boundary_len, indirection, water_alt_pos, _) = get_lakes(/*&/*water_alt*/alt*/|posi| alt[posi] as f32, &mut dh); let flux_old = get_drainage(&water_alt_pos, &dh, boundary_len); @@ -1258,7 +1289,9 @@ impl WorldSim { /* // Find the pass this lake is flowing into (i.e. water at the lake bottom gets // pushed towards the point identified by pass_idx). let neighbor_pass_idx = dh[lake_idx]; */ - let chunk_water_alt = if /*neighbor_pass_idx*/dh[lake_idx] < 0 { + let chunk_water_alt = if + /*neighbor_pass_idx*/ + dh[lake_idx] < 0 { // This is either a boundary node (dh[chunk_idx] == -2, i.e. water is at sea level) // or part of a lake that flows directly into the ocean. In the former case, water // is at sea level so we just return 0.0. In the latter case, the lake bottom must @@ -1293,9 +1326,9 @@ impl WorldSim { let water_alt = fill_sinks(water_height_initial, is_ocean_fn); /* let water_alt = (0..WORLD_SIZE.x * WORLD_SIZE.y) - .into_par_iter() - .map(|posi| water_height_initial(posi)) - .collect::>(); */ + .into_par_iter() + .map(|posi| water_height_initial(posi)) + .collect::>(); */ let rivers = get_rivers(&water_alt_pos, &water_alt, &dh, &indirection, &flux_old); @@ -1312,7 +1345,9 @@ impl WorldSim { /* // Find the pass this lake is flowing into (i.e. water at the lake bottom gets // pushed towards the point identified by pass_idx). let neighbor_pass_idx = dh[lake_idx]; */ - if /*neighbor_pass_idx*/dh[lake_idx] < 0 { + if + /*neighbor_pass_idx*/ + dh[lake_idx] < 0 { // This is either a boundary node (dh[chunk_idx] == -2, i.e. water is at sea level) // or part of a lake that flows directly into the ocean. In the former case, water // is at sea level so we just return 0.0. In the latter case, the lake bottom must @@ -1970,9 +2005,7 @@ impl SimChunk { let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale; let mut alt = CONFIG.sea_level.add(alt_pre.div(height_scale)); let mut basement = CONFIG.sea_level.add(basement_pre.div(height_scale)); - let water_alt = CONFIG - .sea_level - .add(water_alt_pre.div(height_scale)); + let water_alt = CONFIG.sea_level.add(water_alt_pre.div(height_scale)); let downhill = if downhill_pre == -2 { None } else if downhill_pre < 0 { @@ -2030,7 +2063,9 @@ impl SimChunk { } // No trees in the ocean, with zero humidity (currently), or directly on bedrock. - let tree_density = if is_underwater/* || alt - basement.min(alt) < 2.0 */ { + let tree_density = if is_underwater + /* || alt - basement.min(alt) < 2.0 */ + { 0.0 } else { let tree_density = (gen_ctx.tree_nz.get((wposf.div(1024.0)).into_array())) diff --git a/world/src/sim/util.rs b/world/src/sim/util.rs index fd587de445..65bf6886de 100644 --- a/world/src/sim/util.rs +++ b/world/src/sim/util.rs @@ -218,7 +218,7 @@ pub fn local_cells(posi: usize) -> impl Clone + Iterator { } // NOTE: want to keep this such that the chunk index is in ascending order! -pub const NEIGHBOR_DELTA : [(i32, i32); 8] = [ +pub const NEIGHBOR_DELTA: [(i32, i32); 8] = [ (-1, -1), (0, -1), (1, -1), @@ -233,12 +233,12 @@ pub const NEIGHBOR_DELTA : [(i32, i32); 8] = [ pub fn neighbors(posi: usize) -> impl Clone + Iterator { let pos = uniform_idx_as_vec2(posi); NEIGHBOR_DELTA - .iter() - .map(move |&(x, y)| Vec2::new(pos.x + x, pos.y + y)) - .filter(|pos| { - pos.x >= 0 && pos.y >= 0 && pos.x < WORLD_SIZE.x as i32 && pos.y < WORLD_SIZE.y as i32 - }) - .map(vec2_as_uniform_idx) + .iter() + .map(move |&(x, y)| Vec2::new(pos.x + x, pos.y + y)) + .filter(|pos| { + pos.x >= 0 && pos.y >= 0 && pos.x < WORLD_SIZE.x as i32 && pos.y < WORLD_SIZE.y as i32 + }) + .map(vec2_as_uniform_idx) } // Note that we should already have okay cache locality since we have a grid. @@ -249,10 +249,14 @@ pub fn uphill<'a>(dh: &'a [isize], posi: usize) -> impl Clone + Iterator(h: impl Fn(usize) -> F + Sync, is_ocean: impl Fn(usize) -> bool + Sync) -> Box<[isize]> { +pub fn downhill( + h: impl Fn(usize) -> F + Sync, + is_ocean: impl Fn(usize) -> bool + Sync, +) -> Box<[isize]> { // Constructs not only the list of downhill nodes, but also computes an ordering (visiting // nodes in order from roots to leaves). - (0..WORLD_SIZE.x * WORLD_SIZE.y).into_par_iter() + (0..WORLD_SIZE.x * WORLD_SIZE.y) + .into_par_iter() // .enumerate() .map(|(posi/*, &nh*/)| { let nh = h(posi); @@ -702,5 +706,3 @@ impl<'a, F: NoiseFn + 'a, T> NoiseFn for ScaleBias<'a, F> { (self.source.get(point) * self.scale) + self.bias } } - -