function FFTbeta = comp_FFTbeta(k)
% COMP_FFTBETA  Precompute the FFT representation of 1/(pi*k).
%
% This function constructs the Cauchy kernel used in the FFT-based
% convolution in the D-bar equation. The enlarged grid and four-point
% average regularize the singular kernel and prevent wraparound from
% contaminating the central Mk-by-Mk grid.
%
% Inputs:
%   k        - Mk-by-Mk complex scattering-frequency grid.
%
% Output:
%   FFTbeta  - Mtimes2-by-Mtimes2 complex array containing the two-
%              dimensional FFT of the regularized Cauchy kernel.
%
% Assumptions:
%   1. Mk is even, Mdiv2 = Mk/2, and Mtimes2 = 2*Mk.
%   2. The coordinate vector x is uniformly spaced with step size h and
%      spans an interval whose largest absolute value is s.
%
% Author: Melody Alsaker, Gonzaga University
% Date Modified: July 2026

% ==============================================================================================
% ==============================================================================================

Mk = size(k, 1);
Mdiv2 = Mk/2;
Mtimes2 = 2*Mk;

x = real(k(1, :));
h = x(2) - x(1);
s = max(abs(x));

xBig = [-(s + (Mdiv2:-1:1)*h), x, s + (1:Mdiv2)*h];
[K1Big, K2Big] = meshgrid(xBig, xBig);
k_Big = K1Big + 1i*K2Big;

% Evaluate the Cauchy kernel on a staggered grid. Because zero is not a
% grid point, no singular value needs to be removed.
beta = h^2 ./ (pi*k_Big);

% The staggered grid admits two possible centerings in each coordinate.
% The index vectors p and p1 implement these alternative shifts; averaging
% their four row-column combinations centers the Cauchy kernel symmetrically
% at the origin.
p  = circshift(1:Mtimes2, Mk);
p1 = circshift(1:Mtimes2, Mk + 1);

beta_avg = (beta(p, p) + beta(p1, p1) + beta(p1, p) + beta(p, p1))/4;

% The symmetric average is theoretically zero at the origin. Explicitly
% enforce this principal-value value against roundoff.
beta_avg(1, 1) = 0;

FFTbeta = fft2(beta_avg);

end
