function mu0 = solve_Dbar_eq(texp, k, z, tol, maxit, restrt, max_workers)
% SOLVE_DBAR_EQ  Solve the discretized D-bar equation and return mu(z,0).
%
% The D-bar equation is real-linear, rather than complex-linear, because
% its integral operator contains conjugation. MATLAB's GMRES solves the
% resulting real system using a matrix-free FFT implementation of the
% Cauchy convolution.
%
% Inputs:
%   texp   - Mk-by-Mk complex array containing the Born-type scattering
%            transform evaluated on the k-grid.
%   k      - Mk-by-Mk complex scattering-frequency grid.
%   z      - Vector of complex spatial points at which to reconstruct the
%            conductivity.
%   tol    - Relative residual tolerance used by GMRES.
%   maxit  - Maximum number of GMRES restart cycles at each spatial point.
%   restrt - Number of GMRES iterations allowed in each restart cycle.
%   use_parallel - Logical flag. If true, open or reuse a parallel pool. If
%                  false, execute the spatial loop serially.
%
% Output:
%   mu0    - Complex column vector containing mu(z,0) at each point in z,
%            in the same linear order as the input vector z.
%
% Assumptions:
%   1. The k-grid is a square, uniformly spaced Cartesian grid with an even
%      number Mk of points per direction.
%   2. Because this even grid excludes k = 0, mu(z,0) is approximated by
%      averaging the solution at the four central k-grid points.
%
% Author: Melody Alsaker, Gonzaga University
% Date Modified: July 2026

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

[Mk1, ~] = size(k);
assert(mod(Mk1, 2) == 0, 'The k-grid must have an even number of points.');

Mk = Mk1;
numk = numel(k);
numz = numel(z);

% Precompute the FFT representation of the Cauchy kernel 1/(pi*k).
fft_beta = comp_FFTbeta(k);

rhs = [ones(numk, 1); zeros(numk, 1)];
mu0 = zeros(numz, 1);
gmres_flags = zeros(numz, 1);

conjk = conj(k);

% The even k-grid excludes zero. Bilinear interpolation at k = 0 is the
% average of the four central grid values.
mid = Mk/2;
k0idx = [sub2ind([Mk, Mk], mid, mid), ...
         sub2ind([Mk, Mk], mid + 1, mid), ...
         sub2ind([Mk, Mk], mid, mid + 1), ...
         sub2ind([Mk, Mk], mid + 1, mid + 1)];


Tbase = texp ./ (4*pi*conjk);  

% The phase exponential is separable in the two k coordinates.
kx = real(k(1, :));
ky = imag(k(:, 1));


parfor (ii = 1:numz, max_workers)
    phase = exp(2i*ky*imag(z(ii))) * exp(-2i*kx*real(z(ii)));
    T = Tbase .* phase;

    [m, gmres_flags(ii)] = gmres( @(f) DBO(f, T, fft_beta), rhs, restrt, tol, maxit, [], [], rhs);

    mu = m(1:numk) + 1i*m(numk + 1:end);
    mu0(ii) = mean(mu(k0idx));
end

num_failed = nnz(gmres_flags);
if num_failed > 0
    warning('solve_Dbar_eq:GMRESConvergence', ...
        'GMRES did not fully converge at %d of %d spatial points.', ...
        num_failed, numz);
end

end


function result = DBO(f, T, fft_beta)
% DBO  Apply the real-linear D-bar operator using an FFT convolution.
%
% This local helper applies the discretized operator
%
%   f --> f - C_k(T*conj(f)),
%
% where C_k denotes convolution with the Cauchy kernel. Complex vectors are
% packed into real form so that MATLAB's GMRES can solve this real-linear
% system.
%
% Inputs:
%   f         - Length-2*numk real vector containing the real and imaginary
%               parts of the current GMRES iterate, packed as
%               [real(f_complex); imag(f_complex)].
%   T         - Mk-by-Mk complex array containing the z-dependent
%               coefficient in the D-bar integral operator.
%   fft_beta  - Mtimes2-by-Mtimes2 FFT representation of the regularized
%               Cauchy kernel, produced by comp_FFTbeta.
%
% Output:
%   result    - Length-2*numk real vector containing the packed real and
%               imaginary parts of the operator applied to f.
%
% This is a local function called by solve_Dbar_eq and is not intended to be
% called directly from the entry-point script.
%
% Author: Melody Alsaker, Gonzaga University
% Date Modified: July 2026

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

Mk = size(T, 1);
numk = numel(T);
center = (Mk/2 + 1):(3*Mk/2);

f = f(1:numk) + 1i*f(numk + 1:end);
temp_fT = conj(reshape(f, Mk, Mk)) .* T;

temp_fT_Big = zeros(size(fft_beta), 'like', temp_fT);
temp_fT_Big(center, center) = temp_fT;

temp_fT_Big = ifft2(fft_beta .* fft2(temp_fT_Big));
tmp = temp_fT_Big(center, center);

f = f - tmp(:);
result = [real(f); imag(f)];

end
