nx_wind_erosion.nx source
↩ module page · 250 lines · 12253 B
1// nx_wind_erosion.nx -- directional anisotropic erosion generator.
2//
3// Fourth demonstration of the kind-specific-generator cardinal
4// `feedback-kind-specific-generators-not-broad-noise`. Previous
5// shipped kind-generators:
6// nx_forest_layout -- scattered features (Poisson-like)
7// nx_river_carve -- polyline path
8// nx_beach_zone -- banded zone
9// This one:
10// nx_wind_erosion -- directional anisotropic vector field
11//
12// Wind doesn't erode terrain uniformly. A south-facing slope under
13// north wind gets MORE scour than a north-facing slope under the
14// same wind. And on the leeward side, blown-particle deposition
15// builds dunes. This primitive captures that anisotropy.
16//
17// Wind-internal logic in v1:
18// - 4 cardinal wind directions (NORTH/EAST/SOUTH/WEST) for an
19// anisotropic dot-product test. Caller supplies LOCAL SLOPE
20// (dh/dx, dh/dz) at the query point.
21// - Windward factor = dot product of wind direction with slope.
22// Positive = facing into the wind (exposed) -> erode.
23// Negative = facing downwind (sheltered) -> deposit dunes.
24// - Two strength params (scour + dune) so caller tunes erosion
25// and deposition independently.
26//
27// V2 follow-ups (queued):
28// - 8 cardinal directions (add NE/SE/SW/NW with diagonal slopes)
29// - Arbitrary wind angle in degrees + sin/cos lookup via nx_camera_q14
30// - Multi-step wind history (prevailing wind + occasional storms)
31// - Wind speed / strength field varying over space (terrain shadow)
32// - Particle transport modeling (sand grains move from windward to
33// leeward; equilibrium dunes emerge)
34// - Composes with nx_beach_zone v2 dune-barrier model
35// - Composes with nx_weather_pattern for storm-driven episodic events
36//
37// Loss audit: Q14 integer arithmetic; dot product is exact. Signed
38// delta scales the strength params linearly; no information loss
39// beyond Q14 quantisation.
40//
41// genealogy_id: bagnold_1941_physics_of_blown_sand +
42// sherman_1995_wind_blown_sand_modelling +
43// nasa_mars_aeolian_processes_canon
44// lineage_id: nx_wind_erosion_anisotropic_q14_v1
45//
46// nx_safety_envelope:
47// intended_use: "Anisotropic wind-erosion sim -- deep-time
48// terrain weathering in directional regime"
49// sil_target: SIL1
50// asil_target: QM
51// dal_target: NONE
52// evidence: [Q14_fixed_point, anisotropy_axis_documented,
53// bounded_cellular_passes]
54// hazard_register: [bug-tape-wind-direction-coordinate-confusion,
55// bug-tape-deposit-mass-leak]
56// residual_risk: "Procgen quality only."
57// verdict: NOT_YET_EVALUATED
58
59import "nx_syscalls.nx"
60import "nx_tier.nx"
61
62// ===== Q14 ==========================================================
63const NX_WIND_Q: nx_int = 16384
64
65// ===== Wind-direction sealed enum ===================================
66// 4 cardinals for v1. Direction names the FROM-direction (the
67// direction the wind is BLOWING FROM):
68// NORTH = wind from +Z, blowing toward -Z
69// EAST = wind from +X, blowing toward -X
70// SOUTH = wind from -Z, blowing toward +Z
71// WEST = wind from -X, blowing toward +X
72const NX_WIND_DIR_NORTH: nx_int = 0
73const NX_WIND_DIR_EAST: nx_int = 1
74const NX_WIND_DIR_SOUTH: nx_int = 2
75const NX_WIND_DIR_WEST: nx_int = 3
76
77const NX_WIND_DIR_COUNT: nx_int = 4
78
79// ===== Validity predicate ===========================================
80func nx_wind_dir_is_valid(d: nx_int) -> nx_int {
81 if d == NX_WIND_DIR_NORTH { return 1 }
82 if d == NX_WIND_DIR_EAST { return 1 }
83 if d == NX_WIND_DIR_SOUTH { return 1 }
84 if d == NX_WIND_DIR_WEST { return 1 }
85 return 0
86}
87
88// ===== Windward-factor dot product =================================
89// Returns the slope component along the wind direction (= how
90// "windward" this slope is).
91//
92// Input slopes (dh/dx, dh/dz) are Q14 metres-per-Q14-distance, i.e.
93// pure Q14 ratios. Wind direction is "where wind comes from" so
94// "windward face" = uphill in the FROM direction = downhill in the
95// TO direction. A slope rising IN THE WIND'S PATH is windward.
96//
97// For NORTH wind (from +Z, blowing toward -Z):
98// windward = slope rises toward +Z = +dh/dz
99// The dot product of wind_to_direction (0, -1) with slope-gradient
100// (dh/dx, dh/dz) = -dh/dz... wait that's wrong sign for "windward
101// into wind."
102//
103// Convention: "windward" is the side FACING the wind. If wind comes
104// FROM +Z and slope rises toward +Z, that slope FACES the wind
105// (windward). In terms of slope gradient:
106// slope_gradient_z > 0 means height increases as Z increases (so
107// the upward direction has +Z component)
108// wind from +Z means windward-facing-normal points toward +Z
109// Therefore windward factor for NORTH wind = +slope_gradient_z.
110//
111// Symmetry for the other 3:
112// EAST wind (from +X): windward = +slope_gradient_x
113// SOUTH wind (from -Z): windward = -slope_gradient_z
114// WEST wind (from -X): windward = -slope_gradient_x
115//
116// Returns Q14 in [-2Q, +2Q] roughly (depends on caller's slope
117// magnitudes). Sign: positive = windward, negative = leeward.
118func nx_wind_dir_dot(
119 wind_dir: nx_int,
120 slope_dx_q14: nx_int,
121 slope_dz_q14: nx_int
122) -> nx_int {
123 if wind_dir == NX_WIND_DIR_NORTH { return slope_dz_q14 }
124 if wind_dir == NX_WIND_DIR_EAST { return slope_dx_q14 }
125 if wind_dir == NX_WIND_DIR_SOUTH { return 0 - slope_dz_q14 }
126 if wind_dir == NX_WIND_DIR_WEST { return 0 - slope_dx_q14 }
127 return 0
128}
129
130// ===== Erosion delta =================================================
131// Given a windward factor (output of nx_wind_dir_dot) and two
132// strength params, returns signed height delta in Q14 metres:
133// windward (factor > 0): scour, NEGATIVE delta (terrain eroded)
134// delta = -scour_strength * factor / Q
135// leeward (factor < 0): dune, POSITIVE delta (terrain built up)
136// delta = -dune_strength * factor / Q = dune_strength * |factor| / Q
137// flat (factor == 0): no change
138//
139// scour_strength and dune_strength are PEAK magnitudes; the actual
140// delta scales by the windward factor. Strong slopes lose more
141// height; weak slopes lose less.
142//
143// Caller adds the returned delta to their base heightmap.
144func nx_wind_erosion_delta(
145 windward_factor_q14: nx_int,
146 scour_strength_q14_m: nx_int,
147 dune_strength_q14_m: nx_int
148) -> nx_int {
149 if windward_factor_q14 == 0 { return 0 }
150 if windward_factor_q14 > 0 {
151 return 0 - scour_strength_q14_m * windward_factor_q14 / NX_WIND_Q
152 }
153 // Leeward: factor is negative; result is positive.
154 return 0 - dune_strength_q14_m * windward_factor_q14 / NX_WIND_Q
155}
156
157// ===== Combined entry =============================================
158// Convenience: compute windward factor + delta in one call. Useful
159// when the caller has slope gradients directly and just wants the
160// erosion-delta result.
161func nx_wind_erosion_at(
162 wind_dir: nx_int,
163 slope_dx_q14: nx_int,
164 slope_dz_q14: nx_int,
165 scour_strength_q14_m: nx_int,
166 dune_strength_q14_m: nx_int
167) -> nx_int {
168 if nx_wind_dir_is_valid(wind_dir) == 0 { return 0 }
169 let factor: nx_int = nx_wind_dir_dot(wind_dir, slope_dx_q14, slope_dz_q14)
170 return nx_wind_erosion_delta(factor, scour_strength_q14_m, dune_strength_q14_m)
171}
172
173// ===== Self-test ====================================================
174func main() -> i64 {
175 let q: nx_int = NX_WIND_Q
176
177 // T1: Validity predicate.
178 if nx_wind_dir_is_valid(NX_WIND_DIR_NORTH) != 1 { return __syscall(93, 1, 0, 0, 0, 0, 0) }
179 if nx_wind_dir_is_valid(NX_WIND_DIR_EAST) != 1 { return __syscall(93, 2, 0, 0, 0, 0, 0) }
180 if nx_wind_dir_is_valid(NX_WIND_DIR_SOUTH) != 1 { return __syscall(93, 3, 0, 0, 0, 0, 0) }
181 if nx_wind_dir_is_valid(NX_WIND_DIR_WEST) != 1 { return __syscall(93, 4, 0, 0, 0, 0, 0) }
182 if nx_wind_dir_is_valid(99) != 0 { return __syscall(93, 5, 0, 0, 0, 0, 0) }
183
184 // T2: Dot product per direction. Slope rising toward +Z (dz=Q,
185 // dx=0).
186 if nx_wind_dir_dot(NX_WIND_DIR_NORTH, 0, q) != q { return __syscall(93, 10, 0, 0, 0, 0, 0) }
187 if nx_wind_dir_dot(NX_WIND_DIR_SOUTH, 0, q) != (0 - q) { return __syscall(93, 11, 0, 0, 0, 0, 0) }
188 if nx_wind_dir_dot(NX_WIND_DIR_EAST, 0, q) != 0 { return __syscall(93, 12, 0, 0, 0, 0, 0) }
189 if nx_wind_dir_dot(NX_WIND_DIR_WEST, 0, q) != 0 { return __syscall(93, 13, 0, 0, 0, 0, 0) }
190 // Slope rising toward +X (dx=Q, dz=0).
191 if nx_wind_dir_dot(NX_WIND_DIR_EAST, q, 0) != q { return __syscall(93, 14, 0, 0, 0, 0, 0) }
192 if nx_wind_dir_dot(NX_WIND_DIR_WEST, q, 0) != (0 - q) { return __syscall(93, 15, 0, 0, 0, 0, 0) }
193 if nx_wind_dir_dot(NX_WIND_DIR_NORTH, q, 0) != 0 { return __syscall(93, 16, 0, 0, 0, 0, 0) }
194 // Negative slope (dz=-Q) -- north wind sees it as leeward.
195 if nx_wind_dir_dot(NX_WIND_DIR_NORTH, 0, 0 - q) != (0 - q) { return __syscall(93, 17, 0, 0, 0, 0, 0) }
196 // Flat ground (dx=0, dz=0): all directions give 0.
197 if nx_wind_dir_dot(NX_WIND_DIR_NORTH, 0, 0) != 0 { return __syscall(93, 18, 0, 0, 0, 0, 0) }
198 if nx_wind_dir_dot(NX_WIND_DIR_SOUTH, 0, 0) != 0 { return __syscall(93, 19, 0, 0, 0, 0, 0) }
199
200 // T3: Erosion delta -- pure windward (full Q factor) -> -scour.
201 if nx_wind_erosion_delta(q, 100, 50) != (0 - 100) { return __syscall(93, 30, 0, 0, 0, 0, 0) }
202 // Pure leeward (full -Q factor) -> +dune.
203 if nx_wind_erosion_delta(0 - q, 100, 50) != 50 { return __syscall(93, 31, 0, 0, 0, 0, 0) }
204 // Flat -> 0.
205 if nx_wind_erosion_delta(0, 100, 50) != 0 { return __syscall(93, 32, 0, 0, 0, 0, 0) }
206 // Half-windward -> -scour/2.
207 if nx_wind_erosion_delta(q / 2, 100, 50) != (0 - 50) { return __syscall(93, 33, 0, 0, 0, 0, 0) }
208 // Half-leeward -> +dune/2.
209 if nx_wind_erosion_delta(0 - q / 2, 100, 50) != 25 { return __syscall(93, 34, 0, 0, 0, 0, 0) }
210
211 // T4: Combined entry -- N wind on +Z-rising slope -> windward -> scour.
212 let d_n_up: nx_int = nx_wind_erosion_at(NX_WIND_DIR_NORTH, 0, q, 100, 50)
213 if d_n_up != (0 - 100) { return __syscall(93, 40, 0, 0, 0, 0, 0) }
214 // Same slope, S wind -> leeward -> dune build.
215 let d_s_up: nx_int = nx_wind_erosion_at(NX_WIND_DIR_SOUTH, 0, q, 100, 50)
216 if d_s_up != 50 { return __syscall(93, 41, 0, 0, 0, 0, 0) }
217 // E wind: slope is perpendicular -> 0 delta.
218 let d_e_up: nx_int = nx_wind_erosion_at(NX_WIND_DIR_EAST, 0, q, 100, 50)
219 if d_e_up != 0 { return __syscall(93, 42, 0, 0, 0, 0, 0) }
220
221 // T5: Anisotropy verification -- the SAME slope under TWO
222 // different wind directions yields different deltas. A
223 // (dx=Q, dz=0) slope under E wind = full scour; under W wind =
224 // full dune-build. Under N or S wind = 0.
225 let slope_e: *i64 = (sys_mmap(2 * NX_SIZEOF_NX_INT)) as *i64
226 slope_e[0] = q
227 slope_e[1] = 0
228 let d_ew: nx_int = nx_wind_erosion_at(NX_WIND_DIR_EAST, slope_e[0], slope_e[1], 100, 50)
229 let d_ww: nx_int = nx_wind_erosion_at(NX_WIND_DIR_WEST, slope_e[0], slope_e[1], 100, 50)
230 let d_nw: nx_int = nx_wind_erosion_at(NX_WIND_DIR_NORTH, slope_e[0], slope_e[1], 100, 50)
231 if d_ew != (0 - 100) { return __syscall(93, 50, 0, 0, 0, 0, 0) }
232 if d_ww != 50 { return __syscall(93, 51, 0, 0, 0, 0, 0) }
233 if d_nw != 0 { return __syscall(93, 52, 0, 0, 0, 0, 0) }
234 // Inverted slope (dx=-Q) flips windward<->leeward for E/W.
235 let d_e_inv: nx_int = nx_wind_erosion_at(NX_WIND_DIR_EAST, 0 - q, 0, 100, 50)
236 let d_w_inv: nx_int = nx_wind_erosion_at(NX_WIND_DIR_WEST, 0 - q, 0, 100, 50)
237 if d_e_inv != 50 { return __syscall(93, 53, 0, 0, 0, 0, 0) }
238 if d_w_inv != (0 - 100) { return __syscall(93, 54, 0, 0, 0, 0, 0) }
239
240 // T6: Refusal -- invalid wind_dir to nx_wind_erosion_at -> 0.
241 if nx_wind_erosion_at(99, q, q, 100, 50) != 0 { return __syscall(93, 60, 0, 0, 0, 0, 0) }
242
243 // T7: Scaling with strength. Doubling scour doubles the
244 // negative delta.
245 let d_base: nx_int = nx_wind_erosion_delta(q, 100, 50)
246 let d_2x: nx_int = nx_wind_erosion_delta(q, 200, 50)
247 if d_2x != 2 * d_base { return __syscall(93, 70, 0, 0, 0, 0, 0) }
248
249 return 0
250}