mirror of
https://gitlab.com/veloren/veloren.git
synced 2024-08-30 18:12:32 +00:00
Adding hillslope diffusion.
This commit is contained in:
Joshua Yanovski
parent
825ba1b4f4
commit
69a8d5039b
@ -990,6 +990,7 @@ fn handle_debug_column(server: &mut Server, entity: EcsEntity, args: String, act
|
||||
let foo = || {
|
||||
// let sim_chunk = sim.get(chunk_pos)?;
|
||||
let alt = sim.get_interpolated(wpos, |chunk| chunk.alt)?;
|
||||
let basement = sim.get_interpolated(wpos, |chunk| chunk.basement)?;
|
||||
let water_alt = sim.get_interpolated(wpos, |chunk| chunk.water_alt)?;
|
||||
let chaos = sim.get_interpolated(wpos, |chunk| chunk.chaos)?;
|
||||
let temp = sim.get_interpolated(wpos, |chunk| chunk.temp)?;
|
||||
@ -1008,6 +1009,7 @@ fn handle_debug_column(server: &mut Server, entity: EcsEntity, args: String, act
|
||||
r#"wpos: {:?}
|
||||
alt {:?} ({:?})
|
||||
water_alt {:?} ({:?})
|
||||
basement {:?}
|
||||
river {:?}
|
||||
downhill {:?}
|
||||
chaos {:?}
|
||||
@ -1022,6 +1024,7 @@ spawn_rate {:?} "#,
|
||||
col.alt,
|
||||
water_alt,
|
||||
col.water_level,
|
||||
basement,
|
||||
river,
|
||||
downhill,
|
||||
chaos,
|
||||
|
||||
@ -26,6 +26,8 @@ fn main() {
|
||||
let light_direction = Vec3::new(-0.8, -1.0, 0.3).normalized();
|
||||
let light_res = 3;
|
||||
|
||||
let mut is_basement = false;
|
||||
|
||||
while win.is_open() {
|
||||
let mut buf = vec![0; W * H];
|
||||
const QUADRANTS: usize = 4;
|
||||
@ -34,6 +36,8 @@ fn main() {
|
||||
let mut lakes = 0u32;
|
||||
let mut oceans = 0u32;
|
||||
|
||||
// let water_light = (light_direction.z + 1.0) / 2.0 * 0.8 + 0.2;
|
||||
|
||||
for i in 0..W {
|
||||
for j in 0..H {
|
||||
let pos = focus + Vec2::new(i as i32, j as i32) * scale;
|
||||
@ -41,11 +45,12 @@ fn main() {
|
||||
let top_right = focus + Vec2::new(i as i32 + light_res, j as i32) * scale;
|
||||
let bottom_left = focus + Vec2::new(i as i32, j as i32 + light_res) * scale; */
|
||||
|
||||
let (alt, water_alt, humidity, temperature, downhill, river_kind) = sampler
|
||||
let (alt, basement, water_alt, humidity, temperature, downhill, river_kind) = sampler
|
||||
.get(pos)
|
||||
.map(|sample| {
|
||||
(
|
||||
sample.alt,
|
||||
sample.basement,
|
||||
sample.water_alt,
|
||||
sample.humidity,
|
||||
sample.temp,
|
||||
@ -53,7 +58,7 @@ fn main() {
|
||||
sample.river.river_kind,
|
||||
)
|
||||
})
|
||||
.unwrap_or((CONFIG.sea_level, CONFIG.sea_level, 0.0, 0.0, None, None));
|
||||
.unwrap_or((CONFIG.sea_level, CONFIG.sea_level, CONFIG.sea_level, 0.0, 0.0, None, None));
|
||||
let humidity = humidity.min(1.0).max(0.0);
|
||||
let temperature = temperature.min(1.0).max(-1.0) * 0.5 + 0.5;
|
||||
let downhill_pos = (downhill
|
||||
@ -63,8 +68,9 @@ fn main() {
|
||||
+ pos;
|
||||
let downhill_alt = sampler
|
||||
.get(downhill_pos)
|
||||
.map(|s| s.alt)
|
||||
.map(|s| if is_basement { s.basement } else { s.alt })
|
||||
.unwrap_or(CONFIG.sea_level);
|
||||
let alt = if is_basement { basement } else { alt };
|
||||
/* let alt_tl = sampler.get(top_left).map(|s| s.alt)
|
||||
.unwrap_or(CONFIG.sea_level);
|
||||
let alt_tr = sampler.get(top_right).map(|s| s.alt)
|
||||
@ -78,7 +84,7 @@ fn main() {
|
||||
.map(|e| e as i32));
|
||||
let cross_alt = sampler
|
||||
.get(cross_pos)
|
||||
.map(|s| s.alt)
|
||||
.map(|s| if is_basement { s.basement } else { s.alt })
|
||||
.unwrap_or(CONFIG.sea_level);
|
||||
let forward_vec = Vec3::new(
|
||||
(downhill_pos.x - pos.x) as f64,
|
||||
@ -98,7 +104,11 @@ fn main() {
|
||||
let water_alt = ((alt.max(water_alt) - CONFIG.sea_level) as f64 / gain as f64)
|
||||
.min(1.0)
|
||||
.max(0.0);
|
||||
let alt = ((alt - CONFIG.sea_level) as f64 / gain as f64)
|
||||
let true_alt = (alt - CONFIG.sea_level) as f64 / gain as f64;
|
||||
let water_depth = (water_alt - true_alt)
|
||||
.min(1.0)
|
||||
.max(0.0);
|
||||
let alt = true_alt
|
||||
.min(1.0)
|
||||
.max(0.0);
|
||||
let quad =
|
||||
@ -120,10 +130,15 @@ fn main() {
|
||||
}
|
||||
|
||||
buf[j * W + i] = match river_kind {
|
||||
Some(RiverKind::Ocean) => u32::from_le_bytes([64, 32, 0, 255]),
|
||||
Some(RiverKind::Ocean) => u32::from_le_bytes([
|
||||
((64.0 - water_depth * 64.0) * 1.0) as u8,
|
||||
((32.0 - water_depth * 32.0) * 1.0) as u8,
|
||||
0,
|
||||
255,
|
||||
]),
|
||||
Some(RiverKind::Lake { .. }) => u32::from_le_bytes([
|
||||
64 + (water_alt * 191.0) as u8,
|
||||
32 + (water_alt * 95.0) as u8,
|
||||
(((64.0 + water_alt * 191.0) + (- water_depth * 64.0)) * 1.0) as u8,
|
||||
(((32.0 + water_alt * 95.0) + (- water_depth * 32.0)) * 1.0) as u8,
|
||||
0,
|
||||
255,
|
||||
]),
|
||||
@ -183,6 +198,9 @@ fn main() {
|
||||
);
|
||||
}
|
||||
}
|
||||
if win.is_key_down(minifb::Key::B) {
|
||||
is_basement ^= true;
|
||||
}
|
||||
if win.is_key_down(minifb::Key::W) {
|
||||
focus.y -= spd * scale;
|
||||
}
|
||||
|
||||
@ -436,9 +436,12 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
|
||||
(temp < CONFIG.snow_temp ||
|
||||
downhill_alt.sub(CONFIG.sea_level) >= CONFIG.mountain_scale * 0.25)*/
|
||||
is_rocky {
|
||||
sim.get_interpolated_bilinear(wpos, |chunk| chunk.alt)?
|
||||
sim.get_interpolated_monotone(wpos, |chunk| chunk.alt)?
|
||||
// sim.get_interpolated_bilinear(wpos, |chunk| chunk.alt)?
|
||||
// sim.get_interpolated(wpos, |chunk| chunk.alt)?
|
||||
} else {
|
||||
sim.get_interpolated_monotone(wpos, |chunk| chunk.alt)?
|
||||
// sim.get_interpolated(wpos, |chunk| chunk.alt)?
|
||||
};
|
||||
|
||||
// Find the average distance to each neighboring body of water.
|
||||
@ -581,23 +584,6 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
|
||||
.get((wposf_turb.div(128.0)).into_array()) as f32)
|
||||
.mul(24.0);
|
||||
|
||||
let riverless_alt_delta = (sim
|
||||
.gen_ctx
|
||||
.small_nz
|
||||
.get((wposf_turb.div(/*200.0*//*50.0*//*24.0*//*56.0 / (chaos as f64).max(0.05)*/50.0)).into_array()) as f32)
|
||||
.abs()
|
||||
.mul(3.0)
|
||||
/* .mul(chaos.max(0.05))
|
||||
.mul(27.0) */
|
||||
/* + (sim
|
||||
.gen_ctx
|
||||
.small_nz
|
||||
.get((wposf_turb.div(400.0)).into_array()) as f32)
|
||||
.abs()
|
||||
.mul((1.0 - chaos).max(0.3))
|
||||
.mul(1.0 - humidity)
|
||||
.mul(32.0) */;
|
||||
|
||||
let alt_for_river = alt
|
||||
+ if overlap_count == 0.0 {
|
||||
0.0
|
||||
@ -818,6 +804,28 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
|
||||
(false, alt_for_river, downhill_water_alt, 1.0)
|
||||
};
|
||||
|
||||
let riverless_alt_delta = (sim
|
||||
.gen_ctx
|
||||
.small_nz
|
||||
.get((wposf_turb.div(/*200.0*/50.0/*24.0*//*56.0 / (chaos as f64).max(0.05)*//*50.0*/)).into_array()) as f32)
|
||||
.min(1.0).max(-1.0)
|
||||
// .mul(0.5).add(0.5)
|
||||
.abs()
|
||||
.mul(3.0)
|
||||
/* .mul(chaos.max(0.05))
|
||||
.mul(27.0)
|
||||
+ (sim
|
||||
.gen_ctx
|
||||
.small_nz
|
||||
.get((wposf_turb.div(400.0)).into_array()) as f32)
|
||||
.min(1.0).max(-1.0)
|
||||
// .mul(0.5).add(0.5)
|
||||
.abs()
|
||||
.mul((1.0 - chaos).max(0.3))
|
||||
.mul(1.0 - humidity)
|
||||
.mul(32.0) */;
|
||||
|
||||
|
||||
let riverless_alt_delta = Lerp::lerp(0.0, riverless_alt_delta, warp_factor);
|
||||
let alt = alt_ + riverless_alt_delta;
|
||||
|
||||
|
||||
421
world/src/sim/diffusion.rs
Normal file
421
world/src/sim/diffusion.rs
Normal file
@ -0,0 +1,421 @@
|
||||
use rayon::prelude::*;
|
||||
|
||||
/// From https://github.com/fastscape-lem/fastscapelib-fortran/blob/master/src/Diffusion.f90
|
||||
///
|
||||
/// See https://fastscape-lem.github.io/fastscapelib-fortran
|
||||
///
|
||||
/// nx = grid size x
|
||||
///
|
||||
/// ny = grid size y
|
||||
///
|
||||
/// xl = x-dimension of the model topography in meters (double precision)
|
||||
///
|
||||
/// yl = y-dimension of the model topography in meters (double precision)
|
||||
///
|
||||
/// dt = length of the time step in years (double precision)
|
||||
///
|
||||
/// ibc = boundary conditions for grid. For now, we assume all four boundaries are fixed height,
|
||||
/// so this parameter is ignored.
|
||||
///
|
||||
/// h = heights of cells at each cell in the grid.
|
||||
///
|
||||
/// b = basement height at each cell in the grid (see https://en.wikipedia.org/wiki/Basement_(geology)).
|
||||
///
|
||||
/// kd = bedrock transport coefficient (or diffusivity) for hillslope processes in meter squared per year
|
||||
/// (double precision) at each cell in the grid.
|
||||
///
|
||||
/// kdsed = sediment transport coefficient (or diffusivity) for hillslope processes in meter squared per year;
|
||||
/// (double precision;) note that when kdsed < 0, its value is not used, i.e., kd for sediment and bedrock have
|
||||
/// the same value, regardless of sediment thickness
|
||||
/* subroutine Diffusion ()
|
||||
|
||||
! subroutine to solve the diffusion equation by ADI
|
||||
|
||||
use FastScapeContext
|
||||
|
||||
implicit none
|
||||
*/
|
||||
pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (),
|
||||
h: &mut [/*f64*/f32], b: &mut [/*f64*/f32],
|
||||
kd: impl Fn(usize) -> f64, kdsed: f64,
|
||||
) {
|
||||
let aij = |i: usize, j: usize| j * nx + i;
|
||||
/*
|
||||
double precision, dimension(:), allocatable :: f,diag,sup,inf,res
|
||||
double precision, dimension(:,:), allocatable :: zint,kdint,zintp
|
||||
integer i,j,ij
|
||||
double precision factxp,factxm,factyp,factym,dx,dy
|
||||
*/
|
||||
let mut f : Vec<f64>;
|
||||
let mut diag : Vec<f64>;
|
||||
let mut sup : Vec<f64>;
|
||||
let mut inf : Vec<f64>;
|
||||
let mut res : Vec<f64>;
|
||||
let mut zint : Vec<f64>;
|
||||
let mut kdint : Vec<f64>;
|
||||
let mut zintp : Vec<f64>;
|
||||
let mut i : usize;
|
||||
let mut j : usize;
|
||||
let mut ij : usize;
|
||||
let mut factxp : f64;
|
||||
let mut factxm : f64;
|
||||
let mut factyp : f64;
|
||||
let mut factym : f64;
|
||||
let mut dx : f64;
|
||||
let mut dy : f64;
|
||||
/*
|
||||
character cbc*4
|
||||
|
||||
!print*,'Diffusion'
|
||||
|
||||
write (cbc,'(i4)') ibc
|
||||
|
||||
dx=xl/(nx-1)
|
||||
dy=yl/(ny-1)
|
||||
*/
|
||||
// 2048*32/2048/2048
|
||||
// 1 / 64 m
|
||||
dx = xl / /*(nx - 1)*/nx as f64;
|
||||
dy = yl / /*(ny - 1)*/ny as f64;
|
||||
/*
|
||||
! creates 2D internal arrays to store topo and kd
|
||||
|
||||
allocate (zint(nx,ny),kdint(nx,ny),zintp(nx,ny))
|
||||
*/
|
||||
zint = vec![Default::default(); nx * ny];
|
||||
kdint = vec![Default::default(); nx * ny];
|
||||
zintp = vec![Default::default(); nx * ny];
|
||||
/*
|
||||
do j=1,ny
|
||||
do i=1,nx
|
||||
ij=(j-1)*nx+i
|
||||
zint(i,j)=h(ij)
|
||||
kdint(i,j)=kd(ij)
|
||||
if (kdsed.gt.0.d0 .and. (h(ij)-b(ij)).gt.1.d-6) kdint(i,j)=kdsed
|
||||
enddo
|
||||
enddo
|
||||
|
||||
zintp = zint
|
||||
*/
|
||||
for j in 0..ny {
|
||||
for i in 0..nx {
|
||||
// ij = vec2_as_uniform_idx(i, j);
|
||||
ij = aij(i, j);
|
||||
zint[ij] = h[ij] as f64;
|
||||
kdint[ij] = kd(ij);
|
||||
if kdsed > 0.0 && (h[ij] - b[ij]) > 1.0e-6 {
|
||||
kdint[ij] = kdsed;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
zintp = zint.clone();
|
||||
/*
|
||||
! first pass along the x-axis
|
||||
|
||||
allocate (f(nx),diag(nx),sup(nx),inf(nx),res(nx))
|
||||
f=0.d0
|
||||
diag=0.d0
|
||||
sup=0.d0
|
||||
inf=0.d0
|
||||
res=0.d0
|
||||
do j=2,ny-1
|
||||
*/
|
||||
f = vec![0.0; nx];
|
||||
diag = vec![0.0; nx];
|
||||
sup = vec![0.0; nx];
|
||||
inf = vec![0.0; nx];
|
||||
res = vec![0.0; nx];
|
||||
for j in 1..ny-1 {
|
||||
/*
|
||||
do i=2,nx-1
|
||||
factxp=(kdint(i+1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2
|
||||
factxm=(kdint(i-1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2
|
||||
factyp=(kdint(i,j+1)+kdint(i,j))/2.d0*(dt/2.)/dy**2
|
||||
factym=(kdint(i,j-1)+kdint(i,j))/2.d0*(dt/2.)/dy**2
|
||||
diag(i)=1.d0+factxp+factxm
|
||||
sup(i)=-factxp
|
||||
inf(i)=-factxm
|
||||
f(i)=zintp(i,j)+factyp*zintp(i,j+1)-(factyp+factym)*zintp(i,j)+factym*zintp(i,j-1)
|
||||
enddo
|
||||
*/
|
||||
for i in 1..nx-1 {
|
||||
factxp = (kdint[aij(i+1, j)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dx * dx);
|
||||
factxm = (kdint[aij(i-1, j)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dx * dx);
|
||||
factyp = (kdint[aij(i, j+1)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dy * dy);
|
||||
factym = (kdint[aij(i, j-1)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dy * dy);
|
||||
diag[i] = 1.0 + factxp + factxm;
|
||||
sup[i] = -factxp;
|
||||
inf[i] = -factxm;
|
||||
f[i] =
|
||||
zintp[aij(i, j)] + factyp * zintp[aij(i, j+1)] -
|
||||
(factyp + factym) *
|
||||
zintp[aij(i, j)] + factym * zintp[aij(i, j-1)];
|
||||
}
|
||||
/*
|
||||
! left bc
|
||||
if (cbc(4:4).eq.'1') then
|
||||
diag(1)=1.
|
||||
sup(1)=0.
|
||||
f(1)=zintp(1,j)
|
||||
else
|
||||
factxp=(kdint(2,j)+kdint(1,j))/2.d0*(dt/2.)/dx**2
|
||||
factyp=(kdint(1,j+1)+kdint(1,j))/2.d0*(dt/2.)/dy**2
|
||||
factym=(kdint(1,j-1)+kdint(1,j))/2.d0*(dt/2.)/dy**2
|
||||
diag(1)=1.d0+factxp
|
||||
sup(1)=-factxp
|
||||
f(1)=zintp(1,j)+factyp*zintp(1,j+1)-(factyp+factym)*zintp(1,j)+factym*zintp(1,j-1)
|
||||
endif
|
||||
*/
|
||||
if true {
|
||||
diag[0] = 1.0;
|
||||
sup[0] = 0.0;
|
||||
f[0] = zintp[aij(0, j)];
|
||||
} else {
|
||||
// FIXME: implement reflective boundary
|
||||
}
|
||||
/*
|
||||
! right bc
|
||||
if (cbc(2:2).eq.'1') then
|
||||
diag(nx)=1.
|
||||
inf(nx)=0.
|
||||
f(nx)=zintp(nx,j)
|
||||
else
|
||||
factxm=(kdint(nx-1,j)+kdint(nx,j))/2.d0*(dt/2.)/dx**2
|
||||
factyp=(kdint(nx,j+1)+kdint(nx,j))/2.d0*(dt/2.)/dy**2
|
||||
factym=(kdint(nx,j-1)+kdint(nx,j))/2.d0*(dt/2.)/dy**2
|
||||
diag(nx)=1.d0+factxm
|
||||
inf(nx)=-factxm
|
||||
f(nx)=zintp(nx,j)+factyp*zintp(nx,j+1)-(factyp+factym)*zintp(nx,j)+factym*zintp(nx,j-1)
|
||||
endif
|
||||
*/
|
||||
if true {
|
||||
diag[nx - 1] = 1.0;
|
||||
inf[nx - 1] = 0.0;
|
||||
f[nx - 1] = zintp[aij(nx - 1, j)];
|
||||
} else {
|
||||
// FIXME: implement reflective boundary
|
||||
}
|
||||
/*
|
||||
call tridag (inf,diag,sup,f,res,nx)
|
||||
do i=1,nx
|
||||
zint(i,j)=res(i)
|
||||
enddo
|
||||
*/
|
||||
tridag(&inf, &diag, &sup, &f, &mut res, nx);
|
||||
for i in 0..nx {
|
||||
zint[aij(i, j)] = res[i];
|
||||
}
|
||||
/*
|
||||
enddo
|
||||
deallocate (f,diag,sup,inf,res)
|
||||
*/
|
||||
}
|
||||
|
||||
/*
|
||||
! second pass along y-axis
|
||||
|
||||
allocate (f(ny),diag(ny),sup(ny),inf(ny),res(ny))
|
||||
f=0.d0
|
||||
diag=0.d0
|
||||
sup=0.d0
|
||||
inf=0.d0
|
||||
res=0.d0
|
||||
do i=2,nx-1
|
||||
*/
|
||||
f = vec![0.0; ny];
|
||||
diag = vec![0.0; ny];
|
||||
sup = vec![0.0; ny];
|
||||
inf = vec![0.0; ny];
|
||||
res = vec![0.0; ny];
|
||||
for i in 1..nx-1 {
|
||||
/*
|
||||
do j=2,ny-1
|
||||
factxp=(kdint(i+1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2
|
||||
factxm=(kdint(i-1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2
|
||||
factyp=(kdint(i,j+1)+kdint(i,j))/2.d0*(dt/2.)/dy**2
|
||||
factym=(kdint(i,j-1)+kdint(i,j))/2.d0*(dt/2.)/dy**2
|
||||
diag(j)=1.d0+factyp+factym
|
||||
sup(j)=-factyp
|
||||
inf(j)=-factym
|
||||
f(j)=zint(i,j)+factxp*zint(i+1,j)-(factxp+factxm)*zint(i,j)+factxm*zint(i-1,j)
|
||||
enddo
|
||||
*/
|
||||
for j in 1..ny-1 {
|
||||
factxp = (kdint[aij(i+1, j)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dx * dx);
|
||||
factxm = (kdint[aij(i-1, j)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dx * dx);
|
||||
factyp = (kdint[aij(i, j+1)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dy * dy);
|
||||
factym = (kdint[aij(i, j-1)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dy * dy);
|
||||
diag[j] = 1.0 + factyp + factym;
|
||||
sup[j] = -factyp;
|
||||
inf[j] = -factym;
|
||||
f[j] =
|
||||
zint[aij(i, j)] + factxp * zint[aij(i+1, j)] -
|
||||
(factxp + factxm) *
|
||||
zint[aij(i, j)] + factxm * zint[aij(i-1, j)];
|
||||
}
|
||||
/*
|
||||
! bottom bc
|
||||
if (cbc(1:1).eq.'1') then
|
||||
diag(1)=1.
|
||||
sup(1)=0.
|
||||
f(1)=zint(i,1)
|
||||
else
|
||||
factxp=(kdint(i+1,1)+kdint(i,j))/2.d0*(dt/2.)/dx**2
|
||||
factxm=(kdint(i-1,1)+kdint(i,1))/2.d0*(dt/2.)/dx**2
|
||||
factyp=(kdint(i,2)+kdint(i,1))/2.d0*(dt/2.)/dy**2
|
||||
diag(1)=1.d0+factyp
|
||||
sup(1)=-factyp
|
||||
f(1)=zint(i,1)+factxp*zint(i+1,1)-(factxp+factxm)*zint(i,1)+factxm*zint(i-1,1)
|
||||
endif
|
||||
*/
|
||||
if true {
|
||||
diag[0] = 1.0;
|
||||
sup[0] = 0.0;
|
||||
f[0] = zint[aij(i, 0)];
|
||||
} else {
|
||||
// FIXME: implement reflective boundary
|
||||
}
|
||||
/*
|
||||
! top bc
|
||||
if (cbc(3:3).eq.'1') then
|
||||
diag(ny)=1.
|
||||
inf(ny)=0.
|
||||
f(ny)=zint(i,ny)
|
||||
else
|
||||
factxp=(kdint(i+1,ny)+kdint(i,ny))/2.d0*(dt/2.)/dx**2
|
||||
factxm=(kdint(i-1,ny)+kdint(i,ny))/2.d0*(dt/2.)/dx**2
|
||||
factym=(kdint(i,ny-1)+kdint(i,ny))/2.d0*(dt/2.)/dy**2
|
||||
diag(ny)=1.d0+factym
|
||||
inf(ny)=-factym
|
||||
f(ny)=zint(i,ny)+factxp*zint(i+1,ny)-(factxp+factxm)*zint(i,ny)+factxm*zint(i-1,ny)
|
||||
endif
|
||||
*/
|
||||
if true {
|
||||
diag[ny - 1] = 1.0;
|
||||
inf[ny - 1] = 0.0;
|
||||
f[ny - 1] = zint[aij(i, ny - 1)];
|
||||
} else {
|
||||
// FIXME: implement reflective boundary
|
||||
}
|
||||
/*
|
||||
call tridag (inf,diag,sup,f,res,ny)
|
||||
do j=1,ny
|
||||
zintp(i,j)=res(j)
|
||||
enddo
|
||||
*/
|
||||
tridag(&inf, &diag, &sup, &f, &mut res, ny);
|
||||
for j in 0..ny {
|
||||
zintp[aij(i, j)] = res[j];
|
||||
}
|
||||
/*
|
||||
enddo
|
||||
deallocate (f,diag,sup,inf,res)
|
||||
*/
|
||||
}
|
||||
/*
|
||||
! stores result in 1D array
|
||||
|
||||
do j=1,ny
|
||||
do i=1,nx
|
||||
ij=(j-1)*nx+i
|
||||
etot(ij)=etot(ij)+h(ij)-zintp(i,j)
|
||||
erate(ij)=erate(ij)+(h(ij)-zintp(i,j))/dt
|
||||
h(ij)=zintp(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
b=min(h,b)
|
||||
*/
|
||||
for j in 0..ny {
|
||||
for i in 0..nx {
|
||||
ij = aij(i, j);
|
||||
// FIXME: Track total erosion and erosion rate.
|
||||
h[aij(i, j)] = zintp[aij(i, j)] as f32;
|
||||
}
|
||||
}
|
||||
|
||||
b.par_iter_mut().zip(h).for_each(|(mut b, h)| {
|
||||
*b = h.min(*b);
|
||||
});
|
||||
/*
|
||||
deallocate (zint,kdint,zintp)
|
||||
|
||||
return
|
||||
|
||||
end subroutine Diffusion
|
||||
*/
|
||||
}
|
||||
|
||||
/*
|
||||
!----------
|
||||
|
||||
! subroutine to solve a tri-diagonal system of equations (from Numerical Recipes)
|
||||
|
||||
SUBROUTINE tridag(a,b,c,r,u,n)
|
||||
|
||||
implicit none
|
||||
|
||||
INTEGER n
|
||||
double precision a(n),b(n),c(n),r(n),u(n)
|
||||
*/
|
||||
pub fn tridag(a: &[f64], b: &[f64], c: &[f64], r: &[f64], u: &mut [f64], n: usize) {
|
||||
/*
|
||||
INTEGER j
|
||||
double precision bet
|
||||
double precision,dimension(:),allocatable::gam
|
||||
|
||||
allocate (gam(n))
|
||||
|
||||
if(b(1).eq.0.d0) stop 'in tridag'
|
||||
*/
|
||||
let mut j : usize;
|
||||
let mut bet : f64;
|
||||
let mut precision : f64;
|
||||
let mut gam: Vec<f64>;
|
||||
|
||||
gam = vec![Default::default(); n];
|
||||
|
||||
assert!(b[0] != 0.0);
|
||||
/*
|
||||
|
||||
! first pass
|
||||
|
||||
bet=b(1)
|
||||
u(1)=r(1)/bet
|
||||
do 11 j=2,n
|
||||
gam(j)=c(j-1)/bet
|
||||
bet=b(j)-a(j)*gam(j)
|
||||
if(bet.eq.0.) then
|
||||
print*,'tridag failed'
|
||||
stop
|
||||
endif
|
||||
u(j)=(r(j)-a(j)*u(j-1))/bet
|
||||
11 continue
|
||||
*/
|
||||
bet = b[0];
|
||||
u[0] = r[0] / bet;
|
||||
for j in 1..n {
|
||||
gam[j] = c[j - 1] / bet;
|
||||
bet = b[j] - a[j] * gam[j];
|
||||
assert!(bet != 0.0);
|
||||
u[j] = (r[j] - a[j] * u[j - 1]) / bet;
|
||||
}
|
||||
/*
|
||||
! second pass
|
||||
|
||||
do 12 j=n-1,1,-1
|
||||
u(j)=u(j)-gam(j+1)*u(j+1)
|
||||
12 continue
|
||||
*/
|
||||
for j in (0..n - 1).rev() {
|
||||
u[j] = u[j] - gam[j + 1] * u[j + 1];
|
||||
}
|
||||
/*
|
||||
deallocate (gam)
|
||||
|
||||
return
|
||||
|
||||
END
|
||||
*/
|
||||
}
|
||||
@ -1,4 +1,5 @@
|
||||
use super::{downhill, neighbors, uniform_idx_as_vec2, uphill, WORLD_SIZE};
|
||||
use super::{diffusion, downhill, neighbors, uniform_idx_as_vec2, uphill, WORLD_SIZE};
|
||||
use bitvec::prelude::{bitbox, bitvec, BitBox};
|
||||
use crate::{config::CONFIG, util::RandomField};
|
||||
use common::{terrain::TerrainChunkSize, vol::RectVolSize};
|
||||
use noise::{NoiseFn, Point3};
|
||||
@ -157,6 +158,7 @@ pub fn get_rivers(
|
||||
let neighbor_coef = TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
|
||||
// NOTE: This technically makes us discontinuous, so we should be cautious about using this.
|
||||
let derivative_divisor = 1.0;
|
||||
let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale as f64;
|
||||
for &chunk_idx in newh.into_iter().rev() {
|
||||
let chunk_idx = chunk_idx as usize;
|
||||
let downhill_idx = downhill[chunk_idx];
|
||||
@ -351,7 +353,7 @@ pub fn get_rivers(
|
||||
// network).
|
||||
let downhill_water_alt = water_alt[downhill_idx];
|
||||
let neighbor_distance = neighbor_dim.magnitude();
|
||||
let dz = (downhill_water_alt - chunk_water_alt) * CONFIG.mountain_scale;
|
||||
let dz = (downhill_water_alt - chunk_water_alt) / height_scale as f32;// * CONFIG.mountain_scale;
|
||||
let slope = dz.abs() as f64 / neighbor_distance;
|
||||
if slope == 0.0 {
|
||||
// This is not a river--how did this even happen?
|
||||
@ -473,23 +475,28 @@ pub fn get_rivers(
|
||||
/// Precompute the maximum slope at all points.
|
||||
///
|
||||
/// TODO: See if allocating in advance is worthwhile.
|
||||
fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync)) -> Box<[f32]> {
|
||||
const MIN_MAX_ANGLE: f64 = 6.0 / 360.0 * 2.0 * f64::consts::PI;
|
||||
const MAX_MAX_ANGLE: f64 = 54.0 / 360.0 * 2.0 * f64::consts::PI;
|
||||
fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync)) -> Box<[f64]> {
|
||||
const MIN_MAX_ANGLE: f64 = 6.0/*30.0*//*6.0*//*15.0*/ / 360.0 * 2.0 * f64::consts::PI;
|
||||
const MAX_MAX_ANGLE: f64 = 54.0/*50.0*//*54.0*//*45.0*/ / 360.0 * 2.0 * f64::consts::PI;
|
||||
const MAX_ANGLE_RANGE: f64 = MAX_MAX_ANGLE - MIN_MAX_ANGLE;
|
||||
let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale as f64;
|
||||
h.par_iter()
|
||||
.enumerate()
|
||||
.map(|(posi, &z)| {
|
||||
let wposf = (uniform_idx_as_vec2(posi)).map(|e| e as f64);
|
||||
let wposz = z as f64 * CONFIG.mountain_scale as f64;
|
||||
let wposf = uniform_idx_as_vec2(posi).map(|e| e as f64) * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
|
||||
let wposz = z as f64 / height_scale;// * CONFIG.mountain_scale as f64;
|
||||
// Normalized to be between 6 and and 54 degrees.
|
||||
let div_factor = 32.0;
|
||||
let div_factor = /*32.0*//*16.0*//*64.0*/256.0/*1.0*//*4.0*/ * /*1.0*/4.0/* TerrainChunkSize::RECT_SIZE.x as f64 / 8.0 */;
|
||||
let rock_strength = rock_strength_nz
|
||||
.get([
|
||||
wposf.x, /* / div_factor*/
|
||||
wposf.y, /* / div_factor*/
|
||||
wposz,
|
||||
])
|
||||
wposz * div_factor,
|
||||
]);
|
||||
/* if rock_strength < -1.0 || rock_strength > 1.0 {
|
||||
println!("Nooooo: {:?}", rock_strength);
|
||||
} */
|
||||
let rock_strength = rock_strength
|
||||
.max(-1.0)
|
||||
.min(1.0)
|
||||
* 0.5
|
||||
@ -508,20 +515,22 @@ fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync
|
||||
// We do log-odds against center, so that our log odds are 0 when x = 0.25, lower when x is
|
||||
// lower, and higher when x is higher.
|
||||
//
|
||||
// (NOTE: below sea level, we invert it).
|
||||
//
|
||||
// TODO: Make all this stuff configurable... but honestly, it's so complicated that I'm not
|
||||
// sure anyone would be able to usefully tweak them on a per-map basis? Plus it's just a
|
||||
// hacky heuristic anyway.
|
||||
let center = 0.25;
|
||||
let delta = 0.05;
|
||||
let dmin = center - delta;
|
||||
let dmax = center + delta;
|
||||
let dmin = center - 0.15;//0.05;
|
||||
let dmax = center + 0.05;//0.05;
|
||||
let log_odds = |x: f64| logit(x) - logit(center);
|
||||
let rock_strength = logistic_cdf(
|
||||
1.0 * logit(rock_strength.min(1.0f64 - 1e-7).max(1e-7))
|
||||
+ 1.0 * log_odds((z as f64).min(dmax).max(dmin)),
|
||||
+ 1.0 * log_odds((wposz / CONFIG.mountain_scale as f64).abs().min(dmax).max(dmin)),
|
||||
);
|
||||
// let rock_strength = 0.5;
|
||||
let max_slope = (rock_strength * MAX_ANGLE_RANGE + MIN_MAX_ANGLE).tan();
|
||||
max_slope as f32
|
||||
max_slope
|
||||
})
|
||||
.collect::<Vec<_>>()
|
||||
.into_boxed_slice()
|
||||
@ -569,27 +578,73 @@ fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync
|
||||
/// always be the *new* h[j](t + δt), while h[i] will still not have been computed yet, so we
|
||||
/// only need to visit each node once.
|
||||
///
|
||||
/// Afterwards, we also apply a hillslope diffusion process using an ADI (alternating direction
|
||||
/// implicit) method:
|
||||
///
|
||||
/// https://github.com/fastscape-lem/fastscapelib-fortran/blob/master/src/Diffusion.f90
|
||||
///
|
||||
/// [1] Guillaume Cordonnier, Jean Braun, Marie-Paule Cani, Bedrich Benes, Eric Galin, et al..
|
||||
/// Large Scale Terrain Generation from Tectonic Uplift and Fluvial Erosion.
|
||||
/// Computer Graphics Forum, Wiley, 2016, Proc. EUROGRAPHICS 2016, 35 (2), pp.165-175.
|
||||
/// ⟨10.1111/cgf.12820⟩. ⟨hal-01262376⟩
|
||||
///
|
||||
///
|
||||
fn erode(
|
||||
h: &mut [f32],
|
||||
b: &mut [f32],
|
||||
is_done: &mut BitBox,
|
||||
done_val: bool,
|
||||
erosion_base: f32,
|
||||
max_uplift: f32,
|
||||
kdsed: f64,
|
||||
_seed: &RandomField,
|
||||
rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync),
|
||||
uplift: impl Fn(usize) -> f32,
|
||||
kd: impl Fn(usize) -> f64,
|
||||
is_ocean: impl Fn(usize) -> bool + Sync,
|
||||
) {
|
||||
log::debug!("Done draining...");
|
||||
let mmaxh = 1.0;
|
||||
// Landslide constant: ideally scaled to 10 m^-2 / y^-1
|
||||
let l = 200.0 * max_uplift as f64;
|
||||
// Erosion constant.
|
||||
let k = erosion_base as f64 + 2.244 / mmaxh as f64 * max_uplift as f64;
|
||||
let ((dh, indirection, newh, area), max_slope) = rayon::join(
|
||||
let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale as f64;
|
||||
let mmaxh = CONFIG.mountain_scale as f64 * height_scale;
|
||||
// Since maximum uplift rate is expected to be 5.010e-4 m * y^-1, and
|
||||
// 1.0 height units is 1.0 / height_scale m, whatever the
|
||||
// max uplift rate is (in units / y), we can find dt by multiplying by
|
||||
// 1.0 / height_scale m / unit and then dividing by 5.010e-4 m / y
|
||||
// (to get dt in y / unit). More formally:
|
||||
//
|
||||
// max_uplift h_unit / dt y = 5.010e-4 m / y
|
||||
//
|
||||
// 1 h_unit = 1.0 / height_scale m
|
||||
//
|
||||
// max_uplift h_unit / dt * 1.0 / height_scale m / h_unit =
|
||||
// max_uplift / height_scale m / dt =
|
||||
// 5.010e-4 m / y
|
||||
//
|
||||
// max_uplift / height_scale m / dt / (5.010e-4 m / y) = 1
|
||||
// (max_uplift / height_scale / 5.010e-4) y = dt
|
||||
let dt = max_uplift as f64 / height_scale /* * CONFIG.mountain_scale as f64*/ / 5.010e-4;
|
||||
println!("dt={:?}", dt);
|
||||
// Landslide constant: ideally scaled to 10e-2 m / y^-1
|
||||
let l = /*200.0 * max_uplift as f64;*/10.0e-2 /*/ CONFIG.mountain_scale as f64*/ * height_scale * dt;
|
||||
// Net precipitation rate (m / year)
|
||||
let p = 1.0 * height_scale;
|
||||
// Stream power erosion constant (bedrock), in m^(1-2m) / year (times dt).
|
||||
let k_fb = // erosion_base as f64 + 2.244 / mmaxh as f64 * max_uplift as f64;
|
||||
// 2.244*(5.010e-4)/512*5- (1.097e-5)
|
||||
// 2.244*(5.010e-4)/2048*5- (1.097e-5)
|
||||
// 2.244*(5.010e-4)/512- (8e-6)
|
||||
// 2.244*(5.010e-4)/512- (2e-6)
|
||||
// 2e-6 * dt;
|
||||
// 8e-6 * dt
|
||||
// 2e-5 * dt;
|
||||
// 2.244/2048*5*32/(250000/4)*10^6
|
||||
erosion_base as f64 + 2.244 / mmaxh as f64 * 5.0 * max_uplift as f64;
|
||||
// 1e-5 * dt;
|
||||
// (2.244*(5.010e-4)/512)/(2.244*(5.010e-4)/2500) = 4.88...
|
||||
// 2.444 * 5
|
||||
// Stream power erosion constant (sediment), in m^(1-2m) / year (times dt).
|
||||
let k_fs = k_fb * 1.0/* 2.0*//*4.0*/;
|
||||
let ((dh, indirection, newh, area), mut max_slopes) = rayon::join(
|
||||
|| {
|
||||
let mut dh = downhill(h, |posi| is_ocean(posi));
|
||||
log::debug!("Computed downhill...");
|
||||
@ -622,8 +677,10 @@ fn erode(
|
||||
if posj == -1 {
|
||||
panic!("Disconnected lake!");
|
||||
}
|
||||
// max_slopes[posi] = kd(posi);
|
||||
// Egress with no outgoing flows.
|
||||
} else {
|
||||
// *is_done.at(posi) = done_val;
|
||||
nland += 1;
|
||||
let posj = posj as usize;
|
||||
let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posj)).map(|e| e as f64);
|
||||
@ -633,23 +690,53 @@ fn erode(
|
||||
let neighbor_distance = (neighbor_coef * dxy).magnitude();
|
||||
// Since the area is in meters^(2m) and neighbor_distance is in m, so long as m = 0.5,
|
||||
// we have meters^(1) / meters^(1), so they should cancel out. Otherwise, we would
|
||||
// want to divide neighbor_distance by CONFIG.mountain_scale and area[posi] by
|
||||
// CONFIG.mountain_scale^2, to make sure we were working in the correct units for dz
|
||||
// (which has 1.0 height unit = CONFIG.mountain_scale meters).
|
||||
let flux = k * (chunk_area * area[posi] as f64).sqrt() / neighbor_distance;
|
||||
// want to multiply neighbor_distance by height_scale and area[posi] by
|
||||
// height_scale^2, to make sure we were working in the correct units for dz
|
||||
// (which has height_scale height unit = 1.0 meters).
|
||||
let old_h_i = h[posi] as f64;
|
||||
let uplift_i = uplift(posi) as f64;
|
||||
let k = if (old_h_i - b[posi] as f64) > 1.0e-6 {
|
||||
// Sediment
|
||||
k_fs
|
||||
} else {
|
||||
// Bedrock
|
||||
k_fb
|
||||
};
|
||||
// let k = k * uplift_i / max_uplift as f64;
|
||||
let flux = k * (p * chunk_area * area[posi] as f64).sqrt() / neighbor_distance;
|
||||
// h[i](t + dt) = (h[i](t) + δt * (uplift[i] + flux(i) * h[j](t + δt))) / (1 + flux(i) * δt).
|
||||
// NOTE: posj has already been computed since it's downhill from us.
|
||||
let h_j = h[posj] as f64;
|
||||
let old_h_i = h[posi] as f64;
|
||||
let mut new_h_i = (old_h_i + (uplift(posi) as f64 + flux * h_j)) / (1.0 + flux);
|
||||
// hi(t + ∂t) = (hi(t) + ∂t(ui + kp^mAi^m(hj(t + ∂t)/||pi - pj||))) / (1 + ∂t * kp^mAi^m / ||pi - pj||)
|
||||
let mut new_h_i = (old_h_i + (uplift_i + flux * h_j)) / (1.0 + flux);
|
||||
// Prevent erosion from dropping us below our receiver, unless we were already below it.
|
||||
// new_h_i = h_j.min(old_h_i + uplift_i).max(new_h_i);
|
||||
// Find out if this is a lake bottom.
|
||||
let indirection_idx = indirection[posi];
|
||||
let is_lake_bottom = indirection_idx < 0;
|
||||
// Test the slope.
|
||||
let max_slope = max_slope[posi] as f64;
|
||||
let dz = (new_h_i - h_j).max(0.0) * CONFIG.mountain_scale as f64;
|
||||
let mag_slope = dz/*.abs()*/ / neighbor_distance;
|
||||
let _fake_neighbor = is_lake_bottom && dxy.x.abs() > 1.0 && dxy.y.abs() > 1.0;
|
||||
// Test the slope.
|
||||
let max_slope = max_slopes[posi] as f64;
|
||||
// Hacky version of thermal erosion: only consider lowest neighbor, don't redistribute
|
||||
// uplift to other neighbors.
|
||||
let (posk, h_k) = /* neighbors(posi)
|
||||
.filter(|&posk| *is_done.at(posk) == done_val)
|
||||
// .filter(|&posk| *is_done.at(posk) == done_val || is_ocean(posk))
|
||||
.map(|posk| (posk, h[posk] as f64))
|
||||
// .filter(|&(posk, h_k)| *is_done.at(posk) == done_val || h_k < 0.0)
|
||||
.min_by(|&(_, a), &(_, b)| a.partial_cmp(&b).unwrap())
|
||||
.unwrap_or((posj, h_j)); */
|
||||
(posj, h_j);
|
||||
// .max(h_j);
|
||||
let (posk, h_k) = if h_k < h_j {
|
||||
(posk, h_k)
|
||||
} else {
|
||||
(posj, h_j)
|
||||
};
|
||||
let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posk)).map(|e| e as f64);
|
||||
let neighbor_distance = (neighbor_coef * dxy).magnitude();
|
||||
let dz = (new_h_i - /*h_j*/h_k).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
|
||||
let mag_slope = dz/*.abs()*/ / neighbor_distance;
|
||||
// If you're on the lake bottom and not right next to your neighbor, don't compute a
|
||||
// slope.
|
||||
if
|
||||
@ -661,26 +748,57 @@ fn erode(
|
||||
// println!("old slope: {:?}, new slope: {:?}, dz: {:?}, h_j: {:?}, new_h_i: {:?}", mag_slope, max_slope, dz, h_j, new_h_i);
|
||||
// Thermal erosion says this can't happen, so we reduce dh_i to make the slope
|
||||
// exactly max_slope.
|
||||
// max_slope = (old_h_i + dh - h_j) * CONFIG.mountain_scale / NEIGHBOR_DISTANCE
|
||||
// dh = max_slope * NEIGHBOR_DISTANCE / CONFIG.mountain_scale + h_j - old_h_i.
|
||||
let dh = max_slope * neighbor_distance / CONFIG.mountain_scale as f64;
|
||||
new_h_i = (h_j + dh).max(new_h_i - l * (mag_slope - max_slope));
|
||||
let dz = (new_h_i - h_j).max(0.0) * CONFIG.mountain_scale as f64;
|
||||
// max_slope = (old_h_i + dh - h_j) / height_scale/* * CONFIG.mountain_scale */ / NEIGHBOR_DISTANCE
|
||||
// dh = max_slope * NEIGHBOR_DISTANCE * height_scale/* / CONFIG.mountain_scale */ + h_j - old_h_i.
|
||||
let dh = max_slope * neighbor_distance * height_scale/* / CONFIG.mountain_scale as f64*/;
|
||||
new_h_i = /*h_j.max*/(h_k + dh).max(new_h_i - l * (mag_slope - max_slope));
|
||||
let dz = (new_h_i - /*h_j*/h_k).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
|
||||
let slope = dz/*.abs()*/ / neighbor_distance;
|
||||
sums += slope;
|
||||
/* max_slopes[posi] = /*(mag_slope - max_slope) * */kd(posi);
|
||||
sums += mag_slope; */
|
||||
// let slope = dz.signum() * max_slope;
|
||||
// new_h_i = slope * neighbor_distance / CONFIG.mountain_scale as f64 + h_j;
|
||||
// new_h_i = slope * neighbor_distance * height_scale /* / CONFIG.mountain_scale as f64 */ + h_j;
|
||||
// sums += max_slope;
|
||||
} else {
|
||||
// max_slopes[posi] = 0.0;
|
||||
sums += mag_slope;
|
||||
// Just use the computed rate.
|
||||
}
|
||||
h[posi] = new_h_i as f32;
|
||||
// Make sure to update the basement as well!
|
||||
b[posi] = (old_h_i + uplift_i).min(new_h_i) as f32;
|
||||
sumh += new_h_i;
|
||||
}
|
||||
}
|
||||
*is_done.at(posi) = done_val;
|
||||
maxh = h[posi].max(maxh);
|
||||
}
|
||||
log::debug!(
|
||||
"Done applying stream power (max height: {:?}) (avg height: {:?}) (avg slope: {:?})",
|
||||
maxh,
|
||||
if nland == 0 {
|
||||
f64::INFINITY
|
||||
} else {
|
||||
sumh / nland as f64
|
||||
},
|
||||
if nland == 0 {
|
||||
f64::INFINITY
|
||||
} else {
|
||||
sums / nland as f64
|
||||
},
|
||||
);
|
||||
|
||||
// Apply hillslope diffusion.
|
||||
diffusion(WORLD_SIZE.x, WORLD_SIZE.y,
|
||||
WORLD_SIZE.x as f64 * TerrainChunkSize::RECT_SIZE.x as f64 * height_scale/* / CONFIG.mountain_scale as f64*/,
|
||||
WORLD_SIZE.y as f64 * TerrainChunkSize::RECT_SIZE.y as f64 * height_scale/* / CONFIG.mountain_scale as f64*/,
|
||||
dt,
|
||||
(),
|
||||
h, b,
|
||||
/*|posi| max_slopes[posi]*/kd,
|
||||
kdsed,
|
||||
);
|
||||
log::debug!(
|
||||
"Done eroding (max height: {:?}) (avg height: {:?}) (avg slope: {:?})",
|
||||
maxh,
|
||||
@ -708,7 +826,7 @@ pub fn fill_sinks(
|
||||
) -> Box<[f32]> {
|
||||
// NOTE: We are using the "exact" version of depression-filling, which is slower but doesn't
|
||||
// change altitudes.
|
||||
let epsilon = /*1.0 / (1 << 7) as f32 / CONFIG.mountain_scale*/0.0;
|
||||
let epsilon = /*1.0 / (1 << 7) as f32 * height_scale/* / CONFIG.mountain_scale */*/0.0;
|
||||
let infinity = f32::INFINITY;
|
||||
let range = 0..WORLD_SIZE.x * WORLD_SIZE.y;
|
||||
let oldh = range
|
||||
@ -1122,14 +1240,21 @@ pub fn do_erosion(
|
||||
seed: &RandomField,
|
||||
rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync),
|
||||
oldh: impl Fn(usize) -> f32 + Sync,
|
||||
oldb: impl Fn(usize) -> f32 + Sync,
|
||||
is_ocean: impl Fn(usize) -> bool + Sync,
|
||||
uplift: impl Fn(usize) -> f32 + Sync,
|
||||
) -> Box<[f32]> {
|
||||
) -> (Box<[f32]>, Box<[f32]>) {
|
||||
let oldh_ = (0..WORLD_SIZE.x * WORLD_SIZE.y)
|
||||
.into_par_iter()
|
||||
.map(|posi| oldh(posi))
|
||||
.collect::<Vec<_>>()
|
||||
.into_boxed_slice();
|
||||
// Topographic basement (The height of bedrock, not including sediment).
|
||||
let mut b = (0..WORLD_SIZE.x * WORLD_SIZE.y)
|
||||
.into_par_iter()
|
||||
.map(|posi| oldb(posi))
|
||||
.collect::<Vec<_>>()
|
||||
.into_boxed_slice();
|
||||
// TODO: Don't do this, maybe?
|
||||
// (To elaborate, maybe we should have varying uplift or compute it some other way).
|
||||
let uplift = (0..oldh_.len())
|
||||
@ -1150,18 +1275,37 @@ pub fn do_erosion(
|
||||
.max_by(|a, b| a.partial_cmp(&b).unwrap())
|
||||
.unwrap();
|
||||
log::debug!("Max uplift: {:?}", max_uplift);
|
||||
// Height of terrain, including sediment.
|
||||
let mut h = oldh_;
|
||||
// 0.01 / 2e-5 = 500
|
||||
// Bedrock transport coefficients (diffusivity) in m^2 / year. For now, we set them all to be equal
|
||||
// on land, but in theory we probably want to at least differentiate between soil, bedrock, and
|
||||
// sediment.
|
||||
let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale as f64;
|
||||
let kdsed = /*1.5e-2*//*1e-4*/1e-2 * height_scale * height_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */;
|
||||
let kd_bedrock = /*1e-2*//*0.25e-2*/1e-2 * height_scale * height_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */;
|
||||
let kd = |posi: usize| if is_ocean(posi) { /*0.0*/kd_bedrock } else { kd_bedrock };
|
||||
// Hillslope diffusion coefficient for sediment.
|
||||
let mut is_done = bitbox![0; WORLD_SIZE.x * WORLD_SIZE.y];
|
||||
for i in 0..n {
|
||||
log::debug!("Erosion iteration #{:?}", i);
|
||||
erode(
|
||||
&mut h,
|
||||
&mut b,
|
||||
&mut is_done,
|
||||
// The value to use to indicate that erosion is complete on a chunk. Should toggle
|
||||
// once per iteration, to avoid having to reset the bits, and start at true, since
|
||||
// we initialize to 0 (false).
|
||||
i & 1 == 0,
|
||||
erosion_base,
|
||||
max_uplift,
|
||||
kdsed,
|
||||
seed,
|
||||
rock_strength_nz,
|
||||
|posi| uplift[posi],
|
||||
kd,
|
||||
|posi| is_ocean(posi),
|
||||
);
|
||||
}
|
||||
h
|
||||
(h, b)
|
||||
}
|
||||
|
||||
@ -1,16 +1,19 @@
|
||||
mod diffusion;
|
||||
mod erosion;
|
||||
mod location;
|
||||
mod settlement;
|
||||
mod util;
|
||||
|
||||
// Reexports
|
||||
pub use self::diffusion::diffusion;
|
||||
pub use self::erosion::{
|
||||
do_erosion, fill_sinks, get_drainage, get_lakes, get_rivers, RiverData, RiverKind,
|
||||
};
|
||||
pub use self::location::Location;
|
||||
pub use self::settlement::Settlement;
|
||||
pub use self::util::{
|
||||
cdf_irwin_hall, downhill, get_oceans, local_cells, map_edge_factor, neighbors,
|
||||
cdf_irwin_hall, downhill, get_oceans, HybridMulti as HybridMulti_, local_cells, map_edge_factor, neighbors,
|
||||
ScaleBias,
|
||||
uniform_idx_as_vec2, uniform_noise, uphill, vec2_as_uniform_idx, InverseCdf,
|
||||
};
|
||||
|
||||
@ -60,6 +63,7 @@ struct GenCdf {
|
||||
temp_base: InverseCdf,
|
||||
chaos: InverseCdf,
|
||||
alt: Box<[f32]>,
|
||||
basement: Box<[f32]>,
|
||||
water_alt: Box<[f32]>,
|
||||
dh: Box<[isize]>,
|
||||
/// NOTE: Until we hit 4096 × 4096, this should suffice since integers with an absolute value
|
||||
@ -74,7 +78,7 @@ pub(crate) struct GenCtx {
|
||||
pub turb_x_nz: SuperSimplex,
|
||||
pub turb_y_nz: SuperSimplex,
|
||||
pub chaos_nz: RidgedMulti,
|
||||
pub alt_nz: HybridMulti,
|
||||
pub alt_nz: HybridMulti_,
|
||||
pub hill_nz: SuperSimplex,
|
||||
pub temp_nz: Fbm,
|
||||
// Humidity noise
|
||||
@ -111,28 +115,34 @@ pub struct WorldSim {
|
||||
impl WorldSim {
|
||||
pub fn generate(seed: u32) -> Self {
|
||||
let mut rng = ChaChaRng::from_seed(seed_expan::rng_state(seed));
|
||||
let continent_scale = 5_000.0/*32768.0*/;
|
||||
let continent_scale = 5_000.0f64/*32768.0*/.div(32.0).mul(TerrainChunkSize::RECT_SIZE.x as f64);
|
||||
|
||||
let gen_ctx = GenCtx {
|
||||
turb_x_nz: SuperSimplex::new().set_seed(rng.gen()),
|
||||
turb_y_nz: SuperSimplex::new().set_seed(rng.gen()),
|
||||
chaos_nz: RidgedMulti::new()
|
||||
.set_octaves(/*7*//*3*/ 7)
|
||||
.set_octaves(/*7*//*3*/ /*7*//*3*/7)
|
||||
.set_frequency(
|
||||
/*RidgedMulti::DEFAULT_FREQUENCY **/ 3_000.0 * 8.0 / continent_scale,
|
||||
)
|
||||
// .set_persistence(RidgedMulti::DEFAULT_LACUNARITY.powf(-(1.0 - 0.5)))
|
||||
.set_seed(rng.gen()),
|
||||
hill_nz: SuperSimplex::new().set_seed(rng.gen()),
|
||||
alt_nz: HybridMulti::new()
|
||||
.set_octaves(/*3*//*2*/ 8)
|
||||
alt_nz: HybridMulti_::new()
|
||||
.set_octaves(/*3*//*2*/ /*8*//*3*/8)
|
||||
// 1/2048*32*1024 = 16
|
||||
.set_frequency(
|
||||
/*HybridMulti::DEFAULT_FREQUENCY*/
|
||||
// (2^8*(10000/5000/10000))*32 = per-chunk
|
||||
(10_000.0/* * 2.0*/ / continent_scale) as f64,
|
||||
)
|
||||
// .set_frequency(1.0 / ((1 << 0) as f64))
|
||||
// .set_lacunarity(1.0)
|
||||
.set_persistence(/*0.5*//*0.5*/ 0.5)
|
||||
// persistence = lacunarity^(-(1.0 - fractal increment))
|
||||
.set_persistence(HybridMulti_::DEFAULT_LACUNARITY.powf(-(1.0 - /*0.75*/0.75)))
|
||||
// .set_persistence(/*0.5*//*0.5*/0.5 + 1.0 / ((1 << 6) as f64))
|
||||
// .set_offset(/*0.7*//*0.5*//*0.75*/0.7)
|
||||
// .set_offset(/*0.7*//*0.5*//*0.75*/0.0)
|
||||
.set_seed(rng.gen()),
|
||||
//temp_nz: SuperSimplex::new().set_seed(rng.gen()),
|
||||
temp_nz: Fbm::new()
|
||||
@ -178,17 +188,81 @@ impl WorldSim {
|
||||
};
|
||||
|
||||
let river_seed = RandomField::new(rng.gen());
|
||||
let rock_strength_nz = Fbm::new()
|
||||
.set_octaves(8)
|
||||
.set_persistence(/*0.9*/ 2.0)
|
||||
.set_frequency(/*0.9*/ Fbm::DEFAULT_FREQUENCY / (64.0 * 32.0))
|
||||
let rock_lacunarity = 0.5/*HybridMulti::DEFAULT_LACUNARITY*/;
|
||||
let rock_strength_nz = /*Fbm*/HybridMulti_/*BasicMulti*//*Fbm*/::new()
|
||||
.set_octaves(/*6*//*5*//*4*//*5*/4)
|
||||
// persistence = lacunarity^(-(1.0 - fractal increment))
|
||||
.set_persistence(/*0.9*/ /*2.0*//*1.5*//*HybridMulti::DEFAULT_LACUNARITY*/rock_lacunarity.powf(-(1.0 - 0.9)))
|
||||
// 256*32/2^4
|
||||
// (0.5^(-(1.0-0.9)))^4/256/32*2^4*16*32
|
||||
// (0.5^(-(1.0-0.9)))^4/256/32*2^4*256*4
|
||||
// (0.5^(-(1.0-0.9)))^1/256/32*2^4*256*4
|
||||
// (2^(-(1.0-0.9)))^4
|
||||
// 16.0
|
||||
.set_frequency(/*0.9*/ /*Fbm*/HybridMulti_::DEFAULT_FREQUENCY / (/*64.0*/256.0/*1.0*//*16.0*/ * 32.0/* TerrainChunkSize::RECT_SIZE.x as f64 */))
|
||||
// .set_persistence(/*0.9*/ /*2.0*/0.67)
|
||||
// .set_frequency(/*0.9*/ Fbm::DEFAULT_FREQUENCY / (2.0 * 32.0))
|
||||
// .set_lacunarity(0.5)
|
||||
.set_seed(rng.gen());
|
||||
// NOTE: octaves should definitely fit into i32, but we should check anyway to make
|
||||
// sure.
|
||||
/* assert!(rock_strength_nz.persistence > 0.0);
|
||||
let rock_strength_scale = (1..rock_strength_nz.octaves as i32)
|
||||
.map(|octave| rock_strength_nz.persistence.powi(octave + 1))
|
||||
.sum::<f64>()
|
||||
// For some reason, this is "scaled" by 3.0.
|
||||
.mul(3.0);
|
||||
let rock_strength_nz = ScaleBias::new(&rock_strength_nz)
|
||||
.set_scale(1.0 / rock_strength_scale); */
|
||||
|
||||
let max_erosion_per_delta_t = 32.0 / CONFIG.mountain_scale as f64;
|
||||
let height_scale = 1.0f64; // 1.0 / CONFIG.mountain_scale as f64;
|
||||
let max_erosion_per_delta_t = 32.0 * height_scale;
|
||||
let erosion_pow_low = /*0.25*//*1.5*//*2.0*//*0.5*//*4.0*//*0.25*//*1.0*//*2.0*//*1.5*//*1.5*//*0.35*//*0.43*//*0.5*//*0.45*//*0.37*/1.002;
|
||||
let erosion_pow_high = /*1.5*//*1.0*//*0.55*//*0.51*//*2.0*/1.002;
|
||||
let erosion_center = /*0.45*//*0.75*//*0.75*//*0.5*//*0.75*/0.5;
|
||||
let n_steps = 150; //150;//200;
|
||||
let n_steps = 150;//37/*100*/;//50;//50;//37;//50;//37; // /*37*//*29*//*40*//*150*/37; //150;//200;
|
||||
let n_small_steps = 0;//8;//8; // 8
|
||||
|
||||
// fractal dimension should be between 0 and 0.9999...
|
||||
// (lacunarity^octaves)^(-H) = persistence^(octaves)
|
||||
// lacunarity^(octaves*-H) = persistence^(octaves)
|
||||
// e^(-octaves*H*ln(lacunarity)) = e^(octaves * ln(persistence))
|
||||
// -octaves * H * ln(lacunarity) = octaves * ln(persistence)
|
||||
// -H = ln(persistence) / ln(lacunarity)
|
||||
// H = -ln(persistence) / ln(lacunarity)
|
||||
// ln(persistence) = -H * ln(lacunarity)
|
||||
// persistence = lacunarity^(-H)
|
||||
//
|
||||
// -ln(2^(-0.25))/ln(2) = 0.25
|
||||
//
|
||||
// -ln(2^(-0.1))/ln(2)
|
||||
//
|
||||
// 0 = -ln(persistence) / ln(lacunarity)
|
||||
// 0 = ln(persistence) => persistence = e^0 = 1
|
||||
//
|
||||
// 1 = -ln(persistence) / ln(lacunarity)
|
||||
// -ln(lacunarity) = ln(persistence)
|
||||
// e^(-ln(lacunarity)) = e^(ln(persistence))
|
||||
// 1 / lacunarity = persistence
|
||||
//
|
||||
// Ergo, we should not set fractal dimension to anything not between 1 / lacunarity and 1.
|
||||
//
|
||||
// dimension = -ln(0.25)/ln(2*pi/3) = 1.875
|
||||
//
|
||||
// (2*pi/3^1)^(-(-ln(0.25)/ln(2*pi/3))) = 0.25
|
||||
//
|
||||
// Default should be at most 1 / lacunarity.
|
||||
//
|
||||
// (2 * pi / 3)^(-ln(0.25)/ln(2*pi/3))
|
||||
//
|
||||
// -ln(0.25)/ln(2*pi/3) = 1.88
|
||||
//
|
||||
// (2 * pi / 3)^(-ln(0.25)/ln(2*pi/3))
|
||||
//
|
||||
// 2 * pi / 3
|
||||
//
|
||||
// 2.0^(2(-ln(1.5)/ln(2)))
|
||||
// (1 / 1.5)^(2)
|
||||
|
||||
// No NaNs in these uniform vectors, since the original noise value always returns Some.
|
||||
let ((alt_base, _), (chaos, _)) = rayon::join(
|
||||
@ -202,8 +276,13 @@ impl WorldSim {
|
||||
.alt_nz
|
||||
.get((wposf.div(10_000.0)).into_array())
|
||||
.min(1.0)
|
||||
.max(-1.0))
|
||||
.sub(0.05)
|
||||
.max(-1.0)
|
||||
/* .mul(0.25)
|
||||
.add(0.125) */)
|
||||
// .add(0.5)
|
||||
// .sub(0.05)
|
||||
// .add(0.05)
|
||||
.add(0.075)
|
||||
.mul(0.35), /*-0.0175*/
|
||||
)
|
||||
})
|
||||
@ -217,13 +296,13 @@ impl WorldSim {
|
||||
//.add(0.0)
|
||||
+ gen_ctx
|
||||
.hill_nz
|
||||
.get((wposf.div(1_500.0)).into_array())
|
||||
.get((wposf.mul(32.0).div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(1_500.0)).into_array())
|
||||
.min(1.0)
|
||||
.max(-1.0)
|
||||
.mul(1.0)
|
||||
+ gen_ctx
|
||||
.hill_nz
|
||||
.get((wposf.div(400.0)).into_array())
|
||||
.get((wposf.mul(32.0).div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(400.0)).into_array())
|
||||
.min(1.0)
|
||||
.max(-1.0)
|
||||
.mul(0.3))
|
||||
@ -296,6 +375,7 @@ impl WorldSim {
|
||||
.min(1.0)
|
||||
.max(-1.0))
|
||||
.abs()
|
||||
// 0.5
|
||||
.powf(1.35);
|
||||
|
||||
fn spring(x: f64, pow: f64) -> f64 {
|
||||
@ -305,14 +385,14 @@ impl WorldSim {
|
||||
(0.0 + alt_main/*0.4*/
|
||||
+ (gen_ctx
|
||||
.small_nz
|
||||
.get((wposf.div(300.0)).into_array())
|
||||
.get((wposf.mul(32.0).div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(300.0)).into_array())
|
||||
.min(1.0)
|
||||
.max(-1.0))
|
||||
.mul(alt_main.powf(0.8).max(/*0.25*/ 0.15))
|
||||
.mul(0.3)
|
||||
.add(1.0)
|
||||
.mul(0.4)
|
||||
/*0.52*/
|
||||
// 0.52
|
||||
+ spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0).mul(0.045))
|
||||
};
|
||||
|
||||
@ -357,7 +437,11 @@ impl WorldSim {
|
||||
});
|
||||
|
||||
// Calculate oceans.
|
||||
let old_height = |posi: usize| alt_old[posi].1;
|
||||
let old_height = |posi: usize| alt_old[posi].1 * CONFIG.mountain_scale * height_scale as f32;
|
||||
/* let is_ocean = (0..WORLD_SIZE.x * WORLD_SIZE.y)
|
||||
.into_par_iter()
|
||||
.map(|i| map_edge_factor(i) == 0.0)
|
||||
.collect::<Vec<_>>(); */
|
||||
let is_ocean = get_oceans(old_height);
|
||||
let is_ocean_fn = |posi: usize| is_ocean[posi];
|
||||
|
||||
@ -406,6 +490,8 @@ impl WorldSim {
|
||||
let logistic_cdf = |x: f64| (x / logistic_2_base).tanh() * 0.5 + 0.5;
|
||||
|
||||
let exp_inverse_cdf = |x: f64/*, pow: f64*/| -(-x).ln_1p()/* / ln(pow)*/;
|
||||
// 2 / pi * ln(tan(pi/2 * p))
|
||||
let hypsec_inverse_cdf = |x: f64| f64::consts::FRAC_2_PI * ((x * f64::consts::FRAC_PI_2).tan().ln());
|
||||
// 2^((2^10-2)/256) = 15.91...
|
||||
// -ln(1-(1-(2^(-22)*0.5)))
|
||||
// -ln(1-(1-(2^(-53)*0.5)))
|
||||
@ -423,8 +509,9 @@ impl WorldSim {
|
||||
1.0 / (WORLD_SIZE.x as f64 * WORLD_SIZE.y as f64).max(f64::EPSILON as f64 * 0.5);
|
||||
let max_epsilon = (1.0 - 1.0 / (WORLD_SIZE.x as f64 * WORLD_SIZE.y as f64))
|
||||
.min(1.0 - f64::EPSILON as f64 * 0.5);
|
||||
let alt_exp_min_uniform = exp_inverse_cdf(min_epsilon);
|
||||
let alt_exp_max_uniform = exp_inverse_cdf(max_epsilon);
|
||||
let inv_func = /*|x: f64| x*//*exp_inverse_cdf*/logit/*hypsec_inverse_cdf*/;
|
||||
let alt_exp_min_uniform = /*exp_inverse_cdf*//*logit*/inv_func(min_epsilon);
|
||||
let alt_exp_max_uniform = /*exp_inverse_cdf*//*logit*/inv_func(max_epsilon);
|
||||
|
||||
// let erosion_pow = 2.0;
|
||||
// let n_steps = 100;//150;
|
||||
@ -437,58 +524,9 @@ impl WorldSim {
|
||||
};
|
||||
/* let erosion_factor = |x: f64| logistic_cdf(logistic_base * if x <= /*erosion_center*/alt_old_center_uniform/*alt_old_center*/ { erosion_pow_low.ln() } else { erosion_pow_high.ln() } * log_odds(x))/*0.5 + (x - 0.5).signum() * ((x - 0.5).mul(2.0).abs(
|
||||
).powf(erosion_pow).mul(0.5))*/; */
|
||||
let erosion_factor = |x: f64| (/*if x <= /*erosion_center*/alt_old_center_uniform/*alt_old_center*/ { erosion_pow_low.ln() } else { erosion_pow_high.ln() } * */(exp_inverse_cdf(x) - alt_exp_min_uniform) / (alt_exp_max_uniform - alt_exp_min_uniform))/*0.5 + (x - 0.5).signum() * ((x - 0.5).mul(2.0).abs(
|
||||
).powf(erosion_pow).mul(0.5))*/;
|
||||
let alt = do_erosion(
|
||||
0.0,
|
||||
max_erosion_per_delta_t as f32,
|
||||
n_steps,
|
||||
&river_seed,
|
||||
&rock_strength_nz,
|
||||
|posi| {
|
||||
if is_ocean_fn(posi) {
|
||||
old_height(posi)
|
||||
} else {
|
||||
let wposf = (uniform_idx_as_vec2(posi)
|
||||
* TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
|
||||
.map(|e| e as f64);
|
||||
let alt_main = {
|
||||
// Extension upwards from the base. A positive number from 0 to 1 curved to be
|
||||
// maximal at 0. Also to be multiplied by CONFIG.mountain_scale.
|
||||
let alt_main = (gen_ctx
|
||||
.alt_nz
|
||||
.get((wposf.div(2_000.0)).into_array())
|
||||
.min(1.0)
|
||||
.max(-1.0))
|
||||
.abs()
|
||||
.powf(1.35);
|
||||
|
||||
fn spring(x: f64, pow: f64) -> f64 {
|
||||
x.abs().powf(pow) * x.signum()
|
||||
}
|
||||
|
||||
(0.0 + alt_main
|
||||
+ (gen_ctx
|
||||
.small_nz
|
||||
.get((wposf.div(300.0)).into_array())
|
||||
.min(1.0)
|
||||
.max(-1.0))
|
||||
.mul(alt_main.powf(0.8).max(/*0.25*/ 0.15))
|
||||
.mul(0.3)
|
||||
.add(1.0)
|
||||
.mul(0.4)
|
||||
+ spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0)
|
||||
.mul(0.045))
|
||||
};
|
||||
// old_height_uniform(posi) *
|
||||
(/*((old_height(posi) - alt_old_min) as f64 / (alt_old_max - alt_old_min) as f64) **/(((6.0 / 360.0 * 2.0 * f64::consts::PI).tan()
|
||||
* TerrainChunkSize::RECT_SIZE.reduce_partial_min() as f64)
|
||||
.floor()
|
||||
/ CONFIG.mountain_scale as f64)) as f32
|
||||
// 5.0 / CONFIG.mountain_scale
|
||||
}
|
||||
},
|
||||
is_ocean_fn,
|
||||
let erosion_factor = |x: f64| (/*if x <= /*erosion_center*/alt_old_center_uniform/*alt_old_center*/ { erosion_pow_low.ln() } else { erosion_pow_high.ln() } * */(/*exp_inverse_cdf*//*logit*/inv_func(x) - alt_exp_min_uniform) / (alt_exp_max_uniform - alt_exp_min_uniform))/*0.5 + (x - 0.5).signum() * ((x - 0.5).mul(2.0).abs(
|
||||
).powf(erosion_pow).mul(0.5))*//*.powf(0.5)*//*.powf(1.5)*//*.powf(2.0)*/;
|
||||
let uplift_fn =
|
||||
|posi| {
|
||||
if is_ocean_fn(posi) {
|
||||
return 0.0;
|
||||
@ -521,8 +559,8 @@ impl WorldSim {
|
||||
.mul(0.3)
|
||||
.add(1.0)
|
||||
.mul(0.4)
|
||||
+ spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0)
|
||||
.mul(0.045))
|
||||
/* + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0)
|
||||
.mul(0.045)*/)
|
||||
};
|
||||
let height =
|
||||
((old_height_uniform(posi) - alt_old_min_uniform) as f64
|
||||
@ -577,16 +615,98 @@ impl WorldSim {
|
||||
} else {
|
||||
(0.0, 0.0)
|
||||
};
|
||||
let height = 1.0f64;
|
||||
// let height = 1.0 / 7.0f64;
|
||||
let height = height
|
||||
.mul(max_erosion_per_delta_t * 7.0 / 8.0)
|
||||
.add(max_erosion_per_delta_t / 8.0)
|
||||
/* .mul(15.0 / 16.0)
|
||||
.add(1.0 / 16.0) */
|
||||
.mul(7.0 / 8.0)
|
||||
.add(1.0 / 8.0)
|
||||
.mul(max_erosion_per_delta_t)
|
||||
.sub(/*1.0 / CONFIG.mountain_scale as f64*/ bump_max)
|
||||
.add(bump_factor);
|
||||
/* .sub(/*1.0 / CONFIG.mountain_scale as f64*/(f32::EPSILON * 0.5) as f64)
|
||||
.add(bump_factor); */
|
||||
height as f32
|
||||
},
|
||||
};
|
||||
let alt_func = |posi| {
|
||||
if is_ocean_fn(posi) {
|
||||
// -max_erosion_per_delta_t as f32
|
||||
// -1.0 / CONFIG.mountain_scale
|
||||
// -0.75
|
||||
// -CONFIG.sea_level / CONFIG.mountain_scale
|
||||
// 0.0
|
||||
old_height(posi) // 0.0
|
||||
} else {
|
||||
// uplift_fn(posi)
|
||||
let wposf = (uniform_idx_as_vec2(posi)
|
||||
* TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
|
||||
.map(|e| e as f64);
|
||||
let alt_main = {
|
||||
// Extension upwards from the base. A positive number from 0 to 1 curved to be
|
||||
// maximal at 0. Also to be multiplied by CONFIG.mountain_scale.
|
||||
let alt_main = (gen_ctx
|
||||
.alt_nz
|
||||
.get((wposf.div(2_000.0)).into_array())
|
||||
.min(1.0)
|
||||
.max(-1.0))
|
||||
.abs()
|
||||
.powf(1.35);
|
||||
|
||||
fn spring(x: f64, pow: f64) -> f64 {
|
||||
x.abs().powf(pow) * x.signum()
|
||||
}
|
||||
|
||||
(0.0 + alt_main
|
||||
+ (gen_ctx
|
||||
.small_nz
|
||||
.get((wposf.div(300.0)).into_array())
|
||||
.min(1.0)
|
||||
.max(-1.0))
|
||||
.mul(alt_main.powf(0.8).max(/*0.25*/ 0.15))
|
||||
.mul(0.3)
|
||||
.add(1.0)
|
||||
.mul(0.4)
|
||||
/* + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0)
|
||||
.mul(0.045)*/)
|
||||
};
|
||||
old_height_uniform(posi)
|
||||
/* // 0.0
|
||||
// -/*CONFIG.sea_level / CONFIG.mountain_scale*//* 0.75 */1.0
|
||||
// ((old_height(posi) - alt_old_min) as f64 / (alt_old_max - alt_old_min) as f64) as f32
|
||||
// uplift_fn(posi) / max_erosion_per_delta_t as f32
|
||||
// old_height_uniform(posi) *
|
||||
(/*((old_height(posi) - alt_old_min) as f64 / (alt_old_max - alt_old_min) as f64) **/(((6.0 / 360.0 * 2.0 * f64::consts::PI).tan()
|
||||
* TerrainChunkSize::RECT_SIZE.reduce_partial_min() as f64)
|
||||
.floor()
|
||||
* height_scale)) as f32
|
||||
// 5.0 / CONFIG.mountain_scale */
|
||||
}
|
||||
};
|
||||
let (alt, basement) = do_erosion(
|
||||
0.0,
|
||||
max_erosion_per_delta_t as f32,
|
||||
n_steps,
|
||||
&river_seed,
|
||||
&rock_strength_nz,
|
||||
|posi| alt_func(posi),
|
||||
|posi| alt_func(posi),// if is_ocean_fn(posi) { old_height(posi) } else { 0.0 },
|
||||
is_ocean_fn,
|
||||
uplift_fn,
|
||||
);
|
||||
// Quick "small scale" erosion cycle in order to lower extreme angles.
|
||||
let (alt, basement) = do_erosion(
|
||||
0.0,
|
||||
(1.0 * height_scale) as f32,
|
||||
n_small_steps,
|
||||
&river_seed,
|
||||
&rock_strength_nz,
|
||||
|posi| /* if is_ocean_fn(posi) { old_height(posi) } else { alt[posi] } */alt[posi],
|
||||
|posi| basement[posi],
|
||||
is_ocean_fn,
|
||||
|posi| uplift_fn(posi) * (1.0 * height_scale / max_erosion_per_delta_t) as f32,
|
||||
);
|
||||
|
||||
let is_ocean = get_oceans(|posi| alt[posi]);
|
||||
let is_ocean_fn = |posi: usize| is_ocean[posi];
|
||||
let mut dh = downhill(&alt, /*old_height*/ is_ocean_fn);
|
||||
@ -765,6 +885,7 @@ impl WorldSim {
|
||||
temp_base,
|
||||
chaos,
|
||||
alt,
|
||||
basement,
|
||||
water_alt,
|
||||
dh,
|
||||
flux: flux_old,
|
||||
@ -786,7 +907,7 @@ impl WorldSim {
|
||||
rng,
|
||||
};
|
||||
|
||||
this.seed_elements();
|
||||
// this.seed_elements();
|
||||
|
||||
this
|
||||
}
|
||||
@ -1199,6 +1320,7 @@ impl WorldSim {
|
||||
pub struct SimChunk {
|
||||
pub chaos: f32,
|
||||
pub alt: f32,
|
||||
pub basement: f32,
|
||||
pub water_alt: f32,
|
||||
pub downhill: Option<Vec2<i32>>,
|
||||
pub flux: f32,
|
||||
@ -1244,6 +1366,7 @@ impl SimChunk {
|
||||
let _map_edge_factor = map_edge_factor(posi);
|
||||
let (_, chaos) = gen_cdf.chaos[posi];
|
||||
let alt_pre = gen_cdf.alt[posi];
|
||||
let basement_pre = gen_cdf.basement[posi];
|
||||
let water_alt_pre = gen_cdf.water_alt[posi];
|
||||
let downhill_pre = gen_cdf.dh[posi];
|
||||
let flux = gen_cdf.flux[posi];
|
||||
@ -1300,10 +1423,12 @@ impl SimChunk {
|
||||
panic!("Halp!");
|
||||
} */
|
||||
|
||||
let mut alt = CONFIG.sea_level.add(alt_pre.mul(CONFIG.mountain_scale));
|
||||
let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale;
|
||||
let mut alt = CONFIG.sea_level.add(alt_pre.div(height_scale));
|
||||
let mut basement = CONFIG.sea_level.add(basement_pre.div(height_scale));
|
||||
let water_alt = CONFIG
|
||||
.sea_level
|
||||
.add(water_alt_pre.mul(CONFIG.mountain_scale));
|
||||
.add(water_alt_pre.div(height_scale));
|
||||
let downhill = if downhill_pre == -2 {
|
||||
None
|
||||
} else if downhill_pre < 0 {
|
||||
@ -1391,6 +1516,7 @@ impl SimChunk {
|
||||
chaos,
|
||||
flux,
|
||||
alt,
|
||||
basement: basement.min(alt),
|
||||
water_alt,
|
||||
downhill,
|
||||
temp,
|
||||
|
||||
@ -1,9 +1,10 @@
|
||||
use super::WORLD_SIZE;
|
||||
use bitvec::prelude::{bitbox, bitvec, BitBox};
|
||||
use common::{terrain::TerrainChunkSize, vol::RectVolSize};
|
||||
use noise::{MultiFractal, NoiseFn, Perlin, Point2, Point3, Point4, Seedable};
|
||||
use num::Float;
|
||||
use rayon::prelude::*;
|
||||
use std::{f32, f64, u32};
|
||||
use std::{f32, f64, ops::Mul, u32};
|
||||
use vek::*;
|
||||
|
||||
/// Calculates the smallest distance along an axis (x, y) from an edge of
|
||||
@ -301,3 +302,392 @@ pub fn get_oceans(oldh: impl Fn(usize) -> f32 + Sync) -> BitBox {
|
||||
}
|
||||
is_ocean
|
||||
}
|
||||
|
||||
/// A 2-dimensional vector, for internal use.
|
||||
type Vector2<T> = [T; 2];
|
||||
/// A 3-dimensional vector, for internal use.
|
||||
type Vector3<T> = [T; 3];
|
||||
/// A 4-dimensional vector, for internal use.
|
||||
type Vector4<T> = [T; 4];
|
||||
|
||||
#[inline]
|
||||
fn zip_with2<T, U, V, F>(a: Vector2<T>, b: Vector2<U>, f: F) -> Vector2<V>
|
||||
where
|
||||
T: Copy,
|
||||
U: Copy,
|
||||
F: Fn(T, U) -> V,
|
||||
{
|
||||
let (ax, ay) = (a[0], a[1]);
|
||||
let (bx, by) = (b[0], b[1]);
|
||||
[f(ax, bx), f(ay, by)]
|
||||
}
|
||||
|
||||
#[inline]
|
||||
fn zip_with3<T, U, V, F>(a: Vector3<T>, b: Vector3<U>, f: F) -> Vector3<V>
|
||||
where
|
||||
T: Copy,
|
||||
U: Copy,
|
||||
F: Fn(T, U) -> V,
|
||||
{
|
||||
let (ax, ay, az) = (a[0], a[1], a[2]);
|
||||
let (bx, by, bz) = (b[0], b[1], b[2]);
|
||||
[f(ax, bx), f(ay, by), f(az, bz)]
|
||||
}
|
||||
|
||||
#[inline]
|
||||
fn zip_with4<T, U, V, F>(a: Vector4<T>, b: Vector4<U>, f: F) -> Vector4<V>
|
||||
where
|
||||
T: Copy,
|
||||
U: Copy,
|
||||
F: Fn(T, U) -> V,
|
||||
{
|
||||
let (ax, ay, az, aw) = (a[0], a[1], a[2], a[3]);
|
||||
let (bx, by, bz, bw) = (b[0], b[1], b[2], b[3]);
|
||||
[f(ax, bx), f(ay, by), f(az, bz), f(aw, bw)]
|
||||
}
|
||||
|
||||
#[inline]
|
||||
fn mul2<T>(a: Vector2<T>, b: T) -> Vector2<T>
|
||||
where
|
||||
T: Copy + Mul<T, Output = T>,
|
||||
{
|
||||
zip_with2(a, const2(b), Mul::mul)
|
||||
}
|
||||
|
||||
#[inline]
|
||||
fn mul3<T>(a: Vector3<T>, b: T) -> Vector3<T>
|
||||
where
|
||||
T: Copy + Mul<T, Output = T>,
|
||||
{
|
||||
zip_with3(a, const3(b), Mul::mul)
|
||||
}
|
||||
|
||||
#[inline]
|
||||
fn mul4<T>(a: Vector4<T>, b: T) -> Vector4<T>
|
||||
where
|
||||
T: Copy + Mul<T, Output = T>,
|
||||
{
|
||||
zip_with4(a, const4(b), Mul::mul)
|
||||
}
|
||||
|
||||
#[inline]
|
||||
fn const2<T: Copy>(x: T) -> Vector2<T> {
|
||||
[x, x]
|
||||
}
|
||||
|
||||
#[inline]
|
||||
fn const3<T: Copy>(x: T) -> Vector3<T> {
|
||||
[x, x, x]
|
||||
}
|
||||
|
||||
#[inline]
|
||||
fn const4<T: Copy>(x: T) -> Vector4<T> {
|
||||
[x, x, x, x]
|
||||
}
|
||||
|
||||
fn build_sources(seed: u32, octaves: usize) -> Vec<Perlin> {
|
||||
let mut sources = Vec::with_capacity(octaves);
|
||||
for x in 0..octaves {
|
||||
sources.push(Perlin::new().set_seed(seed + x as u32));
|
||||
}
|
||||
sources
|
||||
}
|
||||
|
||||
/// Noise function that outputs hybrid Multifractal noise.
|
||||
///
|
||||
/// The result of this multifractal noise is that valleys in the noise should
|
||||
/// have smooth bottoms at all altitudes.
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct HybridMulti {
|
||||
/// Total number of frequency octaves to generate the noise with.
|
||||
///
|
||||
/// The number of octaves control the _amount of detail_ in the noise
|
||||
/// function. Adding more octaves increases the detail, with the drawback
|
||||
/// of increasing the calculation time.
|
||||
pub octaves: usize,
|
||||
|
||||
/// The number of cycles per unit length that the noise function outputs.
|
||||
pub frequency: f64,
|
||||
|
||||
/// A multiplier that determines how quickly the frequency increases for
|
||||
/// each successive octave in the noise function.
|
||||
///
|
||||
/// The frequency of each successive octave is equal to the product of the
|
||||
/// previous octave's frequency and the lacunarity value.
|
||||
///
|
||||
/// A lacunarity of 2.0 results in the frequency doubling every octave. For
|
||||
/// almost all cases, 2.0 is a good value to use.
|
||||
pub lacunarity: f64,
|
||||
|
||||
/// A multiplier that determines how quickly the amplitudes diminish for
|
||||
/// each successive octave in the noise function.
|
||||
///
|
||||
/// The amplitude of each successive octave is equal to the product of the
|
||||
/// previous octave's amplitude and the persistence value. Increasing the
|
||||
/// persistence produces "rougher" noise.
|
||||
///
|
||||
/// H = 1.0 - fractal increment = -ln(persistence) / ln(lacunarity). For
|
||||
/// a fractal increment between 0 (inclusive) and 1 (exclusive), keep
|
||||
/// persistence between 1 / lacunarity (inclusive, for low fractal
|
||||
/// dimension) and 1 (exclusive, for high fractal dimension).
|
||||
pub persistence: f64,
|
||||
|
||||
/// An offset that is added to the output of each sample of the underlying
|
||||
/// Perlin noise function. Because each successive octave is weighted in
|
||||
/// part by the previous signal's output, increasing the offset will weight
|
||||
/// the output more heavily towards 1.0.
|
||||
pub offset: f64,
|
||||
|
||||
seed: u32,
|
||||
sources: Vec<Perlin>,
|
||||
}
|
||||
|
||||
impl HybridMulti {
|
||||
pub const DEFAULT_SEED: u32 = 0;
|
||||
pub const DEFAULT_OCTAVES: usize = 6;
|
||||
pub const DEFAULT_FREQUENCY: f64 = 2.0;
|
||||
pub const DEFAULT_LACUNARITY: f64 = /* std::f64::consts::PI * 2.0 / 3.0 */2.0;
|
||||
// -ln(2^(-0.25))/ln(2) = 0.25
|
||||
// 2^(-0.25) ~ 13/16
|
||||
pub const DEFAULT_PERSISTENCE: f64 = /* 0.25 *//* 0.5*/ 13.0 / 16.0;
|
||||
pub const DEFAULT_OFFSET: f64 = /* 0.25 *//* 0.5*/ 0.7;
|
||||
pub const MAX_OCTAVES: usize = 32;
|
||||
|
||||
pub fn new() -> Self {
|
||||
Self {
|
||||
seed: Self::DEFAULT_SEED,
|
||||
octaves: Self::DEFAULT_OCTAVES,
|
||||
frequency: Self::DEFAULT_FREQUENCY,
|
||||
lacunarity: Self::DEFAULT_LACUNARITY,
|
||||
persistence: Self::DEFAULT_PERSISTENCE,
|
||||
offset: Self::DEFAULT_OFFSET,
|
||||
sources: build_sources(Self::DEFAULT_SEED, Self::DEFAULT_OCTAVES),
|
||||
}
|
||||
}
|
||||
|
||||
pub fn set_offset(self, offset: f64) -> Self {
|
||||
Self { offset, ..self }
|
||||
}
|
||||
}
|
||||
|
||||
impl Default for HybridMulti {
|
||||
fn default() -> Self {
|
||||
Self::new()
|
||||
}
|
||||
}
|
||||
|
||||
impl MultiFractal for HybridMulti {
|
||||
fn set_octaves(self, mut octaves: usize) -> Self {
|
||||
if self.octaves == octaves {
|
||||
return self;
|
||||
}
|
||||
|
||||
octaves = octaves.max(1).min(Self::MAX_OCTAVES);
|
||||
Self {
|
||||
octaves,
|
||||
sources: build_sources(self.seed, octaves),
|
||||
..self
|
||||
}
|
||||
}
|
||||
|
||||
fn set_frequency(self, frequency: f64) -> Self {
|
||||
Self { frequency, ..self }
|
||||
}
|
||||
|
||||
fn set_lacunarity(self, lacunarity: f64) -> Self {
|
||||
Self { lacunarity, ..self }
|
||||
}
|
||||
|
||||
fn set_persistence(self, persistence: f64) -> Self {
|
||||
Self {
|
||||
persistence,
|
||||
..self
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl Seedable for HybridMulti {
|
||||
fn set_seed(self, seed: u32) -> Self {
|
||||
if self.seed == seed {
|
||||
return self;
|
||||
}
|
||||
|
||||
Self {
|
||||
seed,
|
||||
sources: build_sources(seed, self.octaves),
|
||||
..self
|
||||
}
|
||||
}
|
||||
|
||||
fn seed(&self) -> u32 {
|
||||
self.seed
|
||||
}
|
||||
}
|
||||
|
||||
/// 2-dimensional `HybridMulti` noise
|
||||
impl NoiseFn<Point2<f64>> for HybridMulti {
|
||||
fn get(&self, mut point: Point2<f64>) -> f64 {
|
||||
// First unscaled octave of function; later octaves are scaled.
|
||||
point = mul2(point, self.frequency);
|
||||
// Offset and bias to scale into [offset - 1.0, 1.0 + offset] range.
|
||||
let bias = 1.0;
|
||||
let mut result = (self.sources[0].get(point) + self.offset) * bias * self.persistence;
|
||||
let mut scale = self.persistence;
|
||||
let mut weight = result;
|
||||
|
||||
// Spectral construction inner loop, where the fractal is built.
|
||||
for x in 1..self.octaves {
|
||||
// Prevent divergence.
|
||||
weight = weight.min(1.0);
|
||||
|
||||
// Raise the spatial frequency.
|
||||
point = mul2(point, self.lacunarity);
|
||||
|
||||
// Get noise value, and scale it to the [offset - 1.0, 1.0 + offset] range.
|
||||
let mut signal = (self.sources[x].get(point) + self.offset) * bias;
|
||||
|
||||
// Scale the amplitude appropriately for this frequency.
|
||||
let exp_scale = self.persistence.powi(x as i32);
|
||||
signal *= exp_scale;
|
||||
|
||||
// Add it in, weighted by previous octave's noise value.
|
||||
result += weight * signal;
|
||||
|
||||
// Update the weighting value.
|
||||
weight *= signal;
|
||||
scale += exp_scale;
|
||||
}
|
||||
|
||||
// Scale the result to the [-1,1] range
|
||||
(result / scale) / bias - self.offset
|
||||
}
|
||||
}
|
||||
|
||||
/// 3-dimensional `HybridMulti` noise
|
||||
impl NoiseFn<Point3<f64>> for HybridMulti {
|
||||
fn get(&self, mut point: Point3<f64>) -> f64 {
|
||||
// First unscaled octave of function; later octaves are scaled.
|
||||
point = mul3(point, self.frequency);
|
||||
// Offset and bias to scale into [offset - 1.0, 1.0 + offset] range.
|
||||
let bias = 1.0;
|
||||
let mut result = (self.sources[0].get(point) + self.offset) * bias * self.persistence;
|
||||
let mut scale = self.persistence;
|
||||
let mut weight = result;
|
||||
|
||||
// Spectral construction inner loop, where the fractal is built.
|
||||
for x in 1..self.octaves {
|
||||
// Prevent divergence.
|
||||
weight = weight.min(1.0);
|
||||
|
||||
// Raise the spatial frequency.
|
||||
point = mul3(point, self.lacunarity);
|
||||
|
||||
// Get noise value, and scale it to the [0, 1.0] range.
|
||||
let mut signal = (self.sources[x].get(point) + self.offset) * bias;
|
||||
|
||||
// Scale the amplitude appropriately for this frequency.
|
||||
let exp_scale = self.persistence.powi(x as i32);
|
||||
signal *= exp_scale;
|
||||
|
||||
// Add it in, weighted by previous octave's noise value.
|
||||
result += weight * signal;
|
||||
|
||||
// Update the weighting value.
|
||||
weight *= signal;
|
||||
scale += exp_scale;
|
||||
}
|
||||
|
||||
// Scale the result to the [-1,1] range
|
||||
(result / scale) / bias - self.offset
|
||||
}
|
||||
}
|
||||
|
||||
/// 4-dimensional `HybridMulti` noise
|
||||
impl NoiseFn<Point4<f64>> for HybridMulti {
|
||||
fn get(&self, mut point: Point4<f64>) -> f64 {
|
||||
// First unscaled octave of function; later octaves are scaled.
|
||||
point = mul4(point, self.frequency);
|
||||
// Offset and bias to scale into [offset - 1.0, 1.0 + offset] range.
|
||||
let bias = 1.0;
|
||||
let mut result = (self.sources[0].get(point) + self.offset) * bias * self.persistence;
|
||||
let mut exp_scale = self.persistence;
|
||||
let mut scale = self.persistence;
|
||||
let mut weight = result;
|
||||
|
||||
// Spectral construction inner loop, where the fractal is built.
|
||||
for x in 1..self.octaves {
|
||||
// Prevent divergence.
|
||||
weight = weight.min(1.0);
|
||||
|
||||
// Raise the spatial frequency.
|
||||
point = mul4(point, self.lacunarity);
|
||||
|
||||
// Get noise value, and scale it to the [0, 1.0] range.
|
||||
let mut signal = (self.sources[x].get(point) + self.offset) * bias;
|
||||
|
||||
// Scale the amplitude appropriately for this frequency.
|
||||
exp_scale *= self.persistence;
|
||||
signal *= exp_scale;
|
||||
|
||||
// Add it in, weighted by previous octave's noise value.
|
||||
result += weight * signal;
|
||||
|
||||
// Update the weighting value.
|
||||
weight *= signal;
|
||||
scale += exp_scale;
|
||||
}
|
||||
|
||||
// Scale the result to the [-1,1] range
|
||||
(result / scale) / bias - self.offset
|
||||
}
|
||||
}
|
||||
|
||||
/// Noise function that applies a scaling factor and a bias to the output value
|
||||
/// from the source function.
|
||||
///
|
||||
/// The function retrieves the output value from the source function, multiplies
|
||||
/// it with the scaling factor, adds the bias to it, then outputs the value.
|
||||
pub struct ScaleBias<'a, F: 'a> {
|
||||
/// Outputs a value.
|
||||
pub source: &'a F,
|
||||
|
||||
/// Scaling factor to apply to the output value from the source function.
|
||||
/// The default value is 1.0.
|
||||
pub scale: f64,
|
||||
|
||||
/// Bias to apply to the scaled output value from the source function.
|
||||
/// The default value is 0.0.
|
||||
pub bias: f64,
|
||||
}
|
||||
|
||||
impl<'a, F> ScaleBias<'a, F> {
|
||||
pub fn new(source: &'a F) -> Self {
|
||||
ScaleBias {
|
||||
source,
|
||||
scale: 1.0,
|
||||
bias: 0.0,
|
||||
}
|
||||
}
|
||||
|
||||
pub fn set_scale(self, scale: f64) -> Self {
|
||||
ScaleBias { scale, ..self }
|
||||
}
|
||||
|
||||
pub fn set_bias(self, bias: f64) -> Self {
|
||||
ScaleBias { bias, ..self }
|
||||
}
|
||||
}
|
||||
|
||||
impl<'a, F: NoiseFn<T> + 'a, T> NoiseFn<T> for ScaleBias<'a, F> {
|
||||
#[cfg(not(target_os = "emscripten"))]
|
||||
fn get(&self, point: T) -> f64 {
|
||||
(self.source.get(point)).mul_add(self.scale, self.bias)
|
||||
}
|
||||
|
||||
#[cfg(target_os = "emscripten")]
|
||||
fn get(&self, point: T) -> f64 {
|
||||
(self.source.get(point) * self.scale) + self.bias
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user