% RUNTHIS_DBAR_TEXP_TRIGPATT_CIRCLE  Compute a t^exp D-bar reconstruction.
%
% This script provides a very simple interface for a D-bar reconstruction solver. It loads ACT5 
% target and homogeneous-reference measurements, constructs a difference DN matrix, computes the
% Born-type scattering transform texp, and reconstructs the conductivity on the unit disk.
%
% Assumptions:
%   1. The target is a circular tank.
%   2. The number of electrodes, L, is even. 
%   3. ACT5 applies trigonometric current patterns in the order
%          cos(theta), ..., cos(L*theta/2),
%          sin(theta), ..., sin((L/2 - 1)*theta).
%   4. The reference data were measured in homogeneous saline.
%
% Required files:
%   prepare_ACT5_data.m
%   comp_texp.m
%   solve_Dbar_eq.m
%   comp_FFTbeta.m
%
% Author: Melody Alsaker, Gonzaga University
% Date Modified: July 2026

% ==============================================================================================
% ==============================================================================================
clear;
close all;

use_parallel = true;  % Do you want to parallelize over the z-grid to solve the D-bar equation?

% Open the parallel pool before starting the timer for the main algorithm
if use_parallel
    pool = gcp;
    max_workers = pool.NumWorkers;
else
    max_workers = 0;
end

% This end-to-end timer includes parallel-pool startup if a new pool is opened.
timestart = tic;


% ==============================================================================================
% Select the target and reconstruction parameters

% This script is intended to sit beside the ACT5data directory.
datadir = fullfile(fileparts(mfilename('fullpath')), 'ACT5data');

reffname = 'saline_270_22_11_07_10_16_18';          % File name for homogeneous reference data
datafname = 'two_large_red_22_11_07_10_27_08';      % File name for target data

sigma_ref = 0.278;    % Measured homogeneous saline conductivity [S/m]

Mk = 32;              % Number of points per direction on the spectral k-grid
trunc = 5;            % Scattering-transform truncation radius

hz = 0.0442;          % Coarse solve grid (approximately 1600 points in the unit disk)
hz_plot = 1/128;      % Fine display grid

cmap = 'jet';         % Colormap used for plots


% ==============================================================================================
% Load and prepare the ACT5 data

% The subroutine averages the measured frames, chooses zero-mean electrode
% voltages, and normalizes each current pattern and voltage response together.
[Phi, Vref, V, data] = prepare_ACT5_data(datadir, reffname, datafname);

L = size(Phi, 1);              % Number of electrodes. 
assert(mod(L, 2) == 0, 'This code assumes an even number of electrodes.');

numCP = size(Phi, 2);          % Number of linearly independent current patterns
Ldiv2 = L/2;


% ==============================================================================================
% Construct the difference DN matrix

% These are discrete electrode transfer-impedance (ND) matrices [ohm] in
% the ACT5 trigonometric basis.
refR = Phi.' * Vref;
R = Phi.' * V;

% Invert the ND matrices to obtain the corresponding conductance (DN)
% matrices [S].
refLambda = refR \ eye(numCP);
Lambda = R \ eye(numCP);

dLambda = Lambda - refLambda;

% Convert the measured DN-matrix difference to the dimensionless unit-disk
% scaling used in the D-bar equations.
perim = data.circumference * 0.0254;        % Tank perimeter [m] (stored in inches in the dataset)
rtank = perim/(2*pi);                        % Tank radius [m]
dtheta = 2*pi/L;
eArea = data.elec_height * data.elec_width; % Area of each electrode [m^2]

scaling_factor = rtank*dtheta/(eArea*sigma_ref);
dLambda = scaling_factor*dLambda;


% ==============================================================================================
% Construct the spatial grid on the unit disk

Nz = floor(2/hz) + 1;
xx = ((1:Nz) - (Nz + 1)/2)*hz;
[z1, z2] = meshgrid(xx, xx);

zgrid = z1 + 1i*z2;
zidx = find(abs(zgrid) <= 1);
z = zgrid(zidx);


% ==============================================================================================
% Construct the scattering-frequency grid

assert(mod(Mk, 2) == 0, 'Mk must be even so that the k-grid excludes k = 0.');

x = linspace(-trunc, trunc, Mk);
[K1, K2] = meshgrid(x, x);
k = K1 + 1i*K2;


% ==============================================================================================
% Compute the Born-type scattering transform

texp = comp_texp(dLambda, Ldiv2, k, trunc);

% ==============================================================================================
% Solve the D-bar equation and recover the conductivity

tol = 1e-5;           % GMRES relative residual tolerance
maxit = 10;           % Maximum number of restart cycles
restrt = 5;           % Number of GMRES iterations per restart cycle

% For each z, solve the D-bar equation in the k variable and evaluate the
% resulting CGO solution at k = 0.
mu0 = solve_Dbar_eq(texp, k, z, tol, maxit, restrt, max_workers);

% The recovered conductivity should theoretically be real; discard the small
% imaginary residue introduced by discretization and the numerical D-bar solve.
sigmatemp = real(sigma_ref * mu0.^2);

% Done with actual reconstruction process! Display runtime
total_runtime = toc(timestart);
fprintf('Total runtime for reconstruction: %.1f seconds\n', total_runtime);


% ==============================================================================================
% Interpolate onto a finer grid for display.

xx_plot = -1:hz_plot:1;
[z1_plot, z2_plot] = meshgrid(xx_plot, xx_plot);
zgrid_plot = z1_plot + 1i*z2_plot;

F = scatteredInterpolant(real(z), imag(z), sigmatemp,'linear', 'nearest');

% Nearest-neighbor extrapolation fills only the thin region between the
% convex hull of the coarse points and the circular boundary. This display
% interpolation does not add reconstruction resolution.
sigma_plot = F(z1_plot, z2_plot);
sigma_plot(abs(zgrid_plot) > 1) = nan;

% ==============================================================================================
% Display the scattering transform and reconstruction

figure;
tiledlayout(1, 2, 'Padding', 'compact', 'TileSpacing', 'compact');

nexttile;
imagesc(x, x, real(texp));
set(gca, 'YDir', 'normal');
axis image;
colorbar;
title('Real part of t^{exp}');

nexttile;
imagesc(x, x, imag(texp));
set(gca, 'YDir', 'normal');
axis image;
colorbar;
title('Imaginary part of t^{exp}');

colormap(cmap);


figure;
imagesc(xx_plot, xx_plot, sigma_plot);
set(gca, 'YDir', 'normal');
axis image off;
colorbar;
title('D-bar reconstruction (S/m)');
colormap(cmap);


