1% Smooth random function in 3D — chebfun's randnfun3, as the Fourier modes
2% it is built from rather than as a chebfun3.
3%
4% [K, C] = randnfun3(LAMBDA, DOM) draws a random trig series on the box
5% DOM = [x0 x1 y0 y1 z0 z1] with maximum frequency about 2*pi/LAMBDA in
6% each direction and standard normal distribution N(0,1) at each point.
7% K is nmodes x 3 (angular wavenumbers) and C is nmodes x 2 (real and
8% imaginary parts), defining
9%
10% f(x,y,z) = sum_j C(j,1)*cos(K(j,:)*[x;y;z]) - C(j,2)*sin(K(j,:)*[x;y;z])
11%
12% Seed the draw with rng(...) before calling.
13%
14% chebfun returns a chebfun3 and evaluates it later; this project has no
15% such object, and the sum above is what the GPU evaluates at the surface
16% points (src/mgpu/randnfun3.ts). Splitting it here is also what keeps the
17% draw in MATLAB: randn has no counterpart in the compiled WGSL dialect.
19function [k, c] = randnfun3(lambda, dom)
20 % chebfun's nonperiodic path builds a periodic function on a domain about
21 % 20% larger and restricts it. Restriction is free when evaluating at
22 % points, so we keep the enlarged period and never form the smaller one.
23 m = round(1.2*(dom(2)-dom(1))/lambda + 2);
24 n = round(1.2*(dom(4)-dom(3))/lambda + 2);
25 p = round(1.2*(dom(6)-dom(5))/lambda + 2);
26 m2 = 2*m+1;
27 n2 = 2*n+1;
28 p2 = 2*p+1;
29 N = m2*n2*p2;
31 % chebfun draws the whole cube (column-major) before masking; drawing in
32 % that same order keeps a seed meaning the same thing here as there.
33 cr = randn(N, 1);
34 ci = randn(N, 1);
36 % The cube's integer wavenumbers, -m:m x -n:n x -p:p in column-major order.
37 i = (0:N-1).';
38 jx = mod(i, m2) - m;
39 jy = mod(floor(i/m2), n2) - n;
40 jz = floor(i/(m2*n2)) - p;
42 % Confine to a ball for isotropy.
43 keep = ((jx/m).^2 + (jy/n).^2 + (jz/p).^2) <= 1;
44 jx = jx(keep);
45 jy = jy(keep);
46 jz = jz(keep);
47 cr = cr(keep);
48 ci = ci(keep);
50 % Normalize so the variance is 1 at each point.
51 s = 1/sqrt(numel(cr));
52 cr = s*cr;
53 ci = s*ci;
55 % Angular wavenumbers on the enlarged period, which is a whole number of
56 % wavelengths on each side.
57 kx = 2*pi*jx/(m*lambda);
58 ky = 2*pi*jy/(n*lambda);
59 kz = 2*pi*jz/(p*lambda);
61 % Fold the box's origin into the phase, so evaluating is a plain sum over
62 % cos(k.x) and sin(k.x) with no offset left to carry.
63 ph = -(kx*dom(1) + ky*dom(3) + kz*dom(5));
64 k = [kx, ky, kz];
65 c = [cr.*cos(ph) - ci.*sin(ph), cr.*sin(ph) + ci.*cos(ph)];
66end