veloren_world/sim/
mod.rs

1mod diffusion;
2mod erosion;
3mod location;
4mod map;
5mod util;
6mod way;
7
8// Reexports
9use self::erosion::Compute;
10pub use self::{
11    diffusion::diffusion,
12    location::Location,
13    map::{sample_pos, sample_wpos},
14    util::get_horizon_map,
15    way::{Path, Way},
16};
17pub(crate) use self::{
18    erosion::{
19        Alt, RiverData, RiverKind, do_erosion, fill_sinks, get_lakes, get_multi_drainage,
20        get_multi_rec, get_rivers,
21    },
22    util::{
23        InverseCdf, cdf_irwin_hall, downhill, get_oceans, local_cells, map_edge_factor,
24        uniform_noise, uphill,
25    },
26};
27
28use crate::{
29    CONFIG, IndexRef,
30    all::{Environment, ForestKind, TreeAttr},
31    block::BlockGen,
32    civ::{Place, PointOfInterest},
33    column::ColumnGen,
34    site::Site,
35    util::{
36        CARDINALS, DHashSet, FastNoise, FastNoise2d, LOCALITY, NEIGHBORS, RandomField, Sampler,
37        StructureGen2d, seed_expan,
38    },
39};
40use common::{
41    assets::{self, AssetExt},
42    calendar::Calendar,
43    grid::Grid,
44    lottery::Lottery,
45    resources::MapKind,
46    spiral::Spiral2d,
47    spot::Spot,
48    store::{Id, Store},
49    terrain::{
50        BiomeKind, CoordinateConversions, MapSizeLg, TerrainChunk, TerrainChunkSize,
51        map::MapConfig, uniform_idx_as_vec2, vec2_as_uniform_idx,
52    },
53    vol::RectVolSize,
54};
55use common_base::prof_span;
56use common_net::msg::WorldMapMsg;
57use noise::{
58    BasicMulti, Billow, Fbm, HybridMulti, MultiFractal, NoiseFn, Perlin, RidgedMulti, SuperSimplex,
59    core::worley::distance_functions,
60};
61use num::{Float, Signed, traits::FloatConst};
62use rand::{Rng, SeedableRng};
63use rand_chacha::ChaChaRng;
64use rayon::prelude::*;
65use serde::{Deserialize, Serialize};
66use std::{
67    f32,
68    fs::File,
69    io::{BufReader, BufWriter},
70    ops::{Add, Div, Mul, Neg, Sub},
71    path::PathBuf,
72    sync::Arc,
73};
74use strum::IntoEnumIterator;
75use tracing::{debug, info, warn};
76use vek::*;
77
78/// Default base two logarithm of the world size, in chunks, per dimension.
79///
80/// Currently, our default map dimensions are 2^10 × 2^10 chunks,
81/// mostly for historical reasons.  It is likely that we will increase this
82/// default at some point.
83const DEFAULT_WORLD_CHUNKS_LG: MapSizeLg =
84    if let Ok(map_size_lg) = MapSizeLg::new(Vec2 { x: 10, y: 10 }) {
85        map_size_lg
86    } else {
87        panic!("Default world chunk size does not satisfy required invariants.");
88    };
89
90/// A structure that holds cached noise values and cumulative distribution
91/// functions for the input that led to those values.  See the definition of
92/// InverseCdf for a description of how to interpret the types of its fields.
93struct GenCdf {
94    humid_base: InverseCdf,
95    temp_base: InverseCdf,
96    chaos: InverseCdf,
97    alt: Box<[Alt]>,
98    basement: Box<[Alt]>,
99    water_alt: Box<[f32]>,
100    dh: Box<[isize]>,
101    /// NOTE: Until we hit 4096 × 4096, this should suffice since integers with
102    /// an absolute value under 2^24 can be exactly represented in an f32.
103    flux: Box<[Compute]>,
104    pure_flux: InverseCdf<Compute>,
105    alt_no_water: InverseCdf,
106    rivers: Box<[RiverData]>,
107}
108
109pub(crate) struct GenCtx {
110    pub turb_x_nz: SuperSimplex,
111    pub turb_y_nz: SuperSimplex,
112    pub chaos_nz: RidgedMulti<Perlin>,
113    pub alt_nz: util::HybridMulti<Perlin>,
114    pub hill_nz: SuperSimplex,
115    pub temp_nz: Fbm<Perlin>,
116    // Humidity noise
117    pub humid_nz: Billow<Perlin>,
118    // Small amounts of noise for simulating rough terrain.
119    pub small_nz: BasicMulti<Perlin>,
120    pub rock_nz: HybridMulti<Perlin>,
121    pub tree_nz: BasicMulti<Perlin>,
122
123    // TODO: unused, remove??? @zesterer
124    pub _cave_0_nz: SuperSimplex,
125    pub _cave_1_nz: SuperSimplex,
126
127    pub structure_gen: StructureGen2d,
128    pub _big_structure_gen: StructureGen2d,
129    pub _region_gen: StructureGen2d,
130
131    pub _fast_turb_x_nz: FastNoise,
132    pub _fast_turb_y_nz: FastNoise,
133
134    pub _town_gen: StructureGen2d,
135    pub river_seed: RandomField,
136    pub rock_strength_nz: Fbm<Perlin>,
137    pub uplift_nz: util::Worley,
138}
139
140#[derive(Clone, Debug, Deserialize, Serialize)]
141#[serde(default)]
142pub struct GenOpts {
143    pub x_lg: u32,
144    pub y_lg: u32,
145    pub scale: f64,
146    pub map_kind: MapKind,
147    pub erosion_quality: f32,
148}
149
150impl Default for GenOpts {
151    fn default() -> Self {
152        Self {
153            x_lg: 10,
154            y_lg: 10,
155            scale: 2.0,
156            map_kind: MapKind::Square,
157            erosion_quality: 1.0,
158        }
159    }
160}
161
162#[derive(Clone, Debug, Deserialize, Serialize)]
163pub enum FileOpts {
164    /// If set, generate the world map and do not try to save to or load from
165    /// file (default).
166    Generate(GenOpts),
167    /// If set, generate the world map and save the world file (path is created
168    /// the same way screenshot paths are).
169    Save(PathBuf, GenOpts),
170    /// Combination of Save and Load.
171    /// Load map if exists or generate the world map and save the
172    /// world file.
173    LoadOrGenerate {
174        name: String,
175        #[serde(default)]
176        opts: GenOpts,
177        #[serde(default)]
178        overwrite: bool,
179    },
180    /// If set, load the world file from this path in legacy format (errors if
181    /// path not found).  This option may be removed at some point, since it
182    /// only applies to maps generated before map saving was merged into
183    /// master.
184    LoadLegacy(PathBuf),
185    /// If set, load the world file from this path (errors if path not found).
186    Load(PathBuf),
187    /// If set, look for  the world file at this asset specifier (errors if
188    /// asset is not found).
189    ///
190    /// NOTE: Could stand to merge this with `Load` and construct an enum that
191    /// can handle either a PathBuf or an asset specifier, at some point.
192    LoadAsset(String),
193}
194
195impl Default for FileOpts {
196    fn default() -> Self { Self::Generate(GenOpts::default()) }
197}
198
199impl FileOpts {
200    fn load_content(&self) -> (Option<ModernMap>, MapSizeLg, GenOpts) {
201        let parsed_world_file = self.try_load_map();
202
203        let mut gen_opts = self.gen_opts().unwrap_or_default();
204
205        let map_size_lg = if let Some(map) = &parsed_world_file {
206            MapSizeLg::new(map.map_size_lg)
207                .expect("World size of loaded map does not satisfy invariants.")
208        } else {
209            self.map_size()
210        };
211
212        // NOTE: Change 1.0 to 4.0 for a 4x
213        // improvement in world detail.  We also use this to automatically adjust
214        // grid_scale (multiplying by 4.0) and multiply mins_per_sec by
215        // 1.0 / (4.0 * 4.0) in ./erosion.rs, in order to get a similar rate of river
216        // formation.
217        //
218        // FIXME: This is a hack!  At some point we will have a more principled way of
219        // dealing with this.
220        if let Some(map) = &parsed_world_file {
221            gen_opts.scale = map.continent_scale_hack;
222        };
223
224        (parsed_world_file, map_size_lg, gen_opts)
225    }
226
227    fn gen_opts(&self) -> Option<GenOpts> {
228        match self {
229            Self::Generate(opts) | Self::Save(_, opts) | Self::LoadOrGenerate { opts, .. } => {
230                Some(opts.clone())
231            },
232            _ => None,
233        }
234    }
235
236    // TODO: this should return Option so that caller can choose fallback
237    fn map_size(&self) -> MapSizeLg {
238        match self {
239            Self::Generate(opts) | Self::Save(_, opts) | Self::LoadOrGenerate { opts, .. } => {
240                MapSizeLg::new(Vec2 {
241                    x: opts.x_lg,
242                    y: opts.y_lg,
243                })
244                .unwrap_or_else(|e| {
245                    warn!("World size does not satisfy invariants: {:?}", e);
246                    DEFAULT_WORLD_CHUNKS_LG
247                })
248            },
249            _ => DEFAULT_WORLD_CHUNKS_LG,
250        }
251    }
252
253    // TODO: This should probably return a Result, so that caller can choose
254    // whether to log error
255    fn try_load_map(&self) -> Option<ModernMap> {
256        let map = match self {
257            Self::LoadLegacy(path) => {
258                let file = match File::open(path) {
259                    Ok(file) => file,
260                    Err(e) => {
261                        warn!(?e, ?path, "Couldn't read path for maps");
262                        return None;
263                    },
264                };
265
266                let reader = BufReader::new(file);
267                let map: WorldFileLegacy = match bincode::deserialize_from(reader) {
268                    Ok(map) => map,
269                    Err(e) => {
270                        warn!(
271                            ?e,
272                            "Couldn't parse legacy map.  Maybe you meant to try a regular load?"
273                        );
274                        return None;
275                    },
276                };
277
278                map.into_modern()
279            },
280            Self::Load(path) => {
281                let file = match File::open(path) {
282                    Ok(file) => file,
283                    Err(e) => {
284                        warn!(?e, ?path, "Couldn't read path for maps");
285                        return None;
286                    },
287                };
288
289                let reader = BufReader::new(file);
290                let map: WorldFile = match bincode::deserialize_from(reader) {
291                    Ok(map) => map,
292                    Err(e) => {
293                        warn!(
294                            ?e,
295                            "Couldn't parse modern map.  Maybe you meant to try a legacy load?"
296                        );
297                        return None;
298                    },
299                };
300
301                map.into_modern()
302            },
303            Self::LoadAsset(specifier) => match WorldFile::load_owned(specifier) {
304                Ok(map) => map.into_modern(),
305                Err(err) => {
306                    match err.reason().downcast_ref::<std::io::Error>() {
307                        Some(e) => {
308                            warn!(?e, ?specifier, "Couldn't read asset specifier for maps");
309                        },
310                        None => {
311                            warn!(
312                                ?err,
313                                "Couldn't parse modern map.  Maybe you meant to try a legacy load?"
314                            );
315                        },
316                    }
317                    return None;
318                },
319            },
320            Self::LoadOrGenerate {
321                opts, overwrite, ..
322            } => {
323                // `unwrap` is safe here, because LoadOrGenerate has its path
324                // always defined
325                let path = self.map_path().unwrap();
326
327                let file = match File::open(&path) {
328                    Ok(file) => file,
329                    Err(e) => {
330                        warn!(?e, ?path, "Couldn't find needed map. Generating...");
331                        return None;
332                    },
333                };
334
335                let reader = BufReader::new(file);
336                let map: WorldFile = match bincode::deserialize_from(reader) {
337                    Ok(map) => map,
338                    Err(e) => {
339                        warn!(
340                            ?e,
341                            "Couldn't parse modern map.  Maybe you meant to try a legacy load?"
342                        );
343                        return None;
344                    },
345                };
346
347                // FIXME:
348                // We check if we need to generate new map by comparing gen opts.
349                // But we also have another generation paramater that currently
350                // passed outside and used for both worldsim and worldgen.
351                //
352                // Ideally, we need to figure out how we want to use seed, i. e.
353                // moving worldgen seed to gen opts and use different sim seed from
354                // server config or grab sim seed from world file.
355                //
356                // NOTE: we intentionally use pattern-matching here to get
357                // options, so that when gen opts get another field, compiler
358                // will force you to update following logic
359                let GenOpts {
360                    x_lg, y_lg, scale, ..
361                } = opts;
362                let map = match map {
363                    WorldFile::Veloren0_7_0(map) => map,
364                    WorldFile::Veloren0_5_0(_) => {
365                        panic!("World file v0.5.0 isn't supported with LoadOrGenerate.")
366                    },
367                };
368
369                if map.continent_scale_hack != *scale || map.map_size_lg != Vec2::new(*x_lg, *y_lg)
370                {
371                    if *overwrite {
372                        warn!(
373                            "{}\n{}",
374                            "Specified options don't correspond to these in loaded map.",
375                            "Map will be regenerated and overwritten."
376                        );
377                    } else {
378                        panic!(
379                            "{}\n{}",
380                            "Specified options don't correspond to these in loaded map.",
381                            "Use 'ovewrite' option, if you wish to regenerate map."
382                        );
383                    }
384
385                    return None;
386                }
387
388                map.into_modern()
389            },
390            Self::Generate { .. } | Self::Save { .. } => return None,
391        };
392
393        match map {
394            Ok(map) => Some(map),
395            Err(e) => {
396                match e {
397                    WorldFileError::WorldSizeInvalid => {
398                        warn!("World size of map is invalid.");
399                    },
400                }
401                None
402            },
403        }
404    }
405
406    fn map_path(&self) -> Option<PathBuf> {
407        // TODO: Work out a nice bincode file extension.
408        match self {
409            Self::Save(path, _) => Some(PathBuf::from(&path)),
410            Self::LoadOrGenerate { name, .. } => {
411                const MAP_DIR: &str = "./maps";
412                let file_name = format!("{}.bin", name);
413                Some(std::path::Path::new(MAP_DIR).join(file_name))
414            },
415            _ => None,
416        }
417    }
418
419    fn save(&self, map: &WorldFile) {
420        let path = if let Some(path) = self.map_path() {
421            path
422        } else {
423            return;
424        };
425
426        // Check if folder exists and create it if it does not
427        let map_dir = path.parent().expect("failed to get map directory");
428        if !map_dir.exists() {
429            if let Err(e) = std::fs::create_dir_all(map_dir) {
430                warn!(?e, ?map_dir, "Couldn't create folder for map");
431                return;
432            }
433        }
434
435        let file = match File::create(path.clone()) {
436            Ok(file) => file,
437            Err(e) => {
438                warn!(?e, ?path, "Couldn't create file for maps");
439                return;
440            },
441        };
442
443        let writer = BufWriter::new(file);
444        if let Err(e) = bincode::serialize_into(writer, map) {
445            warn!(?e, "Couldn't write map");
446        }
447        if let Ok(p) = std::fs::canonicalize(path) {
448            info!("Map saved at {}", p.to_string_lossy());
449        }
450    }
451}
452
453pub struct WorldOpts {
454    /// Set to false to disable seeding elements during worldgen.
455    pub seed_elements: bool,
456    pub world_file: FileOpts,
457    pub calendar: Option<Calendar>,
458}
459
460impl Default for WorldOpts {
461    fn default() -> Self {
462        Self {
463            seed_elements: true,
464            world_file: Default::default(),
465            calendar: None,
466        }
467    }
468}
469
470/// LEGACY: Remove when people stop caring.
471#[derive(Serialize, Deserialize)]
472#[repr(C)]
473pub struct WorldFileLegacy {
474    /// Saved altitude height map.
475    pub alt: Box<[Alt]>,
476    /// Saved basement height map.
477    pub basement: Box<[Alt]>,
478}
479
480/// Version of the world map intended for use in Veloren 0.5.0.
481#[derive(Serialize, Deserialize)]
482#[repr(C)]
483pub struct WorldMap_0_5_0 {
484    /// Saved altitude height map.
485    pub alt: Box<[Alt]>,
486    /// Saved basement height map.
487    pub basement: Box<[Alt]>,
488}
489
490/// Version of the world map intended for use in Veloren 0.7.0.
491#[derive(Serialize, Deserialize)]
492#[repr(C)]
493pub struct WorldMap_0_7_0 {
494    /// Saved map size.
495    pub map_size_lg: Vec2<u32>,
496    /// Saved continent_scale hack, to try to better approximate the correct
497    /// seed according to varying map size.
498    ///
499    /// TODO: Remove when generating new maps becomes more principled.
500    pub continent_scale_hack: f64,
501    /// Saved altitude height map.
502    pub alt: Box<[Alt]>,
503    /// Saved basement height map.
504    pub basement: Box<[Alt]>,
505}
506
507/// Errors when converting a map to the most recent type (currently,
508/// shared by the various map types, but at some point we might switch to
509/// version-specific errors if it feels worthwhile).
510#[derive(Debug)]
511pub enum WorldFileError {
512    /// Map size was invalid, and it can't be converted to a valid one.
513    WorldSizeInvalid,
514}
515
516/// WORLD MAP.
517///
518/// A way to store certain components between runs of map generation.  Only
519/// intended for development purposes--no attempt is made to detect map
520/// invalidation or make sure that the map is synchronized with updates to
521/// noise-rs, changes to other parameters, etc.
522///
523/// The map is versioned to enable format detection between versions of Veloren,
524/// so that when we update the map format we don't break existing maps (or at
525/// least, we will try hard not to break maps between versions; if we can't
526/// avoid it, we can at least give a reasonable error message).
527///
528/// NOTE: We rely somewhat heavily on the implementation specifics of bincode
529/// to make sure this is backwards compatible.  When adding new variants here,
530/// Be very careful to make sure tha the old variants are preserved in the
531/// correct order and with the correct names and indices, and make sure to keep
532/// the #[repr(u32)]!
533///
534/// All non-legacy versions of world files should (ideally) fit in this format.
535/// Since the format contains a version and is designed to be extensible
536/// backwards-compatibly, the only reason not to use this forever would be if we
537/// decided to move away from BinCode, or store data across multiple files (or
538/// something else weird I guess).
539///
540/// Update this when you add a new map version.
541#[derive(Serialize, Deserialize)]
542#[repr(u32)]
543pub enum WorldFile {
544    Veloren0_5_0(WorldMap_0_5_0) = 0,
545    Veloren0_7_0(WorldMap_0_7_0) = 1,
546}
547
548impl assets::Asset for WorldFile {
549    type Loader = assets::BincodeLoader;
550
551    const EXTENSION: &'static str = "bin";
552}
553
554/// Data for the most recent map type.  Update this when you add a new map
555/// version.
556pub type ModernMap = WorldMap_0_7_0;
557
558/// The default world map.
559///
560/// TODO: Consider using some naming convention to automatically change this
561/// with changing versions, or at least keep it in a constant somewhere that's
562/// easy to change.
563// Generation parameters:
564//
565// gen_opts: (
566//     erosion_quality: 1.0,
567//     map_kind: Circle,
568//     scale: 2.157574498096227,
569//     x_lg: 10,
570//     y_lg: 10,
571// )
572// seed: 3582734543
573//
574// The biome seed can found below
575pub const DEFAULT_WORLD_MAP: &str = "world.map.veloren_0_18_0_0";
576/// This is *not* the seed used to generate the default map, this seed was used
577/// to generate a better set of biomes on it as the original ones were
578/// unsuitable.
579///
580/// See DEFAULT_WORLD_MAP to get the original worldgen parameters.
581pub const DEFAULT_WORLD_SEED: u32 = 130626853;
582
583impl WorldFileLegacy {
584    #[inline]
585    /// Idea: each map type except the latest knows how to transform
586    /// into the the subsequent map version, and each map type including the
587    /// latest exposes an "into_modern()" method that converts this map type
588    /// to the modern map type.  Thus, to migrate a map from an old format to a
589    /// new format, we just need to transform the old format to the
590    /// subsequent map version, and then call .into_modern() on that--this
591    /// should construct a call chain that ultimately ends up with a modern
592    /// version.
593    pub fn into_modern(self) -> Result<ModernMap, WorldFileError> {
594        // NOTE: At this point, we assume that any remaining legacy maps were 1024 ×
595        // 1024.
596        if self.alt.len() != self.basement.len() || self.alt.len() != 1024 * 1024 {
597            return Err(WorldFileError::WorldSizeInvalid);
598        }
599
600        let map = WorldMap_0_5_0 {
601            alt: self.alt,
602            basement: self.basement,
603        };
604
605        map.into_modern()
606    }
607}
608
609impl WorldMap_0_5_0 {
610    #[inline]
611    pub fn into_modern(self) -> Result<ModernMap, WorldFileError> {
612        let pow_size = (self.alt.len().trailing_zeros()) / 2;
613        let two_coord_size = 1 << (2 * pow_size);
614        if self.alt.len() != self.basement.len() || self.alt.len() != two_coord_size {
615            return Err(WorldFileError::WorldSizeInvalid);
616        }
617
618        // The recommended continent scale for maps from version 0.5.0 is (in all
619        // existing cases) just 1.0 << (f64::from(pow_size) - 10.0).
620        let continent_scale_hack = (f64::from(pow_size) - 10.0).exp2();
621
622        let map = WorldMap_0_7_0 {
623            map_size_lg: Vec2::new(pow_size, pow_size),
624            continent_scale_hack,
625            alt: self.alt,
626            basement: self.basement,
627        };
628
629        map.into_modern()
630    }
631}
632
633impl WorldMap_0_7_0 {
634    #[inline]
635    pub fn into_modern(self) -> Result<ModernMap, WorldFileError> {
636        if self.alt.len() != self.basement.len()
637            || self.alt.len() != (1 << (self.map_size_lg.x + self.map_size_lg.y))
638            || self.continent_scale_hack <= 0.0
639        {
640            return Err(WorldFileError::WorldSizeInvalid);
641        }
642
643        Ok(self)
644    }
645}
646
647impl WorldFile {
648    /// Turns map data from the latest version into a versioned WorldFile ready
649    /// for serialization. Whenever a new map is updated, just change the
650    /// variant we construct here to make sure we're using the latest map
651    /// version.
652    pub fn new(map: ModernMap) -> Self { WorldFile::Veloren0_7_0(map) }
653
654    #[inline]
655    /// Turns a WorldFile into the latest version.  Whenever a new map version
656    /// is added, just add it to this match statement.
657    pub fn into_modern(self) -> Result<ModernMap, WorldFileError> {
658        match self {
659            WorldFile::Veloren0_5_0(map) => map.into_modern(),
660            WorldFile::Veloren0_7_0(map) => map.into_modern(),
661        }
662    }
663}
664
665#[derive(Debug)]
666pub enum WorldSimStage {
667    // TODO: Add more stages
668    Erosion {
669        progress: f64,
670        estimate: Option<std::time::Duration>,
671    },
672}
673
674pub struct WorldSim {
675    pub seed: u32,
676    /// Base 2 logarithm of the map size.
677    map_size_lg: MapSizeLg,
678    /// Maximum height above sea level of any chunk in the map (not including
679    /// post-erosion warping, cliffs, and other things like that).
680    pub max_height: f32,
681    pub(crate) chunks: Vec<SimChunk>,
682    //TODO: remove or use this property
683    pub(crate) _locations: Vec<Location>,
684
685    pub(crate) gen_ctx: GenCtx,
686    pub rng: ChaChaRng,
687
688    pub(crate) calendar: Option<Calendar>,
689}
690
691impl WorldSim {
692    pub fn empty() -> Self {
693        let gen_ctx = GenCtx {
694            turb_x_nz: SuperSimplex::new(0),
695            turb_y_nz: SuperSimplex::new(0),
696            chaos_nz: RidgedMulti::new(0),
697            hill_nz: SuperSimplex::new(0),
698            alt_nz: util::HybridMulti::new(0),
699            temp_nz: Fbm::new(0),
700
701            small_nz: BasicMulti::new(0),
702            rock_nz: HybridMulti::new(0),
703            tree_nz: BasicMulti::new(0),
704            _cave_0_nz: SuperSimplex::new(0),
705            _cave_1_nz: SuperSimplex::new(0),
706
707            structure_gen: StructureGen2d::new(0, 24, 10),
708            _big_structure_gen: StructureGen2d::new(0, 768, 512),
709            _region_gen: StructureGen2d::new(0, 400, 96),
710            humid_nz: Billow::new(0),
711
712            _fast_turb_x_nz: FastNoise::new(0),
713            _fast_turb_y_nz: FastNoise::new(0),
714
715            _town_gen: StructureGen2d::new(0, 2048, 1024),
716            river_seed: RandomField::new(0),
717            rock_strength_nz: Fbm::new(0),
718            uplift_nz: util::Worley::new(0),
719        };
720        Self {
721            seed: 0,
722            map_size_lg: MapSizeLg::new(Vec2::one()).unwrap(),
723            max_height: 0.0,
724            chunks: vec![SimChunk {
725                chaos: 0.0,
726                alt: 0.0,
727                basement: 0.0,
728                water_alt: 0.0,
729                downhill: None,
730                flux: 0.0,
731                temp: 0.0,
732                humidity: 0.0,
733                rockiness: 0.0,
734                tree_density: 0.0,
735                forest_kind: ForestKind::Dead,
736                spawn_rate: 0.0,
737                river: RiverData::default(),
738                surface_veg: 0.0,
739                sites: vec![],
740                place: None,
741                poi: None,
742                path: Default::default(),
743                cliff_height: 0.0,
744                spot: None,
745                contains_waypoint: false,
746            }],
747            _locations: Vec::new(),
748            gen_ctx,
749            rng: rand_chacha::ChaCha20Rng::from_seed([0; 32]),
750            calendar: None,
751        }
752    }
753
754    pub fn generate(
755        seed: u32,
756        opts: WorldOpts,
757        threadpool: &rayon::ThreadPool,
758        stage_report: &dyn Fn(WorldSimStage),
759    ) -> Self {
760        prof_span!("WorldSim::generate");
761        let calendar = opts.calendar; // separate lifetime of elements
762        let world_file = opts.world_file;
763
764        // Parse out the contents of various map formats into the values we need.
765        let (parsed_world_file, map_size_lg, gen_opts) = world_file.load_content();
766        // Currently only used with LoadOrGenerate to know if we need to
767        // overwrite world file
768        let fresh = parsed_world_file.is_none();
769
770        let mut rng = ChaChaRng::from_seed(seed_expan::rng_state(seed));
771        let continent_scale = gen_opts.scale
772            * 5_000.0f64
773                .div(32.0)
774                .mul(TerrainChunkSize::RECT_SIZE.x as f64);
775        let rock_lacunarity = 2.0;
776        let uplift_scale = 128.0;
777        let uplift_turb_scale = uplift_scale / 4.0;
778
779        info!("Starting world generation");
780
781        // NOTE: Changing order will significantly change WorldGen, so try not to!
782        let gen_ctx = GenCtx {
783            turb_x_nz: SuperSimplex::new(rng.gen()),
784            turb_y_nz: SuperSimplex::new(rng.gen()),
785            chaos_nz: RidgedMulti::new(rng.gen()).set_octaves(7).set_frequency(
786                RidgedMulti::<Perlin>::DEFAULT_FREQUENCY * (5_000.0 / continent_scale),
787            ),
788            hill_nz: SuperSimplex::new(rng.gen()),
789            alt_nz: util::HybridMulti::new(rng.gen())
790                .set_octaves(8)
791                .set_frequency(10_000.0 / continent_scale)
792                // persistence = lacunarity^(-(1.0 - fractal increment))
793                .set_lacunarity(util::HybridMulti::<Perlin>::DEFAULT_LACUNARITY)
794                .set_persistence(util::HybridMulti::<Perlin>::DEFAULT_LACUNARITY.powi(-1))
795                .set_offset(0.0),
796            temp_nz: Fbm::new(rng.gen())
797                .set_octaves(6)
798                .set_persistence(0.5)
799                .set_frequency(1.0 / (((1 << 6) * 64) as f64))
800                .set_lacunarity(2.0),
801
802            small_nz: BasicMulti::new(rng.gen()).set_octaves(2),
803            rock_nz: HybridMulti::new(rng.gen()).set_persistence(0.3),
804            tree_nz: BasicMulti::new(rng.gen())
805                .set_octaves(12)
806                .set_persistence(0.75),
807            _cave_0_nz: SuperSimplex::new(rng.gen()),
808            _cave_1_nz: SuperSimplex::new(rng.gen()),
809
810            structure_gen: StructureGen2d::new(rng.gen(), 24, 10),
811            _big_structure_gen: StructureGen2d::new(rng.gen(), 768, 512),
812            _region_gen: StructureGen2d::new(rng.gen(), 400, 96),
813            humid_nz: Billow::new(rng.gen())
814                .set_octaves(9)
815                .set_persistence(0.4)
816                .set_frequency(0.2),
817
818            _fast_turb_x_nz: FastNoise::new(rng.gen()),
819            _fast_turb_y_nz: FastNoise::new(rng.gen()),
820
821            _town_gen: StructureGen2d::new(rng.gen(), 2048, 1024),
822            river_seed: RandomField::new(rng.gen()),
823            rock_strength_nz: Fbm::new(rng.gen())
824                .set_octaves(10)
825                .set_lacunarity(rock_lacunarity)
826                // persistence = lacunarity^(-(1.0 - fractal increment))
827                // NOTE: In paper, fractal increment is roughly 0.25.
828                .set_persistence(rock_lacunarity.powf(-0.75))
829                .set_frequency(
830                    1.0 * (5_000.0 / continent_scale)
831                        / (2.0 * TerrainChunkSize::RECT_SIZE.x as f64 * 2.0.powi(10 - 1)),
832                ),
833            uplift_nz: util::Worley::new(rng.gen())
834                .set_frequency(1.0 / (TerrainChunkSize::RECT_SIZE.x as f64 * uplift_scale))
835                .set_distance_function(distance_functions::euclidean),
836        };
837
838        let river_seed = &gen_ctx.river_seed;
839        let rock_strength_nz = &gen_ctx.rock_strength_nz;
840
841        // Suppose the old world has grid spacing Δx' = Δy', new Δx = Δy.
842        // We define grid_scale such that Δx = height_scale * Δx' ⇒
843        //  grid_scale = Δx / Δx'.
844        let grid_scale = 1.0f64 / (4.0 / gen_opts.scale)/*1.0*/;
845
846        // Now, suppose we want to generate a world with "similar" topography, defined
847        // in this case as having roughly equal slopes at steady state, with the
848        // simulation taking roughly as many steps to get to the point the
849        // previous world was at when it finished being simulated.
850        //
851        // Some computations with our coupled SPL/debris flow give us (for slope S
852        // constant) the following suggested scaling parameters to make this
853        // work:   k_fs_scale ≡ (K𝑓 / K𝑓') = grid_scale^(-2m) =
854        // grid_scale^(-2θn)
855        let k_fs_scale = |theta, n| grid_scale.powf(-2.0 * (theta * n) as f64);
856
857        //   k_da_scale ≡ (K_da / K_da') = grid_scale^(-2q)
858        let k_da_scale = |q| grid_scale.powf(-2.0 * q);
859        //
860        // Some other estimated parameters are harder to come by and *much* more
861        // dubious, not being accurate for the coupled equation. But for the SPL
862        // only one we roughly find, for h the height at steady state and time τ
863        // = time to steady state, with Hack's Law estimated b = 2.0 and various other
864        // simplifying assumptions, the estimate:
865        //   height_scale ≡ (h / h') = grid_scale^(n)
866        let height_scale = |n: f32| grid_scale.powf(n as f64) as Alt;
867        //   time_scale ≡ (τ / τ') = grid_scale^(n)
868        let time_scale = |n: f32| grid_scale.powf(n as f64);
869        //
870        // Based on this estimate, we have:
871        //   delta_t_scale ≡ (Δt / Δt') = time_scale
872        let delta_t_scale = time_scale;
873        //   alpha_scale ≡ (α / α') = height_scale^(-1)
874        let alpha_scale = |n: f32| height_scale(n).recip() as f32;
875        //
876        // Slightly more dubiously (need to work out the math better) we find:
877        //   k_d_scale ≡ (K_d / K_d') = grid_scale^2 / (/*height_scale * */ time_scale)
878        let k_d_scale = |n: f32| grid_scale.powi(2) / (/* height_scale(n) * */time_scale(n));
879        //   epsilon_0_scale ≡ (ε₀ / ε₀') = height_scale(n) / time_scale(n)
880        let epsilon_0_scale = |n| (height_scale(n) / time_scale(n) as Alt) as f32;
881
882        // Approximate n for purposes of computation of parameters above over the whole
883        // grid (when a chunk isn't available).
884        let n_approx = 1.0;
885        let max_erosion_per_delta_t = 64.0 * delta_t_scale(n_approx);
886        let n_steps = (100.0 * gen_opts.erosion_quality) as usize;
887        let n_small_steps = 0;
888        let n_post_load_steps = 0;
889
890        // Logistic regression.  Make sure x ∈ (0, 1).
891        let logit = |x: f64| x.ln() - (-x).ln_1p();
892        // 0.5 + 0.5 * tanh(ln(1 / (1 - 0.1) - 1) / (2 * (sqrt(3)/pi)))
893        let logistic_2_base = 3.0f64.sqrt() * std::f64::consts::FRAC_2_PI;
894        // Assumes μ = 0, σ = 1
895        let logistic_cdf = |x: f64| (x / logistic_2_base).tanh() * 0.5 + 0.5;
896
897        let map_size_chunks_len_f64 = map_size_lg.chunks().map(f64::from).product();
898        let min_epsilon = 1.0 / map_size_chunks_len_f64.max(f64::EPSILON * 0.5);
899        let max_epsilon = (1.0 - 1.0 / map_size_chunks_len_f64).min(1.0 - f64::EPSILON * 0.5);
900
901        // No NaNs in these uniform vectors, since the original noise value always
902        // returns Some.
903        let ((alt_base, _), (chaos, _)) = threadpool.join(
904            || {
905                uniform_noise(map_size_lg, |_, wposf| {
906                    match gen_opts.map_kind {
907                        MapKind::Square => {
908                            // "Base" of the chunk, to be multiplied by CONFIG.mountain_scale
909                            // (multiplied value is from -0.35 *
910                            // (CONFIG.mountain_scale * 1.05) to
911                            // 0.35 * (CONFIG.mountain_scale * 0.95), but value here is from -0.3675
912                            // to 0.3325).
913                            Some(
914                                (gen_ctx
915                                    .alt_nz
916                                    .get((wposf.div(10_000.0)).into_array())
917                                    .clamp(-1.0, 1.0))
918                                .sub(0.05)
919                                .mul(0.35),
920                            )
921                        },
922                        MapKind::Circle => {
923                            let world_sizef = map_size_lg.chunks().map(|e| e as f64)
924                                * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
925                            Some(
926                                (gen_ctx
927                                    .alt_nz
928                                    .get((wposf.div(5_000.0 * gen_opts.scale)).into_array())
929                                    .clamp(-1.0, 1.0))
930                                .add(
931                                    0.2 - ((wposf / world_sizef) * 2.0 - 1.0)
932                                        .magnitude_squared()
933                                        .powf(0.75)
934                                        .clamped(0.0, 1.0)
935                                        .powf(1.0)
936                                        * 0.6,
937                                )
938                                .mul(0.5),
939                            )
940                        },
941                    }
942                })
943            },
944            || {
945                uniform_noise(map_size_lg, |_, wposf| {
946                    // From 0 to 1.6, but the distribution before the max is from -1 and 1.6, so
947                    // there is a 50% chance that hill will end up at 0.3 or
948                    // lower, and probably a very high change it will be exactly
949                    // 0.
950                    let hill = (0.0f64
951                        + gen_ctx
952                            .hill_nz
953                            .get(
954                                (wposf
955                                    .mul(32.0)
956                                    .div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64))
957                                    .div(1_500.0))
958                                .into_array(),
959                            )
960                            .clamp(-1.0, 1.0)
961                            .mul(1.0)
962                        + gen_ctx
963                            .hill_nz
964                            .get(
965                                (wposf
966                                    .mul(32.0)
967                                    .div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64))
968                                    .div(400.0))
969                                .into_array(),
970                            )
971                            .clamp(-1.0, 1.0)
972                            .mul(0.3))
973                    .add(0.3)
974                    .max(0.0);
975
976                    // chaos produces a value in [0.12, 1.32].  It is a meta-level factor intended
977                    // to reflect how "chaotic" the region is--how much weird
978                    // stuff is going on on this terrain.
979                    Some(
980                        ((gen_ctx
981                            .chaos_nz
982                            .get((wposf.div(3_000.0)).into_array())
983                            .clamp(-1.0, 1.0))
984                        .add(1.0)
985                        .mul(0.5)
986                        // [0, 1] * [0.4, 1] = [0, 1] (but probably towards the lower end)
987                        .mul(
988                            (gen_ctx
989                                .chaos_nz
990                                .get((wposf.div(6_000.0)).into_array())
991                                .clamp(-1.0, 1.0))
992                            .abs()
993                                .clamp(0.4, 1.0),
994                        )
995                        // Chaos is always increased by a little when we're on a hill (but remember
996                        // that hill is 0.3 or less about 50% of the time).
997                        // [0, 1] + 0.2 * [0, 1.6] = [0, 1.32]
998                        .add(0.2 * hill)
999                        // We can't have *no* chaos!
1000                        .max(0.12)) as f32,
1001                    )
1002                })
1003            },
1004        );
1005
1006        // We ignore sea level because we actually want to be relative to sea level here
1007        // and want things in CONFIG.mountain_scale units, but otherwise this is
1008        // a correct altitude calculation.  Note that this is using the
1009        // "unadjusted" temperature.
1010        //
1011        // No NaNs in these uniform vectors, since the original noise value always
1012        // returns Some.
1013        let (alt_old, _) = uniform_noise(map_size_lg, |posi, wposf| {
1014            // This is the extension upwards from the base added to some extra noise from -1
1015            // to 1.
1016            //
1017            // The extra noise is multiplied by alt_main (the mountain part of the
1018            // extension) powered to 0.8 and clamped to [0.15, 1], to get a
1019            // value between [-1, 1] again.
1020            //
1021            // The sides then receive the sequence (y * 0.3 + 1.0) * 0.4, so we have
1022            // [-1*1*(1*0.3+1)*0.4, 1*(1*0.3+1)*0.4] = [-0.52, 0.52].
1023            //
1024            // Adding this to alt_main thus yields a value between -0.4 (if alt_main = 0 and
1025            // gen_ctx = -1, 0+-1*(0*.3+1)*0.4) and 1.52 (if alt_main = 1 and gen_ctx = 1).
1026            // Most of the points are above 0.
1027            //
1028            // Next, we add again by a sin of alt_main (between [-1, 1])^pow, getting
1029            // us (after adjusting for sign) another value between [-1, 1], and then this is
1030            // multiplied by 0.045 to get [-0.045, 0.045], which is added to [-0.4, 0.52] to
1031            // get [-0.445, 0.565].
1032            let alt_main = {
1033                // Extension upwards from the base.  A positive number from 0 to 1 curved to be
1034                // maximal at 0.  Also to be multiplied by CONFIG.mountain_scale.
1035                let alt_main = (gen_ctx
1036                    .alt_nz
1037                    .get((wposf.div(2_000.0)).into_array())
1038                    .clamp(-1.0, 1.0))
1039                .abs()
1040                .powf(1.35);
1041
1042                fn spring(x: f64, pow: f64) -> f64 { x.abs().powf(pow) * x.signum() }
1043
1044                0.0 + alt_main
1045                    + (gen_ctx
1046                        .small_nz
1047                        .get(
1048                            (wposf
1049                                .mul(32.0)
1050                                .div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64))
1051                                .div(300.0))
1052                            .into_array(),
1053                        )
1054                        .clamp(-1.0, 1.0))
1055                    .mul(alt_main.powf(0.8).max(/* 0.25 */ 0.15))
1056                    .mul(0.3)
1057                    .add(1.0)
1058                    .mul(0.4)
1059                    + spring(alt_main.abs().sqrt().min(0.75).mul(60.0).sin(), 4.0).mul(0.045)
1060            };
1061
1062            // Now we can compute the final altitude using chaos.
1063            // We multiply by chaos clamped to [0.1, 1.32] to get a value between [0.03,
1064            // 2.232] for alt_pre, then multiply by CONFIG.mountain_scale and
1065            // add to the base and sea level to get an adjusted value, then
1066            // multiply the whole thing by map_edge_factor (TODO: compute final
1067            // bounds).
1068            //
1069            // [-.3675, .3325] + [-0.445, 0.565] * [0.12, 1.32]^1.2
1070            // ~ [-.3675, .3325] + [-0.445, 0.565] * [0.07, 1.40]
1071            // = [-.3675, .3325] + ([-0.5785, 0.7345])
1072            // = [-0.946, 1.067]
1073            Some(
1074                ((alt_base[posi].1 + alt_main.mul((chaos[posi].1 as f64).powf(1.2)))
1075                    .mul(map_edge_factor(map_size_lg, posi) as f64)
1076                    .add(
1077                        (CONFIG.sea_level as f64)
1078                            .div(CONFIG.mountain_scale as f64)
1079                            .mul(map_edge_factor(map_size_lg, posi) as f64),
1080                    )
1081                    .sub((CONFIG.sea_level as f64).div(CONFIG.mountain_scale as f64)))
1082                    as f32,
1083            )
1084        });
1085
1086        // Calculate oceans.
1087        let is_ocean = get_oceans(map_size_lg, |posi: usize| alt_old[posi].1);
1088        // NOTE: Uncomment if you want oceans to exclusively be on the border of the
1089        // map.
1090        /* let is_ocean = (0..map_size_lg.chunks())
1091        .into_par_iter()
1092        .map(|i| map_edge_factor(map_size_lg, i) == 0.0)
1093        .collect::<Vec<_>>(); */
1094        let is_ocean_fn = |posi: usize| is_ocean[posi];
1095
1096        let turb_wposf_div = 8.0;
1097        let n_func = |posi| {
1098            if is_ocean_fn(posi) {
1099                return 1.0;
1100            }
1101            1.0
1102        };
1103        let old_height = |posi: usize| {
1104            alt_old[posi].1 * CONFIG.mountain_scale * height_scale(n_func(posi)) as f32
1105        };
1106
1107        // NOTE: Needed if you wish to use the distance to the point defining the Worley
1108        // cell, not just the value within that cell.
1109        // let uplift_nz_dist = gen_ctx.uplift_nz.clone().enable_range(true);
1110
1111        // Recalculate altitudes without oceans.
1112        // NaNs in these uniform vectors wherever is_ocean_fn returns true.
1113        let (alt_old_no_ocean, _) = uniform_noise(map_size_lg, |posi, _| {
1114            if is_ocean_fn(posi) {
1115                None
1116            } else {
1117                Some(old_height(posi))
1118            }
1119        });
1120        let (uplift_uniform, _) = uniform_noise(map_size_lg, |posi, _wposf| {
1121            if is_ocean_fn(posi) {
1122                None
1123            } else {
1124                let oheight = alt_old_no_ocean[posi].0 as f64 - 0.5;
1125                let height = (oheight + 0.5).powi(2);
1126                Some(height)
1127            }
1128        });
1129
1130        let alt_old_min_uniform = 0.0;
1131        let alt_old_max_uniform = 1.0;
1132
1133        let inv_func = |x: f64| x;
1134        let alt_exp_min_uniform = inv_func(min_epsilon);
1135        let alt_exp_max_uniform = inv_func(max_epsilon);
1136
1137        let erosion_factor = |x: f64| {
1138            (inv_func(x) - alt_exp_min_uniform) / (alt_exp_max_uniform - alt_exp_min_uniform)
1139        };
1140        let rock_strength_div_factor = (2.0 * TerrainChunkSize::RECT_SIZE.x as f64) / 8.0;
1141        let theta_func = |_posi| 0.4;
1142        let kf_func = {
1143            |posi| {
1144                let kf_scale_i = k_fs_scale(theta_func(posi), n_func(posi));
1145                if is_ocean_fn(posi) {
1146                    return 1.0e-4 * kf_scale_i;
1147                }
1148
1149                let kf_i = // kf = 1.5e-4: high-high (plateau [fan sediment])
1150                // kf = 1e-4: high (plateau)
1151                // kf = 2e-5: normal (dike [unexposed])
1152                // kf = 1e-6: normal-low (dike [exposed])
1153                // kf = 2e-6: low (mountain)
1154                // --
1155                // kf = 2.5e-7 to 8e-7: very low (Cordonnier papers on plate tectonics)
1156                // ((1.0 - uheight) * (1.5e-4 - 2.0e-6) + 2.0e-6) as f32
1157                //
1158                // ACTUAL recorded values worldwide: much lower...
1159                1.0e-6
1160                ;
1161                kf_i * kf_scale_i
1162            }
1163        };
1164        let kd_func = {
1165            |posi| {
1166                let n = n_func(posi);
1167                let kd_scale_i = k_d_scale(n);
1168                if is_ocean_fn(posi) {
1169                    let kd_i = 1.0e-2 / 4.0;
1170                    return kd_i * kd_scale_i;
1171                }
1172                // kd = 1e-1: high (mountain, dike)
1173                // kd = 1.5e-2: normal-high (plateau [fan sediment])
1174                // kd = 1e-2: normal (plateau)
1175                let kd_i = 1.0e-2 / 4.0;
1176                kd_i * kd_scale_i
1177            }
1178        };
1179        let g_func = |posi| {
1180            if map_edge_factor(map_size_lg, posi) == 0.0 {
1181                return 0.0;
1182            }
1183            // G = d* v_s / p_0, where
1184            //  v_s is the settling velocity of sediment grains
1185            //  p_0 is the mean precipitation rate
1186            //  d* is the sediment concentration ratio (between concentration near riverbed
1187            //  interface, and average concentration over the water column).
1188            //  d* varies with Rouse number which defines relative contribution of bed,
1189            // suspended,  and washed loads.
1190            //
1191            // G is typically on the order of 1 or greater.  However, we are only guaranteed
1192            // to converge for G ≤ 1, so we keep it in the chaos range of [0.12,
1193            // 1.32].
1194            1.0
1195        };
1196        let epsilon_0_func = |posi| {
1197            // epsilon_0_scale is roughly [using Hack's Law with b = 2 and SPL without
1198            // debris flow or hillslopes] equal to the ratio of the old to new
1199            // area, to the power of -n_i.
1200            let epsilon_0_scale_i = epsilon_0_scale(n_func(posi));
1201            if is_ocean_fn(posi) {
1202                // marine: ε₀ = 2.078e-3
1203                let epsilon_0_i = 2.078e-3 / 4.0;
1204                return epsilon_0_i * epsilon_0_scale_i;
1205            }
1206            let wposf = (uniform_idx_as_vec2(map_size_lg, posi)
1207                * TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
1208            .map(|e| e as f64);
1209            let turb_wposf = wposf
1210                .mul(5_000.0 / continent_scale)
1211                .div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64))
1212                .div(turb_wposf_div);
1213            let turb = Vec2::new(
1214                gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
1215                gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
1216            ) * uplift_turb_scale
1217                * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
1218            let turb_wposf = wposf + turb;
1219            let uheight = gen_ctx
1220                .uplift_nz
1221                .get(turb_wposf.into_array())
1222                .clamp(-1.0, 1.0)
1223                .mul(0.5)
1224                .add(0.5);
1225            let wposf3 = Vec3::new(
1226                wposf.x,
1227                wposf.y,
1228                uheight * CONFIG.mountain_scale as f64 * rock_strength_div_factor,
1229            );
1230            let rock_strength = gen_ctx
1231                .rock_strength_nz
1232                .get(wposf3.into_array())
1233                .clamp(-1.0, 1.0)
1234                .mul(0.5)
1235                .add(0.5);
1236            let center = 0.4;
1237            let dmin = center - 0.05;
1238            let dmax = center + 0.05;
1239            let log_odds = |x: f64| logit(x) - logit(center);
1240            let ustrength = logistic_cdf(
1241                1.0 * logit(rock_strength.clamp(1e-7, 1.0f64 - 1e-7))
1242                    + 1.0 * log_odds(uheight.clamp(dmin, dmax)),
1243            );
1244            // marine: ε₀ = 2.078e-3
1245            // San Gabriel Mountains: ε₀ = 3.18e-4
1246            // Oregon Coast Range: ε₀ = 2.68e-4
1247            // Frogs Hollow (peak production = 0.25): ε₀ = 1.41e-4
1248            // Point Reyes: ε₀ = 8.1e-5
1249            // Nunnock River (fractured granite, least weathered?): ε₀ = 5.3e-5
1250            let epsilon_0_i = ((1.0 - ustrength) * (2.078e-3 - 5.3e-5) + 5.3e-5) as f32 / 4.0;
1251            epsilon_0_i * epsilon_0_scale_i
1252        };
1253        let alpha_func = |posi| {
1254            let alpha_scale_i = alpha_scale(n_func(posi));
1255            if is_ocean_fn(posi) {
1256                // marine: α = 3.7e-2
1257                return 3.7e-2 * alpha_scale_i;
1258            }
1259            let wposf = (uniform_idx_as_vec2(map_size_lg, posi)
1260                * TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
1261            .map(|e| e as f64);
1262            let turb_wposf = wposf
1263                .mul(5_000.0 / continent_scale)
1264                .div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64))
1265                .div(turb_wposf_div);
1266            let turb = Vec2::new(
1267                gen_ctx.turb_x_nz.get(turb_wposf.into_array()),
1268                gen_ctx.turb_y_nz.get(turb_wposf.into_array()),
1269            ) * uplift_turb_scale
1270                * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
1271            let turb_wposf = wposf + turb;
1272            let uheight = gen_ctx
1273                .uplift_nz
1274                .get(turb_wposf.into_array())
1275                .clamp(-1.0, 1.0)
1276                .mul(0.5)
1277                .add(0.5);
1278            let wposf3 = Vec3::new(
1279                wposf.x,
1280                wposf.y,
1281                uheight * CONFIG.mountain_scale as f64 * rock_strength_div_factor,
1282            );
1283            let rock_strength = gen_ctx
1284                .rock_strength_nz
1285                .get(wposf3.into_array())
1286                .clamp(-1.0, 1.0)
1287                .mul(0.5)
1288                .add(0.5);
1289            let center = 0.4;
1290            let dmin = center - 0.05;
1291            let dmax = center + 0.05;
1292            let log_odds = |x: f64| logit(x) - logit(center);
1293            let ustrength = logistic_cdf(
1294                1.0 * logit(rock_strength.clamp(1e-7, 1.0f64 - 1e-7))
1295                    + 1.0 * log_odds(uheight.clamp(dmin, dmax)),
1296            );
1297            // Frog Hollow (peak production = 0.25): α = 4.2e-2
1298            // San Gabriel Mountains: α = 3.8e-2
1299            // marine: α = 3.7e-2
1300            // Oregon Coast Range: α = 3e-2
1301            // Nunnock river (fractured granite, least weathered?): α = 2e-3
1302            // Point Reyes: α = 1.6e-2
1303            // The stronger  the rock, the faster the decline in soil production.
1304            let alpha_i = (ustrength * (4.2e-2 - 1.6e-2) + 1.6e-2) as f32;
1305            alpha_i * alpha_scale_i
1306        };
1307        let uplift_fn = |posi| {
1308            if is_ocean_fn(posi) {
1309                return 0.0;
1310            }
1311            let height = (uplift_uniform[posi].1 - alt_old_min_uniform)
1312                / (alt_old_max_uniform - alt_old_min_uniform);
1313
1314            let height = height.mul(max_epsilon - min_epsilon).add(min_epsilon);
1315            let height = erosion_factor(height);
1316            assert!(height >= 0.0);
1317            assert!(height <= 1.0);
1318
1319            // u = 1e-3: normal-high (dike, mountain)
1320            // u = 5e-4: normal (mid example in Yuan, average mountain uplift)
1321            // u = 2e-4: low (low example in Yuan; known that lagoons etc. may have u ~
1322            // 0.05). u = 0: low (plateau [fan, altitude = 0.0])
1323
1324            height.mul(max_erosion_per_delta_t)
1325        };
1326        let alt_func = |posi| {
1327            if is_ocean_fn(posi) {
1328                old_height(posi)
1329            } else {
1330                (old_height(posi) as f64 / CONFIG.mountain_scale as f64) as f32 - 0.5
1331            }
1332        };
1333
1334        // Perform some erosion.
1335
1336        let mut last = None;
1337        let mut all_samples = std::time::Duration::default();
1338        let mut sample_count = 0;
1339        let report_erosion: &mut dyn FnMut(f64) = &mut move |progress: f64| {
1340            let now = std::time::Instant::now();
1341            let estimate = if let Some((last_instant, last_progress)) = last {
1342                if last_progress > progress {
1343                    None
1344                } else {
1345                    if last_progress < progress {
1346                        let sample = now
1347                            .duration_since(last_instant)
1348                            .div_f64(progress - last_progress);
1349                        all_samples += sample;
1350                        sample_count += 1;
1351                    }
1352
1353                    Some((all_samples / sample_count).mul_f64(100.0 - progress))
1354                }
1355            } else {
1356                None
1357            };
1358            last = Some((now, progress));
1359            stage_report(WorldSimStage::Erosion { progress, estimate })
1360        };
1361
1362        let (alt, basement) = if let Some(map) = parsed_world_file {
1363            (map.alt, map.basement)
1364        } else {
1365            let (alt, basement) = do_erosion(
1366                map_size_lg,
1367                max_erosion_per_delta_t as f32,
1368                n_steps,
1369                river_seed,
1370                // varying conditions
1371                &rock_strength_nz,
1372                // initial conditions
1373                alt_func,
1374                alt_func,
1375                is_ocean_fn,
1376                // empirical constants
1377                uplift_fn,
1378                n_func,
1379                theta_func,
1380                kf_func,
1381                kd_func,
1382                g_func,
1383                epsilon_0_func,
1384                alpha_func,
1385                // scaling factors
1386                height_scale,
1387                k_d_scale(n_approx),
1388                k_da_scale,
1389                threadpool,
1390                report_erosion,
1391            );
1392
1393            // Quick "small scale" erosion cycle in order to lower extreme angles.
1394            do_erosion(
1395                map_size_lg,
1396                1.0f32,
1397                n_small_steps,
1398                river_seed,
1399                &rock_strength_nz,
1400                |posi| alt[posi] as f32,
1401                |posi| basement[posi] as f32,
1402                is_ocean_fn,
1403                |posi| uplift_fn(posi) * (1.0 / max_erosion_per_delta_t),
1404                n_func,
1405                theta_func,
1406                kf_func,
1407                kd_func,
1408                g_func,
1409                epsilon_0_func,
1410                alpha_func,
1411                height_scale,
1412                k_d_scale(n_approx),
1413                k_da_scale,
1414                threadpool,
1415                report_erosion,
1416            )
1417        };
1418
1419        // Save map, if necessary.
1420        // NOTE: We wll always save a map with latest version.
1421        let map = WorldFile::new(ModernMap {
1422            continent_scale_hack: gen_opts.scale,
1423            map_size_lg: map_size_lg.vec(),
1424            alt,
1425            basement,
1426        });
1427        if fresh {
1428            world_file.save(&map);
1429        }
1430
1431        // Skip validation--we just performed a no-op conversion for this map, so it had
1432        // better be valid!
1433        let ModernMap {
1434            continent_scale_hack: _,
1435            map_size_lg: _,
1436            alt,
1437            basement,
1438        } = map.into_modern().unwrap();
1439
1440        // Additional small-scale erosion after map load, only used during testing.
1441        let (alt, basement) = if n_post_load_steps == 0 {
1442            (alt, basement)
1443        } else {
1444            do_erosion(
1445                map_size_lg,
1446                1.0f32,
1447                n_post_load_steps,
1448                river_seed,
1449                &rock_strength_nz,
1450                |posi| alt[posi] as f32,
1451                |posi| basement[posi] as f32,
1452                is_ocean_fn,
1453                |posi| uplift_fn(posi) * (1.0 / max_erosion_per_delta_t),
1454                n_func,
1455                theta_func,
1456                kf_func,
1457                kd_func,
1458                g_func,
1459                epsilon_0_func,
1460                alpha_func,
1461                height_scale,
1462                k_d_scale(n_approx),
1463                k_da_scale,
1464                threadpool,
1465                report_erosion,
1466            )
1467        };
1468
1469        let is_ocean = get_oceans(map_size_lg, |posi| alt[posi]);
1470        let is_ocean_fn = |posi: usize| is_ocean[posi];
1471        let mut dh = downhill(map_size_lg, |posi| alt[posi], is_ocean_fn);
1472        let (boundary_len, indirection, water_alt_pos, maxh) =
1473            get_lakes(map_size_lg, |posi| alt[posi], &mut dh);
1474        debug!(?maxh, "Max height");
1475        let (mrec, mstack, mwrec) = {
1476            let mut wh = vec![0.0; map_size_lg.chunks_len()];
1477            get_multi_rec(
1478                map_size_lg,
1479                |posi| alt[posi],
1480                &dh,
1481                &water_alt_pos,
1482                &mut wh,
1483                usize::from(map_size_lg.chunks().x),
1484                usize::from(map_size_lg.chunks().y),
1485                TerrainChunkSize::RECT_SIZE.x as Compute,
1486                TerrainChunkSize::RECT_SIZE.y as Compute,
1487                maxh,
1488                threadpool,
1489            )
1490        };
1491        let flux_old = get_multi_drainage(map_size_lg, &mstack, &mrec, &mwrec, boundary_len);
1492        // let flux_rivers = get_drainage(map_size_lg, &water_alt_pos, &dh,
1493        // boundary_len); TODO: Make rivers work with multi-direction flux as
1494        // well.
1495        let flux_rivers = flux_old.clone();
1496
1497        let water_height_initial = |chunk_idx| {
1498            let indirection_idx = indirection[chunk_idx];
1499            // Find the lake this point is flowing into.
1500            let lake_idx = if indirection_idx < 0 {
1501                chunk_idx
1502            } else {
1503                indirection_idx as usize
1504            };
1505            let chunk_water_alt = if dh[lake_idx] < 0 {
1506                // This is either a boundary node (dh[chunk_idx] == -2, i.e. water is at sea
1507                // level) or part of a lake that flows directly into the ocean.
1508                // In the former case, water is at sea level so we just return
1509                // 0.0.  In the latter case, the lake bottom must have been a
1510                // boundary node in the first place--meaning this node flows directly
1511                // into the ocean.  In that case, its lake bottom is ocean, meaning its water is
1512                // also at sea level.  Thus, we return 0.0 in both cases.
1513                0.0
1514            } else {
1515                // This chunk is draining into a body of water that isn't the ocean (i.e., a
1516                // lake). Then we just need to find the pass height of the
1517                // surrounding lake in order to figure out the initial water
1518                // height (which fill_sinks will then extend to make
1519                // sure it fills the entire basin).
1520
1521                // Find the height of "our" side of the pass (the part of it that drains into
1522                // this chunk's lake).
1523                let pass_idx = -indirection[lake_idx] as usize;
1524                let pass_height_i = alt[pass_idx];
1525                // Find the pass this lake is flowing into (i.e. water at the lake bottom gets
1526                // pushed towards the point identified by pass_idx).
1527                let neighbor_pass_idx = dh[pass_idx/*lake_idx*/];
1528                // Find the height of the pass into which our lake is flowing.
1529                let pass_height_j = alt[neighbor_pass_idx as usize];
1530                // Find the maximum of these two heights.
1531                // Use the pass height as the initial water altitude.
1532                pass_height_i.max(pass_height_j) /*pass_height*/
1533            };
1534            // Use the maximum of the pass height and chunk height as the parameter to
1535            // fill_sinks.
1536            let chunk_alt = alt[chunk_idx];
1537            chunk_alt.max(chunk_water_alt)
1538        };
1539
1540        // NOTE: If for for some reason you need to avoid the expensive `fill_sinks`
1541        // step here, and we haven't yet replaced it with a faster version, you
1542        // may comment out this line and replace it with the commented-out code
1543        // below; however, there are no guarantees that this
1544        // will work correctly.
1545        let water_alt = fill_sinks(map_size_lg, water_height_initial, is_ocean_fn);
1546        /* let water_alt = (0..map_size_lg.chunks_len())
1547        .into_par_iter()
1548        .map(|posi| water_height_initial(posi))
1549        .collect::<Vec<_>>(); */
1550
1551        let rivers = get_rivers(
1552            map_size_lg,
1553            gen_opts.scale,
1554            &water_alt_pos,
1555            &water_alt,
1556            &dh,
1557            &indirection,
1558            &flux_rivers,
1559        );
1560
1561        let water_alt = indirection
1562            .par_iter()
1563            .enumerate()
1564            .map(|(chunk_idx, &indirection_idx)| {
1565                // Find the lake this point is flowing into.
1566                let lake_idx = if indirection_idx < 0 {
1567                    chunk_idx
1568                } else {
1569                    indirection_idx as usize
1570                };
1571                if dh[lake_idx] < 0 {
1572                    // This is either a boundary node (dh[chunk_idx] == -2, i.e. water is at sea
1573                    // level) or part of a lake that flows directly into the
1574                    // ocean.  In the former case, water is at sea level so we
1575                    // just return 0.0.  In the latter case, the lake bottom must
1576                    // have been a boundary node in the first place--meaning this node flows
1577                    // directly into the ocean.  In that case, its lake bottom
1578                    // is ocean, meaning its water is also at sea level.  Thus,
1579                    // we return 0.0 in both cases.
1580                    0.0
1581                } else {
1582                    // This is not flowing into the ocean, so we can use the existing water_alt.
1583                    water_alt[chunk_idx] as f32
1584                }
1585            })
1586            .collect::<Vec<_>>()
1587            .into_boxed_slice();
1588
1589        let is_underwater = |chunk_idx: usize| match rivers[chunk_idx].river_kind {
1590            Some(RiverKind::Ocean) | Some(RiverKind::Lake { .. }) => true,
1591            Some(RiverKind::River { .. }) => false, // TODO: inspect width
1592            None => false,
1593        };
1594
1595        // Check whether any tiles around this tile are not water (since Lerp will
1596        // ensure that they are included).
1597        let pure_water = |posi: usize| {
1598            let pos = uniform_idx_as_vec2(map_size_lg, posi);
1599            for x in pos.x - 1..(pos.x + 1) + 1 {
1600                for y in pos.y - 1..(pos.y + 1) + 1 {
1601                    if x >= 0
1602                        && y >= 0
1603                        && x < map_size_lg.chunks().x as i32
1604                        && y < map_size_lg.chunks().y as i32
1605                    {
1606                        let posi = vec2_as_uniform_idx(map_size_lg, Vec2::new(x, y));
1607                        if !is_underwater(posi) {
1608                            return false;
1609                        }
1610                    }
1611                }
1612            }
1613            true
1614        };
1615
1616        // NaNs in these uniform vectors wherever pure_water() returns true.
1617        let (((alt_no_water, _), (pure_flux, _)), ((temp_base, _), (humid_base, _))) = threadpool
1618            .join(
1619                || {
1620                    threadpool.join(
1621                        || {
1622                            uniform_noise(map_size_lg, |posi, _| {
1623                                if pure_water(posi) {
1624                                    None
1625                                } else {
1626                                    // A version of alt that is uniform over *non-water* (or
1627                                    // land-adjacent water) chunks.
1628                                    Some(alt[posi] as f32)
1629                                }
1630                            })
1631                        },
1632                        || {
1633                            uniform_noise(map_size_lg, |posi, _| {
1634                                if pure_water(posi) {
1635                                    None
1636                                } else {
1637                                    Some(flux_old[posi])
1638                                }
1639                            })
1640                        },
1641                    )
1642                },
1643                || {
1644                    threadpool.join(
1645                        || {
1646                            uniform_noise(map_size_lg, |posi, wposf| {
1647                                if pure_water(posi) {
1648                                    None
1649                                } else {
1650                                    // -1 to 1.
1651                                    Some(gen_ctx.temp_nz.get((wposf).into_array()) as f32)
1652                                }
1653                            })
1654                        },
1655                        || {
1656                            uniform_noise(map_size_lg, |posi, wposf| {
1657                                // Check whether any tiles around this tile are water.
1658                                if pure_water(posi) {
1659                                    None
1660                                } else {
1661                                    // 0 to 1, hopefully.
1662                                    Some(
1663                                        (gen_ctx.humid_nz.get(wposf.div(1024.0).into_array())
1664                                            as f32)
1665                                            .add(1.0)
1666                                            .mul(0.5),
1667                                    )
1668                                }
1669                            })
1670                        },
1671                    )
1672                },
1673            );
1674
1675        let gen_cdf = GenCdf {
1676            humid_base,
1677            temp_base,
1678            chaos,
1679            alt,
1680            basement,
1681            water_alt,
1682            dh,
1683            flux: flux_old,
1684            pure_flux,
1685            alt_no_water,
1686            rivers,
1687        };
1688
1689        let chunks = (0..map_size_lg.chunks_len())
1690            .into_par_iter()
1691            .map(|i| SimChunk::generate(map_size_lg, i, &gen_ctx, &gen_cdf))
1692            .collect::<Vec<_>>();
1693
1694        let mut this = Self {
1695            seed,
1696            map_size_lg,
1697            max_height: maxh as f32,
1698            chunks,
1699            _locations: Vec::new(),
1700            gen_ctx,
1701            rng,
1702            calendar,
1703        };
1704
1705        this.generate_cliffs();
1706
1707        if opts.seed_elements {
1708            this.seed_elements();
1709        }
1710
1711        this
1712    }
1713
1714    #[inline(always)]
1715    pub const fn map_size_lg(&self) -> MapSizeLg { self.map_size_lg }
1716
1717    pub fn get_size(&self) -> Vec2<u32> { self.map_size_lg().chunks().map(u32::from) }
1718
1719    pub fn get_aabr(&self) -> Aabr<i32> {
1720        let size = self.get_size();
1721        Aabr {
1722            min: Vec2 { x: 0, y: 0 },
1723            max: Vec2 {
1724                x: size.x as i32,
1725                y: size.y as i32,
1726            },
1727        }
1728    }
1729
1730    pub fn generate_oob_chunk(&self) -> TerrainChunk {
1731        TerrainChunk::water(CONFIG.sea_level as i32)
1732    }
1733
1734    pub fn approx_chunk_terrain_normal(&self, chunk_pos: Vec2<i32>) -> Option<Vec3<f32>> {
1735        let curr_chunk = self.get(chunk_pos)?;
1736        let downhill_chunk_pos = curr_chunk.downhill?.wpos_to_cpos();
1737        let downhill_chunk = self.get(downhill_chunk_pos)?;
1738        // special case if chunks are flat
1739        if (curr_chunk.alt - downhill_chunk.alt) == 0. {
1740            return Some(Vec3::unit_z());
1741        }
1742        let curr = chunk_pos.cpos_to_wpos_center().as_().with_z(curr_chunk.alt);
1743        let down = downhill_chunk_pos
1744            .cpos_to_wpos_center()
1745            .as_()
1746            .with_z(downhill_chunk.alt);
1747        let downwards = curr - down;
1748        let flat = downwards.with_z(down.z);
1749        let mut res = downwards.cross(flat).cross(downwards);
1750        res.normalize();
1751        Some(res)
1752    }
1753
1754    /// Draw a map of the world based on chunk information.  Returns a buffer of
1755    /// u32s.
1756    pub fn get_map(&self, index: IndexRef, calendar: Option<&Calendar>) -> WorldMapMsg {
1757        prof_span!("WorldSim::get_map");
1758        let mut map_config = MapConfig::orthographic(
1759            self.map_size_lg(),
1760            core::ops::RangeInclusive::new(CONFIG.sea_level, CONFIG.sea_level + self.max_height),
1761        );
1762        // Build a horizon map.
1763        let scale_angle = |angle: Alt| {
1764            (/* 0.0.max( */angle /* ) */
1765                .atan()
1766                * <Alt as FloatConst>::FRAC_2_PI()
1767                * 255.0)
1768                .floor() as u8
1769        };
1770        let scale_height = |height: Alt| {
1771            (/* 0.0.max( */height/*)*/ as Alt * 255.0 / self.max_height as Alt).floor() as u8
1772        };
1773
1774        let samples_data = {
1775            prof_span!("samples data");
1776            let column_sample = ColumnGen::new(self);
1777            (0..self.map_size_lg().chunks_len())
1778                .into_par_iter()
1779                .map_init(
1780                    || Box::new(BlockGen::new(ColumnGen::new(self))),
1781                    |_block_gen, posi| {
1782                        let sample = column_sample.get(
1783                            (
1784                                uniform_idx_as_vec2(self.map_size_lg(), posi) * TerrainChunkSize::RECT_SIZE.map(|e| e as i32),
1785                                index,
1786                                calendar,
1787                            )
1788                        )?;
1789                        // sample.water_level = CONFIG.sea_level.max(sample.water_level);
1790
1791                        Some(sample)
1792                    },
1793                )
1794                /* .map(|posi| {
1795                    let mut sample = column_sample.get(
1796                        uniform_idx_as_vec2(self.map_size_lg(), posi) * TerrainChunkSize::RECT_SIZE.map(|e| e as i32),
1797                    );
1798                }) */
1799                .collect::<Vec<_>>()
1800                .into_boxed_slice()
1801        };
1802
1803        let horizons = get_horizon_map(
1804            self.map_size_lg(),
1805            Aabr {
1806                min: Vec2::zero(),
1807                max: self.map_size_lg().chunks().map(|e| e as i32),
1808            },
1809            CONFIG.sea_level,
1810            CONFIG.sea_level + self.max_height,
1811            |posi| {
1812                /* let chunk = &self.chunks[posi];
1813                chunk.alt.max(chunk.water_alt) as Alt */
1814                let sample = samples_data[posi].as_ref();
1815                sample
1816                    .map(|s| s.alt.max(s.water_level))
1817                    .unwrap_or(CONFIG.sea_level)
1818            },
1819            |a| scale_angle(a.into()),
1820            |h| scale_height(h.into()),
1821        )
1822        .unwrap();
1823
1824        let mut v = vec![0u32; self.map_size_lg().chunks_len()];
1825        let mut alts = vec![0u32; self.map_size_lg().chunks_len()];
1826        // TODO: Parallelize again.
1827        map_config.is_shaded = false;
1828
1829        map_config.generate(
1830            |pos| sample_pos(&map_config, self, index, Some(&samples_data), pos),
1831            |pos| sample_wpos(&map_config, self, pos),
1832            |pos, (r, g, b, _a)| {
1833                // We currently ignore alpha and replace it with the height at pos, scaled to
1834                // u8.
1835                let alt = sample_wpos(
1836                    &map_config,
1837                    self,
1838                    pos.map(|e| e as i32) * TerrainChunkSize::RECT_SIZE.map(|e| e as i32),
1839                );
1840                let a = 0; //(alt.min(1.0).max(0.0) * 255.0) as u8;
1841
1842                // NOTE: Safe by invariants on map_size_lg.
1843                let posi = (pos.y << self.map_size_lg().vec().x) | pos.x;
1844                v[posi] = u32::from_le_bytes([r, g, b, a]);
1845                alts[posi] = (((alt.clamp(0.0, 1.0) * 8191.0) as u32) & 0x1FFF) << 3;
1846            },
1847        );
1848        WorldMapMsg {
1849            dimensions_lg: self.map_size_lg().vec(),
1850            max_height: self.max_height,
1851            rgba: Grid::from_raw(self.get_size().map(|e| e as i32), v),
1852            alt: Grid::from_raw(self.get_size().map(|e| e as i32), alts),
1853            horizons,
1854            sites: Vec::new(),                   // Will be substituted later
1855            pois: Vec::new(),                    // Will be substituted later
1856            possible_starting_sites: Vec::new(), // Will be substituted later
1857            default_chunk: Arc::new(self.generate_oob_chunk()),
1858        }
1859    }
1860
1861    pub fn generate_cliffs(&mut self) {
1862        let mut rng = self.rng.clone();
1863
1864        for _ in 0..self.get_size().product() / 10 {
1865            let mut pos = self.get_size().map(|e| rng.gen_range(0..e) as i32);
1866
1867            let mut cliffs = DHashSet::default();
1868            let mut cliff_path = Vec::new();
1869
1870            for _ in 0..64 {
1871                if self.get_gradient_approx(pos).is_some_and(|g| g > 1.5) {
1872                    if !cliffs.insert(pos) {
1873                        break;
1874                    }
1875                    cliff_path.push((pos, 0.0));
1876
1877                    pos += CARDINALS
1878                        .iter()
1879                        .copied()
1880                        .max_by_key(|rpos| {
1881                            self.get_gradient_approx(pos + rpos)
1882                                .map_or(0, |g| (g * 1000.0) as i32)
1883                        })
1884                        .unwrap(); // Can't fail
1885                } else {
1886                    break;
1887                }
1888            }
1889
1890            for cliff in cliffs {
1891                Spiral2d::new()
1892                    .take((4usize * 2 + 1).pow(2))
1893                    .for_each(|rpos| {
1894                        let dist = rpos.map(|e| e as f32).magnitude();
1895                        if let Some(c) = self.get_mut(cliff + rpos) {
1896                            let warp = 1.0 / (1.0 + dist);
1897                            if !c.river.near_water() {
1898                                c.tree_density *= 1.0 - warp;
1899                                c.cliff_height = Lerp::lerp(44.0, 0.0, -1.0 + dist / 3.5);
1900                            }
1901                        }
1902                    });
1903            }
1904        }
1905    }
1906
1907    /// Prepare the world for simulation
1908    pub fn seed_elements(&mut self) {
1909        let mut rng = self.rng.clone();
1910
1911        let cell_size = 16;
1912        let grid_size = self.map_size_lg().chunks().map(usize::from) / cell_size;
1913        let loc_count = 100;
1914
1915        let mut loc_grid = vec![None; grid_size.product()];
1916        let mut locations = Vec::new();
1917
1918        // Seed the world with some locations
1919        (0..loc_count).for_each(|_| {
1920            let cell_pos = Vec2::new(
1921                self.rng.gen::<usize>() % grid_size.x,
1922                self.rng.gen::<usize>() % grid_size.y,
1923            );
1924            let wpos = (cell_pos * cell_size + cell_size / 2)
1925                .map2(TerrainChunkSize::RECT_SIZE, |e, sz: u32| {
1926                    e as i32 * sz as i32 + sz as i32 / 2
1927                });
1928
1929            locations.push(Location::generate(wpos, &mut rng));
1930
1931            loc_grid[cell_pos.y * grid_size.x + cell_pos.x] = Some(locations.len() - 1);
1932        });
1933
1934        // Find neighbours
1935        let mut loc_clone = locations
1936            .iter()
1937            .map(|l| l.center)
1938            .enumerate()
1939            .collect::<Vec<_>>();
1940        // NOTE: We assume that usize is 8 or fewer bytes.
1941        (0..locations.len()).for_each(|i| {
1942            let pos = locations[i].center.map(|e| e as i64);
1943
1944            loc_clone.sort_by_key(|(_, l)| l.map(|e| e as i64).distance_squared(pos));
1945
1946            loc_clone.iter().skip(1).take(2).for_each(|(j, _)| {
1947                locations[i].neighbours.insert(*j as u64);
1948                locations[*j].neighbours.insert(i as u64);
1949            });
1950        });
1951
1952        // Simulate invasion!
1953        let invasion_cycles = 25;
1954        (0..invasion_cycles).for_each(|_| {
1955            (0..grid_size.y).for_each(|j| {
1956                (0..grid_size.x).for_each(|i| {
1957                    if loc_grid[j * grid_size.x + i].is_none() {
1958                        const R_COORDS: [i32; 5] = [-1, 0, 1, 0, -1];
1959                        let idx = self.rng.gen::<usize>() % 4;
1960                        let new_i = i as i32 + R_COORDS[idx];
1961                        let new_j = j as i32 + R_COORDS[idx + 1];
1962                        if new_i >= 0 && new_j >= 0 {
1963                            let loc = Vec2::new(new_i as usize, new_j as usize);
1964                            loc_grid[j * grid_size.x + i] =
1965                                loc_grid.get(loc.y * grid_size.x + loc.x).cloned().flatten();
1966                        }
1967                    }
1968                });
1969            });
1970        });
1971
1972        // Place the locations onto the world
1973        /*
1974        let gen = StructureGen2d::new(self.seed, cell_size as u32, cell_size as u32 / 2);
1975
1976        self.chunks
1977            .par_iter_mut()
1978            .enumerate()
1979            .for_each(|(ij, chunk)| {
1980                let chunk_pos = uniform_idx_as_vec2(self.map_size_lg(), ij);
1981                let i = chunk_pos.x as usize;
1982                let j = chunk_pos.y as usize;
1983                let block_pos = Vec2::new(
1984                    chunk_pos.x * TerrainChunkSize::RECT_SIZE.x as i32,
1985                    chunk_pos.y * TerrainChunkSize::RECT_SIZE.y as i32,
1986                );
1987                let _cell_pos = Vec2::new(i / cell_size, j / cell_size);
1988
1989                // Find the distance to each region
1990                let near = gen.get(chunk_pos);
1991                let mut near = near
1992                    .iter()
1993                    .map(|(pos, seed)| RegionInfo {
1994                        chunk_pos: *pos,
1995                        block_pos: pos
1996                            .map2(TerrainChunkSize::RECT_SIZE, |e, sz: u32| e * sz as i32),
1997                        dist: (pos - chunk_pos).map(|e| e as f32).magnitude(),
1998                        seed: *seed,
1999                    })
2000                    .collect::<Vec<_>>();
2001
2002                // Sort regions based on distance
2003                near.sort_by(|a, b| a.dist.partial_cmp(&b.dist).unwrap());
2004
2005                let nearest_cell_pos = near[0].chunk_pos;
2006                if nearest_cell_pos.x >= 0 && nearest_cell_pos.y >= 0 {
2007                    let nearest_cell_pos = nearest_cell_pos.map(|e| e as usize) / cell_size;
2008                    chunk.location = loc_grid
2009                        .get(nearest_cell_pos.y * grid_size.x + nearest_cell_pos.x)
2010                        .cloned()
2011                        .unwrap_or(None)
2012                        .map(|loc_idx| LocationInfo { loc_idx, near });
2013                }
2014            });
2015        */
2016
2017        // Create waypoints
2018        const WAYPOINT_EVERY: usize = 16;
2019        let this = &self;
2020        let waypoints = (0..this.map_size_lg().chunks().x)
2021            .step_by(WAYPOINT_EVERY)
2022            .flat_map(|i| {
2023                (0..this.map_size_lg().chunks().y)
2024                    .step_by(WAYPOINT_EVERY)
2025                    .map(move |j| (i, j))
2026            })
2027            .collect::<Vec<_>>()
2028            .into_par_iter()
2029            .filter_map(|(i, j)| {
2030                let mut pos = Vec2::new(i as i32, j as i32);
2031                let mut chunk = this.get(pos)?;
2032
2033                if chunk.is_underwater() {
2034                    return None;
2035                }
2036                // Slide the waypoints down hills
2037                const MAX_ITERS: usize = 64;
2038                for _ in 0..MAX_ITERS {
2039                    let downhill_pos = match chunk.downhill {
2040                        Some(downhill) => {
2041                            downhill.map2(TerrainChunkSize::RECT_SIZE, |e, sz: u32| e / (sz as i32))
2042                        },
2043                        None => return Some(pos),
2044                    };
2045
2046                    let new_chunk = this.get(downhill_pos)?;
2047                    const SLIDE_THRESHOLD: f32 = 5.0;
2048                    if new_chunk.river.near_water() || new_chunk.alt + SLIDE_THRESHOLD < chunk.alt {
2049                        break;
2050                    } else {
2051                        chunk = new_chunk;
2052                        pos = downhill_pos;
2053                    }
2054                }
2055                Some(pos)
2056            })
2057            .collect::<Vec<_>>();
2058
2059        for waypoint in waypoints {
2060            self.get_mut(waypoint).map(|sc| sc.contains_waypoint = true);
2061        }
2062
2063        self.rng = rng;
2064        self._locations = locations;
2065    }
2066
2067    pub fn get(&self, chunk_pos: Vec2<i32>) -> Option<&SimChunk> {
2068        if chunk_pos
2069            .map2(self.map_size_lg().chunks(), |e, sz| e >= 0 && e < sz as i32)
2070            .reduce_and()
2071        {
2072            Some(&self.chunks[vec2_as_uniform_idx(self.map_size_lg(), chunk_pos)])
2073        } else {
2074            None
2075        }
2076    }
2077
2078    pub fn get_gradient_approx(&self, chunk_pos: Vec2<i32>) -> Option<f32> {
2079        let a = self.get(chunk_pos)?;
2080        if let Some(downhill) = a.downhill {
2081            let b = self.get(downhill.wpos_to_cpos())?;
2082            Some((a.alt - b.alt).abs() / TerrainChunkSize::RECT_SIZE.x as f32)
2083        } else {
2084            Some(0.0)
2085        }
2086    }
2087
2088    /// Get the altitude of the surface, could be water or ground.
2089    pub fn get_surface_alt_approx(&self, wpos: Vec2<i32>) -> f32 {
2090        self.get_interpolated(wpos, |chunk| chunk.alt)
2091            .zip(self.get_interpolated(wpos, |chunk| chunk.water_alt))
2092            .map(|(alt, water_alt)| alt.max(water_alt))
2093            .unwrap_or(CONFIG.sea_level)
2094    }
2095
2096    pub fn get_alt_approx(&self, wpos: Vec2<i32>) -> Option<f32> {
2097        self.get_interpolated(wpos, |chunk| chunk.alt)
2098    }
2099
2100    pub fn get_wpos(&self, wpos: Vec2<i32>) -> Option<&SimChunk> {
2101        self.get(wpos.map2(TerrainChunkSize::RECT_SIZE, |e, sz: u32| {
2102            e.div_euclid(sz as i32)
2103        }))
2104    }
2105
2106    pub fn get_mut(&mut self, chunk_pos: Vec2<i32>) -> Option<&mut SimChunk> {
2107        let map_size_lg = self.map_size_lg();
2108        if chunk_pos
2109            .map2(map_size_lg.chunks(), |e, sz| e >= 0 && e < sz as i32)
2110            .reduce_and()
2111        {
2112            Some(&mut self.chunks[vec2_as_uniform_idx(map_size_lg, chunk_pos)])
2113        } else {
2114            None
2115        }
2116    }
2117
2118    pub fn get_base_z(&self, chunk_pos: Vec2<i32>) -> Option<f32> {
2119        let in_bounds = chunk_pos
2120            .map2(self.map_size_lg().chunks(), |e, sz| {
2121                e > 0 && e < sz as i32 - 2
2122            })
2123            .reduce_and();
2124        if !in_bounds {
2125            return None;
2126        }
2127
2128        let chunk_idx = vec2_as_uniform_idx(self.map_size_lg(), chunk_pos);
2129        local_cells(self.map_size_lg(), chunk_idx)
2130            .flat_map(|neighbor_idx| {
2131                let neighbor_pos = uniform_idx_as_vec2(self.map_size_lg(), neighbor_idx);
2132                let neighbor_chunk = self.get(neighbor_pos);
2133                let river_kind = neighbor_chunk.and_then(|c| c.river.river_kind);
2134                let has_water = river_kind.is_some() && river_kind != Some(RiverKind::Ocean);
2135                if (neighbor_pos - chunk_pos).reduce_partial_max() <= 1 || has_water {
2136                    neighbor_chunk.map(|c| c.get_base_z())
2137                } else {
2138                    None
2139                }
2140            })
2141            .fold(None, |a: Option<f32>, x| a.map(|a| a.min(x)).or(Some(x)))
2142    }
2143
2144    pub fn get_interpolated<T, F>(&self, pos: Vec2<i32>, mut f: F) -> Option<T>
2145    where
2146        T: Copy + Default + Add<Output = T> + Mul<f32, Output = T>,
2147        F: FnMut(&SimChunk) -> T,
2148    {
2149        let pos = pos.as_::<f64>().wpos_to_cpos();
2150
2151        let cubic = |a: T, b: T, c: T, d: T, x: f32| -> T {
2152            let x2 = x * x;
2153
2154            // Catmull-Rom splines
2155            let co0 = a * -0.5 + b * 1.5 + c * -1.5 + d * 0.5;
2156            let co1 = a + b * -2.5 + c * 2.0 + d * -0.5;
2157            let co2 = a * -0.5 + c * 0.5;
2158            let co3 = b;
2159
2160            co0 * x2 * x + co1 * x2 + co2 * x + co3
2161        };
2162
2163        let mut x = [T::default(); 4];
2164
2165        for (x_idx, j) in (-1..3).enumerate() {
2166            let y0 = f(self.get(pos.map2(Vec2::new(j, -1), |e, q| e.max(0.0) as i32 + q))?);
2167            let y1 = f(self.get(pos.map2(Vec2::new(j, 0), |e, q| e.max(0.0) as i32 + q))?);
2168            let y2 = f(self.get(pos.map2(Vec2::new(j, 1), |e, q| e.max(0.0) as i32 + q))?);
2169            let y3 = f(self.get(pos.map2(Vec2::new(j, 2), |e, q| e.max(0.0) as i32 + q))?);
2170
2171            x[x_idx] = cubic(y0, y1, y2, y3, pos.y.fract() as f32);
2172        }
2173
2174        Some(cubic(x[0], x[1], x[2], x[3], pos.x.fract() as f32))
2175    }
2176
2177    /// M. Steffen splines.
2178    ///
2179    /// A more expensive cubic interpolation function that can preserve
2180    /// monotonicity between points.  This is useful if you rely on relative
2181    /// differences between endpoints being preserved at all interior
2182    /// points.  For example, we use this with riverbeds (and water
2183    /// height on along rivers) to maintain the invariant that the rivers always
2184    /// flow downhill at interior points (not just endpoints), without
2185    /// needing to flatten out the river.
2186    pub fn get_interpolated_monotone<T, F>(&self, pos: Vec2<i32>, mut f: F) -> Option<T>
2187    where
2188        T: Copy + Default + Signed + Float + Add<Output = T> + Mul<f32, Output = T>,
2189        F: FnMut(&SimChunk) -> T,
2190    {
2191        // See http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?1990A%26A...239..443S&defaultprint=YES&page_ind=0&filetype=.pdf
2192        //
2193        // Note that these are only guaranteed monotone in one dimension; fortunately,
2194        // that is sufficient for our purposes.
2195        let pos = pos.as_::<f64>().wpos_to_cpos();
2196
2197        let secant = |b: T, c: T| c - b;
2198
2199        let parabola = |a: T, c: T| -a * 0.5 + c * 0.5;
2200
2201        let slope = |_a: T, _b: T, _c: T, s_a: T, s_b: T, p_b: T| {
2202            // ((b - a).signum() + (c - b).signum()) * s
2203            (s_a.signum() + s_b.signum()) * (s_a.abs().min(s_b.abs()).min(p_b.abs() * 0.5))
2204        };
2205
2206        let cubic = |a: T, b: T, c: T, d: T, x: f32| -> T {
2207            // Compute secants.
2208            let s_a = secant(a, b);
2209            let s_b = secant(b, c);
2210            let s_c = secant(c, d);
2211            // Computing slopes from parabolas.
2212            let p_b = parabola(a, c);
2213            let p_c = parabola(b, d);
2214            // Get slopes (setting distance between neighbors to 1.0).
2215            let slope_b = slope(a, b, c, s_a, s_b, p_b);
2216            let slope_c = slope(b, c, d, s_b, s_c, p_c);
2217            let x2 = x * x;
2218
2219            // Interpolating splines.
2220            let co0 = slope_b + slope_c - s_b * 2.0;
2221            // = a * -0.5 + c * 0.5 + b * -0.5 + d * 0.5 - 2 * (c - b)
2222            // = a * -0.5 + b * 1.5 - c * 1.5 + d * 0.5;
2223            let co1 = s_b * 3.0 - slope_b * 2.0 - slope_c;
2224            // = (3.0 * (c - b) - 2.0 * (a * -0.5 + c * 0.5) - (b * -0.5 + d * 0.5))
2225            // = a + b * -2.5 + c * 2.0 + d * -0.5;
2226            let co2 = slope_b;
2227            // = a * -0.5 + c * 0.5;
2228            let co3 = b;
2229
2230            co0 * x2 * x + co1 * x2 + co2 * x + co3
2231        };
2232
2233        let mut x = [T::default(); 4];
2234
2235        for (x_idx, j) in (-1..3).enumerate() {
2236            let y0 = f(self.get(pos.map2(Vec2::new(j, -1), |e, q| e.max(0.0) as i32 + q))?);
2237            let y1 = f(self.get(pos.map2(Vec2::new(j, 0), |e, q| e.max(0.0) as i32 + q))?);
2238            let y2 = f(self.get(pos.map2(Vec2::new(j, 1), |e, q| e.max(0.0) as i32 + q))?);
2239            let y3 = f(self.get(pos.map2(Vec2::new(j, 2), |e, q| e.max(0.0) as i32 + q))?);
2240
2241            x[x_idx] = cubic(y0, y1, y2, y3, pos.y.fract() as f32);
2242        }
2243
2244        Some(cubic(x[0], x[1], x[2], x[3], pos.x.fract() as f32))
2245    }
2246
2247    /// Bilinear interpolation.
2248    ///
2249    /// Linear interpolation in both directions (i.e. quadratic interpolation).
2250    pub fn get_interpolated_bilinear<T, F>(&self, pos: Vec2<i32>, mut f: F) -> Option<T>
2251    where
2252        T: Copy + Default + Signed + Float + Add<Output = T> + Mul<f32, Output = T>,
2253        F: FnMut(&SimChunk) -> T,
2254    {
2255        // (i) Find downhill for all four points.
2256        // (ii) Compute distance from each downhill point and do linear interpolation on
2257        // their heights. (iii) Compute distance between each neighboring point
2258        // and do linear interpolation on       their distance-interpolated
2259        // heights.
2260
2261        // See http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?1990A%26A...239..443S&defaultprint=YES&page_ind=0&filetype=.pdf
2262        //
2263        // Note that these are only guaranteed monotone in one dimension; fortunately,
2264        // that is sufficient for our purposes.
2265        let pos = pos.as_::<f64>().wpos_to_cpos();
2266
2267        // Orient the chunk in the direction of the most downhill point of the four.  If
2268        // there is no "most downhill" point, then we don't care.
2269        let x0 = pos.map2(Vec2::new(0, 0), |e, q| e.max(0.0) as i32 + q);
2270        let p0 = self.get(x0)?;
2271        let y0 = f(p0);
2272
2273        let x1 = pos.map2(Vec2::new(1, 0), |e, q| e.max(0.0) as i32 + q);
2274        let p1 = self.get(x1)?;
2275        let y1 = f(p1);
2276
2277        let x2 = pos.map2(Vec2::new(0, 1), |e, q| e.max(0.0) as i32 + q);
2278        let p2 = self.get(x2)?;
2279        let y2 = f(p2);
2280
2281        let x3 = pos.map2(Vec2::new(1, 1), |e, q| e.max(0.0) as i32 + q);
2282        let p3 = self.get(x3)?;
2283        let y3 = f(p3);
2284
2285        let z0 = y0
2286            .mul(1.0 - pos.x.fract() as f32)
2287            .mul(1.0 - pos.y.fract() as f32);
2288        let z1 = y1.mul(pos.x.fract() as f32).mul(1.0 - pos.y.fract() as f32);
2289        let z2 = y2.mul(1.0 - pos.x.fract() as f32).mul(pos.y.fract() as f32);
2290        let z3 = y3.mul(pos.x.fract() as f32).mul(pos.y.fract() as f32);
2291
2292        Some(z0 + z1 + z2 + z3)
2293    }
2294
2295    #[expect(impl_trait_overcaptures)]
2296    pub fn get_nearest_ways<'a, M: Clone + Lerp<Output = M>>(
2297        &'a self,
2298        wpos: Vec2<i32>,
2299        get_way: &'a impl Fn(&SimChunk) -> Option<(Way, M)>,
2300    ) -> impl Iterator<Item = NearestWaysData<M, impl FnOnce() -> Vec2<f32>>> + 'a {
2301        let chunk_pos = wpos.map2(TerrainChunkSize::RECT_SIZE, |e, sz: u32| {
2302            e.div_euclid(sz as i32)
2303        });
2304        let get_chunk_centre = |chunk_pos: Vec2<i32>| {
2305            chunk_pos.map2(TerrainChunkSize::RECT_SIZE, |e, sz: u32| {
2306                e * sz as i32 + sz as i32 / 2
2307            })
2308        };
2309
2310        LOCALITY
2311            .iter()
2312            .filter_map(move |ctrl| {
2313                let (way, meta) = get_way(self.get(chunk_pos + *ctrl)?)?;
2314                let ctrl_pos = get_chunk_centre(chunk_pos + *ctrl).map(|e| e as f32)
2315                    + way.offset.map(|e| e as f32);
2316
2317                let chunk_connections = way.neighbors.count_ones();
2318                if chunk_connections == 0 {
2319                    return None;
2320                }
2321
2322                let (start_pos, start_idx, start_meta) = if chunk_connections != 2 {
2323                    (ctrl_pos, None, meta.clone())
2324                } else {
2325                    let (start_idx, start_rpos) = NEIGHBORS
2326                        .iter()
2327                        .copied()
2328                        .enumerate()
2329                        .find(|(i, _)| way.neighbors & (1 << *i as u8) != 0)
2330                        .unwrap();
2331                    let start_pos_chunk = chunk_pos + *ctrl + start_rpos;
2332                    let (start_way, start_meta) = get_way(self.get(start_pos_chunk)?)?;
2333                    (
2334                        get_chunk_centre(start_pos_chunk).map(|e| e as f32)
2335                            + start_way.offset.map(|e| e as f32),
2336                        Some(start_idx),
2337                        start_meta,
2338                    )
2339                };
2340
2341                Some(
2342                    NEIGHBORS
2343                        .iter()
2344                        .enumerate()
2345                        .filter(move |(i, _)| {
2346                            way.neighbors & (1 << *i as u8) != 0 && Some(*i) != start_idx
2347                        })
2348                        .filter_map(move |(i, end_rpos)| {
2349                            let end_pos_chunk = chunk_pos + *ctrl + end_rpos;
2350                            let (end_way, end_meta) = get_way(self.get(end_pos_chunk)?)?;
2351                            let end_pos = get_chunk_centre(end_pos_chunk).map(|e| e as f32)
2352                                + end_way.offset.map(|e| e as f32);
2353
2354                            let bez = QuadraticBezier2 {
2355                                start: (start_pos + ctrl_pos) / 2.0,
2356                                ctrl: ctrl_pos,
2357                                end: (end_pos + ctrl_pos) / 2.0,
2358                            };
2359                            let nearest_interval = bez
2360                                .binary_search_point_by_steps(wpos.map(|e| e as f32), 16, 0.001)
2361                                .0
2362                                .clamped(0.0, 1.0);
2363                            let pos = bez.evaluate(nearest_interval);
2364                            let dist_sqrd = pos.distance_squared(wpos.map(|e| e as f32));
2365                            let meta = if nearest_interval < 0.5 {
2366                                Lerp::lerp(start_meta.clone(), meta.clone(), 0.5 + nearest_interval)
2367                            } else {
2368                                Lerp::lerp(meta.clone(), end_meta, nearest_interval - 0.5)
2369                            };
2370                            Some(NearestWaysData {
2371                                i,
2372                                dist_sqrd,
2373                                pos,
2374                                meta,
2375                                bezier: bez,
2376                                calc_tangent: move || {
2377                                    bez.evaluate_derivative(nearest_interval).normalized()
2378                                },
2379                            })
2380                        }),
2381                )
2382            })
2383            .flatten()
2384    }
2385
2386    /// Return the distance to the nearest way in blocks, along with the
2387    /// closest point on the way, the way metadata, and the tangent vector
2388    /// of that way.
2389    pub fn get_nearest_way<M: Clone + Lerp<Output = M>>(
2390        &self,
2391        wpos: Vec2<i32>,
2392        get_way: impl Fn(&SimChunk) -> Option<(Way, M)>,
2393    ) -> Option<(f32, Vec2<f32>, M, Vec2<f32>)> {
2394        let get_way = &get_way;
2395        self.get_nearest_ways(wpos, get_way)
2396            .min_by_key(|NearestWaysData { dist_sqrd, .. }| (dist_sqrd * 1024.0) as i32)
2397            .map(
2398                |NearestWaysData {
2399                     dist_sqrd,
2400                     pos,
2401                     meta,
2402                     calc_tangent,
2403                     ..
2404                 }| (dist_sqrd.sqrt(), pos, meta, calc_tangent()),
2405            )
2406    }
2407
2408    pub fn get_nearest_path(&self, wpos: Vec2<i32>) -> Option<(f32, Vec2<f32>, Path, Vec2<f32>)> {
2409        self.get_nearest_way(wpos, |chunk| Some(chunk.path))
2410    }
2411
2412    /// Create a [`Lottery<Option<ForestKind>>`] that generates [`ForestKind`]s
2413    /// according to the conditions at the given position. If no or fewer
2414    /// trees are appropriate for the conditions, `None` may be generated.
2415    pub fn make_forest_lottery(&self, wpos: Vec2<i32>) -> Lottery<Option<ForestKind>> {
2416        let chunk = if let Some(chunk) = self.get_wpos(wpos) {
2417            chunk
2418        } else {
2419            return Lottery::from(vec![(1.0, None)]);
2420        };
2421        let env = chunk.get_environment();
2422        Lottery::from(
2423            ForestKind::iter()
2424                .enumerate()
2425                .map(|(i, fk)| {
2426                    const CLUSTER_SIZE: f64 = 48.0;
2427                    let nz = (FastNoise2d::new(i as u32 * 37)
2428                        .get(wpos.map(|e| e as f64) / CLUSTER_SIZE)
2429                        + 1.0)
2430                        / 2.0;
2431                    (fk.proclivity(&env) * nz, Some(fk))
2432                })
2433                .chain(std::iter::once((0.001, None)))
2434                .collect::<Vec<_>>(),
2435        )
2436    }
2437
2438    /// WARNING: Not currently used by the tree layer. Needs to be reworked.
2439    /// Return an iterator over candidate tree positions (note that only some of
2440    /// these will become trees since environmental parameters may forbid
2441    /// them spawning).
2442    pub fn get_near_trees(&self, wpos: Vec2<i32>) -> impl Iterator<Item = TreeAttr> + '_ {
2443        // Deterministic based on wpos
2444        self.gen_ctx
2445            .structure_gen
2446            .get(wpos)
2447            .into_iter()
2448            .filter_map(move |(wpos, seed)| {
2449                let lottery = self.make_forest_lottery(wpos);
2450                Some(TreeAttr {
2451                    pos: wpos,
2452                    seed,
2453                    scale: 1.0,
2454                    forest_kind: *lottery.choose_seeded(seed).as_ref()?,
2455                    inhabited: false,
2456                })
2457            })
2458    }
2459
2460    pub fn get_area_trees(
2461        &self,
2462        wpos_min: Vec2<i32>,
2463        wpos_max: Vec2<i32>,
2464    ) -> impl Iterator<Item = TreeAttr> + '_ {
2465        self.gen_ctx
2466            .structure_gen
2467            .iter(wpos_min, wpos_max)
2468            .filter_map(move |(wpos, seed)| {
2469                let lottery = self.make_forest_lottery(wpos);
2470                Some(TreeAttr {
2471                    pos: wpos,
2472                    seed,
2473                    scale: 1.0,
2474                    forest_kind: *lottery.choose_seeded(seed).as_ref()?,
2475                    inhabited: false,
2476                })
2477            })
2478    }
2479}
2480
2481#[derive(Debug)]
2482pub struct SimChunk {
2483    pub chaos: f32,
2484    pub alt: f32,
2485    pub basement: f32,
2486    pub water_alt: f32,
2487    pub downhill: Option<Vec2<i32>>,
2488    pub flux: f32,
2489    pub temp: f32,
2490    pub humidity: f32,
2491    pub rockiness: f32,
2492    pub tree_density: f32,
2493    pub forest_kind: ForestKind,
2494    pub spawn_rate: f32,
2495    pub river: RiverData,
2496    pub surface_veg: f32,
2497
2498    pub sites: Vec<Id<Site>>,
2499    pub place: Option<Id<Place>>,
2500    pub poi: Option<Id<PointOfInterest>>,
2501
2502    pub path: (Way, Path),
2503    pub cliff_height: f32,
2504    pub spot: Option<Spot>,
2505
2506    pub contains_waypoint: bool,
2507}
2508
2509#[derive(Copy, Clone)]
2510pub struct RegionInfo {
2511    pub chunk_pos: Vec2<i32>,
2512    pub block_pos: Vec2<i32>,
2513    pub dist: f32,
2514    pub seed: u32,
2515}
2516
2517pub struct NearestWaysData<M, F: FnOnce() -> Vec2<f32>> {
2518    pub i: usize,
2519    pub dist_sqrd: f32,
2520    pub pos: Vec2<f32>,
2521    pub meta: M,
2522    pub bezier: QuadraticBezier2<f32>,
2523    pub calc_tangent: F,
2524}
2525
2526impl SimChunk {
2527    fn generate(map_size_lg: MapSizeLg, posi: usize, gen_ctx: &GenCtx, gen_cdf: &GenCdf) -> Self {
2528        let pos = uniform_idx_as_vec2(map_size_lg, posi);
2529        let wposf = (pos * TerrainChunkSize::RECT_SIZE.map(|e| e as i32)).map(|e| e as f64);
2530
2531        let (_, chaos) = gen_cdf.chaos[posi];
2532        let alt_pre = gen_cdf.alt[posi] as f32;
2533        let basement_pre = gen_cdf.basement[posi] as f32;
2534        let water_alt_pre = gen_cdf.water_alt[posi];
2535        let downhill_pre = gen_cdf.dh[posi];
2536        let flux = gen_cdf.flux[posi] as f32;
2537        let river = gen_cdf.rivers[posi].clone();
2538
2539        // Can have NaNs in non-uniform part where pure_water returned true.  We just
2540        // test one of the four in order to find out whether this is the case.
2541        let (flux_uniform, /* flux_non_uniform */ _) = gen_cdf.pure_flux[posi];
2542        let (alt_uniform, _) = gen_cdf.alt_no_water[posi];
2543        let (temp_uniform, _) = gen_cdf.temp_base[posi];
2544        let (humid_uniform, _) = gen_cdf.humid_base[posi];
2545
2546        /* // Vertical difference from the equator (NOTE: "uniform" with much lower granularity than
2547        // other uniform quantities, but hopefully this doesn't matter *too* much--if it does, we
2548        // can always add a small x component).
2549        //
2550        // Not clear that we want this yet, let's see.
2551        let latitude_uniform = (pos.y as f32 / f32::from(self.map_size_lg().chunks().y)).sub(0.5).mul(2.0);
2552
2553        // Even less granular--if this matters we can make the sign affect the quantity slightly.
2554        let abs_lat_uniform = latitude_uniform.abs(); */
2555
2556        // We also correlate temperature negatively with altitude and absolute latitude,
2557        // using different weighting than we use for humidity.
2558        const TEMP_WEIGHTS: [f32; 3] = [/* 1.5, */ 1.0, 2.0, 1.0];
2559        let temp = cdf_irwin_hall(
2560            &TEMP_WEIGHTS,
2561            [
2562                temp_uniform,
2563                1.0 - alt_uniform, /* 1.0 - abs_lat_uniform*/
2564                (gen_ctx.rock_nz.get((wposf.div(50000.0)).into_array()) as f32 * 2.5 + 1.0) * 0.5,
2565            ],
2566        )
2567        // Convert to [-1, 1]
2568        .sub(0.5)
2569        .mul(2.0);
2570
2571        // Take the weighted average of our randomly generated base humidity, and the
2572        // calculated water flux over this point in order to compute humidity.
2573        const HUMID_WEIGHTS: [f32; 3] = [1.0, 1.0, 0.75];
2574        let humidity = cdf_irwin_hall(&HUMID_WEIGHTS, [humid_uniform, flux_uniform, 1.0]);
2575        // Moisture evaporates more in hot places
2576        let humidity = humidity
2577            * (1.0
2578                - (temp - CONFIG.tropical_temp)
2579                    .max(0.0)
2580                    .div(1.0 - CONFIG.tropical_temp))
2581            .max(0.0);
2582
2583        let mut alt = CONFIG.sea_level.add(alt_pre);
2584        let basement = CONFIG.sea_level.add(basement_pre);
2585        let water_alt = CONFIG.sea_level.add(water_alt_pre);
2586        let (downhill, _gradient) = if downhill_pre == -2 {
2587            (None, 0.0)
2588        } else if downhill_pre < 0 {
2589            panic!("Uh... shouldn't this never, ever happen?");
2590        } else {
2591            (
2592                Some(
2593                    uniform_idx_as_vec2(map_size_lg, downhill_pre as usize)
2594                        * TerrainChunkSize::RECT_SIZE.map(|e| e as i32)
2595                        + TerrainChunkSize::RECT_SIZE.map(|e| e as i32 / 2),
2596                ),
2597                (alt_pre - gen_cdf.alt[downhill_pre as usize] as f32).abs()
2598                    / TerrainChunkSize::RECT_SIZE.x as f32,
2599            )
2600        };
2601
2602        // Logistic regression.  Make sure x ∈ (0, 1).
2603        let logit = |x: f64| x.ln() - x.neg().ln_1p();
2604        // 0.5 + 0.5 * tanh(ln(1 / (1 - 0.1) - 1) / (2 * (sqrt(3)/pi)))
2605        let logistic_2_base = 3.0f64.sqrt().mul(std::f64::consts::FRAC_2_PI);
2606        // Assumes μ = 0, σ = 1
2607        let logistic_cdf = |x: f64| x.div(logistic_2_base).tanh().mul(0.5).add(0.5);
2608
2609        let is_underwater = match river.river_kind {
2610            Some(RiverKind::Ocean) | Some(RiverKind::Lake { .. }) => true,
2611            Some(RiverKind::River { .. }) => false, // TODO: inspect width
2612            None => false,
2613        };
2614        let river_xy = Vec2::new(river.velocity.x, river.velocity.y).magnitude();
2615        let river_slope = river.velocity.z / river_xy;
2616        match river.river_kind {
2617            Some(RiverKind::River { cross_section }) => {
2618                if cross_section.x >= 0.5 && cross_section.y >= CONFIG.river_min_height {
2619                    /* println!(
2620                        "Big area! Pos area: {:?}, River data: {:?}, slope: {:?}",
2621                        wposf, river, river_slope
2622                    ); */
2623                }
2624                if river_slope.abs() >= 0.25 && cross_section.x >= 1.0 {
2625                    let pos_area = wposf;
2626                    let river_data = &river;
2627                    debug!(?pos_area, ?river_data, ?river_slope, "Big waterfall!",);
2628                }
2629            },
2630            Some(RiverKind::Lake { .. }) => {
2631                // Forces lakes to be downhill from the land around them, and adds some noise to
2632                // the lake bed to make sure it's not too flat.
2633                let lake_bottom_nz = (gen_ctx.small_nz.get((wposf.div(20.0)).into_array()) as f32)
2634                    .clamp(-1.0, 1.0)
2635                    .mul(3.0);
2636                alt = alt.min(water_alt - 5.0) + lake_bottom_nz;
2637            },
2638            _ => {},
2639        }
2640
2641        // No trees in the ocean, with zero humidity (currently), or directly on
2642        // bedrock.
2643        let tree_density = if is_underwater {
2644            0.0
2645        } else {
2646            let tree_density = Lerp::lerp(
2647                -1.5,
2648                2.5,
2649                gen_ctx.tree_nz.get((wposf.div(1024.0)).into_array()) * 0.5 + 0.5,
2650            )
2651            .clamp(0.0, 1.0);
2652            // Tree density should go (by a lot) with humidity.
2653            if humidity <= 0.0 || tree_density <= 0.0 {
2654                0.0
2655            } else if humidity >= 1.0 || tree_density >= 1.0 {
2656                1.0
2657            } else {
2658                // Weighted logit sum.
2659                logistic_cdf(logit(tree_density))
2660            }
2661            // rescale to (-0.95, 0.95)
2662            .sub(0.5)
2663            .add(0.5)
2664        } as f32;
2665        const MIN_TREE_HUM: f32 = 0.15;
2666        let tree_density = tree_density
2667            // Tree density increases exponentially with humidity...
2668            .mul((humidity - MIN_TREE_HUM).max(0.0).mul(1.0 + MIN_TREE_HUM) / temp.max(0.75))
2669            // Places that are *too* wet (like marshes) also get fewer trees because the ground isn't stable enough for
2670            // them.
2671            //.mul((1.0 - flux * 0.05/*(humidity - 0.9).max(0.0) / 0.1*/).max(0.0))
2672            .mul(0.25 + flux * 0.05)
2673            // ...but is ultimately limited by available sunlight (and our tree generation system)
2674            .min(1.0);
2675
2676        // Add geologically short timescale undulation to the world for various reasons
2677        let alt =
2678            // Don't add undulation to rivers, mainly because this could accidentally result in rivers flowing uphill
2679            if river.near_water() {
2680                alt
2681            } else {
2682                // Sand dunes (formed over a short period of time, so we don't care about erosion sim)
2683                let warp = Vec2::new(
2684                    gen_ctx.turb_x_nz.get(wposf.div(350.0).into_array()) as f32,
2685                    gen_ctx.turb_y_nz.get(wposf.div(350.0).into_array()) as f32,
2686                ) * 200.0;
2687                const DUNE_SCALE: f32 = 24.0;
2688                const DUNE_LEN: f32 = 96.0;
2689                const DUNE_DIR: Vec2<f32> = Vec2::new(1.0, 1.0);
2690                let dune_dist = (wposf.map(|e| e as f32) + warp)
2691                    .div(DUNE_LEN)
2692                    .mul(DUNE_DIR.normalized())
2693                    .sum();
2694                let dune_nz = 0.5 - dune_dist.sin().abs() + 0.5 * (dune_dist + 0.5).sin().abs();
2695                let dune = dune_nz * DUNE_SCALE * (temp - 0.75).clamped(0.0, 0.25) * 4.0;
2696
2697                // Trees bind to soil and their roots result in small accumulating undulations over geologically short
2698                // periods of time. Forest floors are generally significantly bumpier than that of deforested areas.
2699                // This is particularly pronounced in high-humidity areas.
2700                let soil_nz = gen_ctx.hill_nz.get(wposf.div(96.0).into_array()) as f32;
2701                let soil_nz = (soil_nz + 1.0) * 0.5;
2702                const SOIL_SCALE: f32 = 16.0;
2703                let soil = soil_nz * SOIL_SCALE * tree_density.sqrt() * humidity.sqrt();
2704
2705                let warp_factor = ((alt - CONFIG.sea_level) / 16.0).clamped(0.0, 1.0);
2706
2707                let warp = (dune + soil) * warp_factor;
2708
2709                // Prevent warping pushing the altitude underwater
2710                if alt + warp < water_alt {
2711                    alt
2712                } else {
2713                    alt + warp
2714                }
2715            };
2716
2717        Self {
2718            chaos,
2719            flux,
2720            alt,
2721            basement: basement.min(alt),
2722            water_alt,
2723            downhill,
2724            temp,
2725            humidity,
2726            rockiness: if true {
2727                (gen_ctx.rock_nz.get((wposf.div(1024.0)).into_array()) as f32)
2728                    //.add(if river.near_river() { 20.0 } else { 0.0 })
2729                    .sub(0.1)
2730                    .mul(1.3)
2731                    .max(0.0)
2732            } else {
2733                0.0
2734            },
2735            tree_density,
2736            forest_kind: {
2737                let env = Environment {
2738                    humid: humidity,
2739                    temp,
2740                    near_water: if river.is_lake() || river.near_river() {
2741                        1.0
2742                    } else {
2743                        0.0
2744                    },
2745                };
2746
2747                ForestKind::iter()
2748                    .max_by_key(|fk| (fk.proclivity(&env) * 10000.0) as u32)
2749                    .unwrap() // Can't fail
2750            },
2751            spawn_rate: 1.0,
2752            river,
2753            surface_veg: 1.0,
2754
2755            sites: Vec::new(),
2756            place: None,
2757            poi: None,
2758            path: Default::default(),
2759            cliff_height: 0.0,
2760            spot: None,
2761
2762            contains_waypoint: false,
2763        }
2764    }
2765
2766    pub fn is_underwater(&self) -> bool {
2767        self.water_alt > self.alt || self.river.river_kind.is_some()
2768    }
2769
2770    pub fn get_base_z(&self) -> f32 { self.alt - self.chaos * 50.0 - 16.0 }
2771
2772    pub fn get_biome(&self) -> BiomeKind {
2773        let savannah_hum_temp = [0.05..0.55, 0.3..1.6];
2774        let taiga_hum_temp = [0.2..1.4, -0.7..-0.3];
2775        if self.river.is_ocean() {
2776            BiomeKind::Ocean
2777        } else if self.river.is_lake() {
2778            BiomeKind::Lake
2779        } else if self.temp < CONFIG.snow_temp {
2780            BiomeKind::Snowland
2781        } else if self.alt > 500.0 && self.chaos > 0.3 && self.tree_density < 0.6 {
2782            BiomeKind::Mountain
2783        } else if self.temp > CONFIG.desert_temp && self.humidity < CONFIG.desert_hum {
2784            BiomeKind::Desert
2785        } else if self.tree_density > 0.65 && self.humidity > 0.65 && self.temp > 0.45 {
2786            BiomeKind::Jungle
2787        } else if savannah_hum_temp[0].contains(&self.humidity)
2788            && savannah_hum_temp[1].contains(&self.temp)
2789        {
2790            BiomeKind::Savannah
2791        } else if taiga_hum_temp[0].contains(&self.humidity)
2792            && taiga_hum_temp[1].contains(&self.temp)
2793        {
2794            BiomeKind::Taiga
2795        } else if self.tree_density > 0.4 {
2796            BiomeKind::Forest
2797        // } else if self.humidity > 0.8 {
2798        //    BiomeKind::Swamp
2799        //      Swamps don't really exist yet.
2800        } else {
2801            BiomeKind::Grassland
2802        }
2803    }
2804
2805    pub fn near_cliffs(&self) -> bool { self.cliff_height > 0.0 }
2806
2807    pub fn get_environment(&self) -> Environment {
2808        Environment {
2809            humid: self.humidity,
2810            temp: self.temp,
2811            near_water: if self.river.is_lake()
2812                || self.river.near_river()
2813                || self.alt < CONFIG.sea_level + 6.0
2814            // Close to sea in altitude
2815            {
2816                1.0
2817            } else {
2818                0.0
2819            },
2820        }
2821    }
2822
2823    pub fn get_location_name(
2824        &self,
2825        index_sites: &Store<crate::site::Site>,
2826        civs_pois: &Store<PointOfInterest>,
2827        wpos2d: Vec2<i32>,
2828    ) -> Option<String> {
2829        self.sites
2830            .iter()
2831            .filter(|id| {
2832                index_sites[**id].get_origin().distance_squared(wpos2d) as f32
2833                    <= index_sites[**id].radius().powi(2)
2834            })
2835            .min_by_key(|id| index_sites[**id].get_origin().distance_squared(wpos2d))
2836            .map(|id| index_sites[*id].name().to_string())
2837            .or_else(|| self.poi.map(|poi| civs_pois[poi].name.clone()))
2838    }
2839}