function texp = comp_texp(dLambda, Ldiv2, k, trunc)
% COMP_TEXP  Compute the Born-type scattering transform on the k-grid.
%
% This function computes the t^exp approximation to the scattering
% transform directly from the difference DN matrix. On the unit circle,
%
%   texp(k) = integral exp(i*conj(k)*conj(z)) * dLambda[exp(i*k*z)] ds.
%
% Inputs:
%   dLambda - Dimensionless difference DN matrix in the normalized
%             trigonometric basis. This is an (L - 1)-by-(L - 1) matrix.
%   Ldiv2   - Half the number of electrodes, L/2.
%   k       - Mk-by-Mk complex scattering-frequency grid.
%   trunc   - Scattering-transform truncation radius. The transform is set
%             to zero at grid points for which abs(k) >= trunc.
%
% Output:
%   texp    - Mk-by-Mk complex array containing the Born-type scattering
%             transform evaluated on the k-grid.
%
% Assumptions:
%   1. The target is the unit disk.
%   2. The number of electrodes, L, is even. 
%   3. The rows and columns of dLambda are ordered as
%          cos(theta), ..., cos(L*theta/2),
%          sin(theta), ..., sin((L/2 - 1)*theta).
%   4. The trigonometric current patterns have been normalized.
%
% For each retained scattering frequency k, a_k and a_kbar contain the
% trigonometric-basis coefficients of the two boundary functions appearing
% in
%
%   texp(k) = integral exp(i*conj(k)*conj(zeta)) * dLambda[exp(i*k*zeta)] ds.
%
% On the unit circle, zeta = exp(i*theta), so
%
%   exp(i*k*zeta)             = sum_j (i*k)^j/j! * exp(i*j*theta),
%   exp(i*conj(k)*conj(zeta)) = sum_j (i*conj(k))^j/j! * exp(-i*j*theta).
%
% Thus a_k represents exp(i*k*zeta), which is acted on by dLambda, while
% a_kbar represents exp(i*conj(k)*conj(zeta)), which serves as the test
% function in the boundary pairing. The factors +i and -i combine the
% cosine and sine basis functions into exp(i*j*theta) and exp(-i*j*theta),
% respectively.
%
% The constant Fourier mode is omitted because dLambda acts on mean-free
% boundary data. The Nyquist cosine mode receives sqrt(2) because it has no
% independent sine partner.
%
% Author: Melody Alsaker, Gonzaga University
% Date Modified: July 2026

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

numCP = size(dLambda, 1);
Mk = size(k, 1);

assert(size(dLambda, 2) == numCP, 'dLambda must be square.');
assert(2*Ldiv2 - 1 == numCP, 'dLambda has the wrong size for the specified number of electrodes.');
assert(isequal(size(k), [Mk, Mk]), 'The k-grid must be square.');

numk = numel(k);
kidx = find(abs(k) < trunc);
ktrunc = k(kidx);

numktrunc = numel(ktrunc);
conjktrunc = conj(ktrunc);

a_k = zeros(numCP, numktrunc);
a_kbar = zeros(numCP, numktrunc);

for j = 1:Ldiv2 - 1
    ak_k = ((1i*ktrunc).^j)/factorial(j);
    akbar_j = ((1i*conjktrunc).^j)/factorial(j);

    a_k(j, :) = ak_k.';
    a_k(Ldiv2 + j, :) = (1i*ak_k).';

    a_kbar(j, :) = akbar_j.';
    a_kbar(Ldiv2 + j, :) = (-1i*akbar_j).';
end

ak_L2 = ((1i*ktrunc).^Ldiv2)/factorial(Ldiv2);
akbar_L2 = ((1i*conjktrunc).^Ldiv2)/factorial(Ldiv2);

a_k(Ldiv2, :) = (sqrt(2)*ak_L2).';
a_kbar(Ldiv2, :) = (sqrt(2)*akbar_L2).';

% Convert the raw sine/cosine coefficients to coefficients in the
% Euclidean-normalized ACT5 trigonometric basis.
basis_scale = sqrt(Ldiv2);
a_k = basis_scale*a_k;
a_kbar = basis_scale*a_kbar;

texp = zeros(numk, 1);

% The boundary pairing is bilinear, not a Hermitian inner product, so
% a_kbar is intentionally not conjugated in this expression.
texp(kidx) = sum(a_kbar .* (dLambda*a_k), 1).';
texp = reshape(texp, Mk, Mk);

end
