0ae15cfSeed runs from smooth random fields, and add the blob geometryDan Fortunato 1% Smooth random function on the unit sphere — chebfun's randnfunsphere,
2% evaluated at the given (theta, phi) instead of returned as a spherefun.
3%
4% F = randnfunsphere(LAMBDA, THETA, PHI) is a combination of all spherical
5% harmonics up to degree floor(2*pi/LAMBDA) with independent N(0,1)
6% coefficients, normalized so the variance is 1 at each point.
7%
8% randnfunsphere(LAMBDA, THETA, PHI, 'monochromatic') uses only the
9% harmonics of that one degree, so every component has the same wave
10% number — chebfun's 'monochrome' option.
11%
12% Seed the draw with rng(...) before calling. This project has no chebfun
13% objects: what would be a spherefun there is returned here as values on the
14% grid the caller passes in.
16function f = randnfunsphere(lambda, theta, phi, type)
17 if ( nargin < 4 )
18 type = 'white';
19 end
20 % The unit sphere has circumference 2*pi, matching randnfun's deg = L/lambda.
21 deg = floor(2*pi/lambda);
22 if ( strncmpi(type, 'm', 1) )
23 c = randn(2*deg+1, 1);
24 c = sqrt(4*pi/numel(c)) * c; % normalize so the variance is 1
25 f = sphHarmSumFixedDeg(theta, phi, deg, c);
26 else
27 c = randn((deg+1)^2, 1);
28 c = sqrt(4*pi/numel(c)) * c; % normalize so the variance is 1
29 f = sphHarmSum(theta, phi, deg, c);
30 end
31end
33% All spherical harmonics up to degree deg, with coefficients ordered by
34% degree and order (0, -1,0,1, -2,-1,0,1,2, ...). Order +m carries
35% cos(m*phi), order -m carries sin(m*phi).
36function f = sphHarmSum(theta, phi, deg, c)
37 f = 1/sqrt(4*pi) * c(1) * ones(size(theta));
38 k = 1; % coefficients consumed so far
39 for l = 1:deg
40 cl = c(k+1 : k+2*l+1); % this degree's orders, -l..l
41 k = k + 2*l + 1;
42 f = f + sphHarmSumFixedDeg(theta, phi, l, cl);
43 end
44end
46% All spherical harmonics of the single degree l.
47function f = sphHarmSumFixedDeg(theta, phi, l, c)
48 m = (0:l).';
49 a = (-1).^m ./ sqrt((1 + double(m==0)) * pi);
50 costh = cos(theta(:)).'; % legendre wants cos(theta), in a row
51 G = legendre(l, costh, 'norm'); % (l+1) x npts
52 f = 0 * theta;
53 for mm = 0:l
54 f = f + a(mm+1) * c(l+1+mm) * (G(mm+1,:).' .* cos(mm*phi));
55 if mm > 0
56 f = f + a(mm+1) * c(l+1-mm) * (G(mm+1,:).' .* sin(mm*phi));
57 end
58 end
59end