1function G = cymbal2_grid(P)
2% CYMBAL2_GRID Polar grid, radial difference matrices, and the azimuthal
3% transform. The radial grid carries l = 1..Nr+2 lines: real nodes at
4% l = 3..Nr spanning [washer radius, rim radius], and two virtual lines
5% beyond each boundary, eliminated by the boundary-condition rows of the
6% assembled systems (Samejima 2021, Fig. 4). Azimuthal operators are
7% circulant, so every linear solve block-diagonalizes over the real
8% cos/sin basis; F maps physical theta samples to coefficients (U * F)
9% and Gi maps back (Uhat * Gi), with mcol giving each column's wavenumber.
11a = P.radius_m;
12c = P.washer_radius_m;
13Nr = P.n_radial;
14nth = P.n_theta;
15if mod(nth, 2) ~= 0
16 error('cymbal2: n_theta must be even');
17end
18if Nr < 8
19 error('cymbal2: n_radial must be at least 8');
20end
21nrt = Nr + 2;
22hr = (a - c) / (Nr - 3);
23if c <= hr
24 error('cymbal2: washer_radius too small for n_radial (need c > (a-c)/(n_radial-3))');
25end
27r = c + ((1:nrt)' - 3) * hr; % r(1) = c - 2hr may be tiny on coarse grids;
28r(1) = max(r(1), hr / 2); % row 1 of every operator is never used
30% Radial stencils (dense; nrt is small). Rows outside the valid range are
31% left zero and are overwritten by boundary-condition rows downstream.
32D1 = zeros(nrt, nrt);
33D2 = zeros(nrt, nrt);
34D3 = zeros(nrt, nrt);
35for l = 2:nrt - 1
36 D1(l, l - 1) = -1 / (2 * hr); D1(l, l + 1) = 1 / (2 * hr);
37 D2(l, l - 1) = 1 / hr^2; D2(l, l) = -2 / hr^2; D2(l, l + 1) = 1 / hr^2;
38end
39for l = 3:nrt - 2
40 D3(l, l - 2) = -1 / (2 * hr^3);
41 D3(l, l - 1) = 2 / (2 * hr^3);
42 D3(l, l + 1) = -2 / (2 * hr^3);
43 D3(l, l + 2) = 1 / (2 * hr^3);
44end
46% Azimuthal stencils as dense circulant matrices, applied on the right
47% (U * T2 etc., so they carry a transpose; T2 is symmetric, T1 is
48% antisymmetric and T1t below is arranged so Uth = U * T1t).
49hth = 2 * pi / nth;
50T2 = zeros(nth, nth);
51T1t = zeros(nth, nth);
52for j = 1:nth
53 jm = 1 + mod(j - 2, nth);
54 jp = 1 + mod(j, nth);
55 T2(j, j) = -2 / hth^2; T2(j, jm) = 1 / hth^2; T2(j, jp) = 1 / hth^2;
56 % row j of T1 has -1 at jm, +1 at jp (scaled); transpose: column j
57 T1t(jm, j) = -1 / (2 * hth);
58 T1t(jp, j) = 1 / (2 * hth);
59end
61% Real transform: columns ordered m = 0, (1c,1s), (2c,2s), ..., nth/2
62th = ((1:nth)' - 1) * hth;
63F = zeros(nth, nth);
64Gi = zeros(nth, nth);
65mcol = zeros(1, nth);
66F(:, 1) = 1 / nth;
67Gi(1, :) = 1;
68col = 2;
69for m = 1:nth / 2 - 1
70 F(:, col) = 2 * cos(m * th) / nth; Gi(col, :) = cos(m * th)'; mcol(col) = m; col = col + 1;
71 F(:, col) = 2 * sin(m * th) / nth; Gi(col, :) = sin(m * th)'; mcol(col) = m; col = col + 1;
72end
73F(:, col) = cos((nth / 2) * th) / nth;
74Gi(col, :) = cos((nth / 2) * th)';
75mcol(col) = nth / 2;
77% Eigenvalue of the circulant second-difference for wavenumber m
78lam2 = (2 * cos(mcol * hth) - 2) / hth^2;
80G.a = a; G.c = c; G.Nr = Nr; G.nrt = nrt; G.nth = nth;
81G.hr = hr; G.hth = hth;
82G.r = r; G.th = th;
83G.D1 = D1; G.D2 = D2; G.D3 = D3;
84G.T1t = T1t; G.T2 = T2;
85G.F = F; G.Gi = Gi; G.mcol = mcol; G.lam2 = lam2;
86G.Rd = 1.0 ./ r; % 1/r column, broadcast over theta
87G.R2d = 1.0 ./ r.^2;
88end