function G = cymbal2_grid(P)
% CYMBAL2_GRID Polar grid, radial difference matrices, and the azimuthal
% transform. The radial grid carries l = 1..Nr+2 lines: real nodes at
% l = 3..Nr spanning [washer radius, rim radius], and two virtual lines
% beyond each boundary, eliminated by the boundary-condition rows of the
% assembled systems (Samejima 2021, Fig. 4). Azimuthal operators are
% circulant, so every linear solve block-diagonalizes over the real
% cos/sin basis; F maps physical theta samples to coefficients (U * F)
% and Gi maps back (Uhat * Gi), with mcol giving each column's wavenumber.

a = P.radius_m;
c = P.washer_radius_m;
Nr = P.n_radial;
nth = P.n_theta;
if mod(nth, 2) ~= 0
  error('cymbal2: n_theta must be even');
end
if Nr < 8
  error('cymbal2: n_radial must be at least 8');
end
nrt = Nr + 2;
hr = (a - c) / (Nr - 3);
if c <= hr
  error('cymbal2: washer_radius too small for n_radial (need c > (a-c)/(n_radial-3))');
end

r = c + ((1:nrt)' - 3) * hr;   % r(1) = c - 2hr may be tiny on coarse grids;
r(1) = max(r(1), hr / 2);      % row 1 of every operator is never used

% Radial stencils (dense; nrt is small). Rows outside the valid range are
% left zero and are overwritten by boundary-condition rows downstream.
D1 = zeros(nrt, nrt);
D2 = zeros(nrt, nrt);
D3 = zeros(nrt, nrt);
for l = 2:nrt - 1
  D1(l, l - 1) = -1 / (2 * hr); D1(l, l + 1) = 1 / (2 * hr);
  D2(l, l - 1) = 1 / hr^2; D2(l, l) = -2 / hr^2; D2(l, l + 1) = 1 / hr^2;
end
for l = 3:nrt - 2
  D3(l, l - 2) = -1 / (2 * hr^3);
  D3(l, l - 1) = 2 / (2 * hr^3);
  D3(l, l + 1) = -2 / (2 * hr^3);
  D3(l, l + 2) = 1 / (2 * hr^3);
end

% Azimuthal stencils as dense circulant matrices, applied on the right
% (U * T2 etc., so they carry a transpose; T2 is symmetric, T1 is
% antisymmetric and T1t below is arranged so Uth = U * T1t).
hth = 2 * pi / nth;
T2 = zeros(nth, nth);
T1t = zeros(nth, nth);
for j = 1:nth
  jm = 1 + mod(j - 2, nth);
  jp = 1 + mod(j, nth);
  T2(j, j) = -2 / hth^2; T2(j, jm) = 1 / hth^2; T2(j, jp) = 1 / hth^2;
  % row j of T1 has -1 at jm, +1 at jp (scaled); transpose: column j
  T1t(jm, j) = -1 / (2 * hth);
  T1t(jp, j) = 1 / (2 * hth);
end

% Real transform: columns ordered m = 0, (1c,1s), (2c,2s), ..., nth/2
th = ((1:nth)' - 1) * hth;
F = zeros(nth, nth);
Gi = zeros(nth, nth);
mcol = zeros(1, nth);
F(:, 1) = 1 / nth;
Gi(1, :) = 1;
col = 2;
for m = 1:nth / 2 - 1
  F(:, col) = 2 * cos(m * th) / nth; Gi(col, :) = cos(m * th)'; mcol(col) = m; col = col + 1;
  F(:, col) = 2 * sin(m * th) / nth; Gi(col, :) = sin(m * th)'; mcol(col) = m; col = col + 1;
end
F(:, col) = cos((nth / 2) * th) / nth;
Gi(col, :) = cos((nth / 2) * th)';
mcol(col) = nth / 2;

% Eigenvalue of the circulant second-difference for wavenumber m
lam2 = (2 * cos(mcol * hth) - 2) / hth^2;

G.a = a; G.c = c; G.Nr = Nr; G.nrt = nrt; G.nth = nth;
G.hr = hr; G.hth = hth;
G.r = r; G.th = th;
G.D1 = D1; G.D2 = D2; G.D3 = D3;
G.T1t = T1t; G.T2 = T2;
G.F = F; G.Gi = Gi; G.mcol = mcol; G.lam2 = lam2;
G.Rd = 1.0 ./ r;               % 1/r column, broadcast over theta
G.R2d = 1.0 ./ r.^2;
end