function [Phi, Vref, V, data] = prepare_ACT5_data(datadir, reffname, datafname)
% PREPARE_ACT5_DATA  Load and normalize ACT5 trigonometric-pattern data.
%
% The ACT5 files store current patterns in columns and voltage responses in
% rows. This routine retains the L - 1 independent patterns, averages the
% measured frames, chooses zero-mean electrode voltages, and divides each
% current/voltage pair by the norm of its applied current pattern.
%
% Inputs:
%   datadir  - Path to the directory containing the ACT5 data files.
%   reffname - Base filename of the homogeneous-reference data file,
%              without the .mat extension.
%   datafname - Base filename of the target data file, without the .mat
%              extension.
%
% Outputs:
%   Phi        - L-by-(L - 1) matrix of normalized trigonometric current
%              patterns. Each column has Euclidean norm one.
%   Vref     - L-by-(L - 1) matrix of averaged homogeneous-reference
%              voltages. Each column has zero mean and has been divided by
%              the norm of its corresponding applied current pattern.
%   V        - L-by-(L - 1) matrix of averaged target voltages, processed
%              and normalized in the same way as Vref.
%   data     - Structure loaded from the target .mat file. It contains the
%              original target measurements and experimental metadata; the
%              entry-point script uses its tank and electrode dimensions.
%
% Assumptions:
%   1. The target and reference files use the same even number L of
%      electrodes and the same trigonometric current patterns.
%   2. Current patterns are stored in columns of cur_pattern.
%   3. Voltage frames are stored in frame_voltage, with current-pattern
%      index in the first dimension, electrode index in the second, and
%      repeated measurement frames in the third.
%
% Author: Melody Alsaker, Gonzaga University
% Date Modified: July 2026

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

refdata = load(fullfile(datadir, [reffname, '.mat']));
data = load(fullfile(datadir, [datafname, '.mat']));

J0 = data.cur_pattern;
L = size(J0, 1);
numCP = L - 1;

assert(mod(L, 2) == 0, 'This code assumes an even number of electrodes.');
assert(size(J0, 2) >= numCP, ...
    'The target data do not contain L - 1 current patterns.');
assert(size(refdata.cur_pattern, 2) >= numCP, ...
    'The reference data do not contain L - 1 current patterns.');
assert(size(data.frame_voltage, 1) >= numCP && size(data.frame_voltage, 2) == L, ...
    'The target voltage array has dimensions incompatible with the current patterns.');
assert(size(refdata.frame_voltage, 1) >= numCP && size(refdata.frame_voltage, 2) == L, ...
    'The reference voltage array has dimensions incompatible with the current patterns.');
assert(size(data.frame_voltage, 1) == size(refdata.frame_voltage, 1) && ...
       size(data.frame_voltage, 2) == size(refdata.frame_voltage, 2), ...
    'The target and reference voltage arrays have incompatible dimensions.');

% The target and reference data must describe the same tank and electrodes.
target_geometry = [data.circumference, data.elec_height, data.elec_width];
ref_geometry = [refdata.circumference, refdata.elec_height, refdata.elec_width];
assert(all(isfinite(target_geometry)) && all(target_geometry > 0), ...
    'The target tank and electrode dimensions must be positive and finite.');
assert(norm(target_geometry - ref_geometry, inf) <= ...
       1e-12*max(1, norm(target_geometry, inf)), ...
    'The target and reference tank or electrode dimensions do not match.');

% The target and reference measurements must use the same current patterns.
refJ0 = refdata.cur_pattern(:, 1:numCP);
J0 = J0(:, 1:numCP);
assert(norm(J0 - refJ0, 'fro') <= 1e-12 * max(1, norm(J0, 'fro')), ...
    'The target and reference current patterns do not match.');

% Retain the real (in-phase) component of the measured phasor voltages,
% consistent with the purely resistive conductivity model. The discarded
% imaginary component is 90 degrees out of phase with the applied current
% and is small in the supplied datasets (about 1.15% in RMS magnitude).
Vref = real(mean(refdata.frame_voltage, 3));
V = real(mean(data.frame_voltage, 3));

% Retain the L - 1 independent patterns and transpose so that each column
% is one electrode-voltage vector.
Vref = Vref(1:numCP, :).';
V = V(1:numCP, :).';

% Electrode potentials are defined only up to an additive constant, establish choice of ground.
Vref = Vref - mean(Vref, 1);
V = V - mean(V, 1);

% Normalize each applied current and its voltage response together.
current_norms = vecnorm(J0, 2, 1);
assert(all(current_norms > 0), 'A selected current pattern has zero norm.');
Phi = J0 ./ current_norms;
Vref = Vref ./ current_norms;
V = V ./ current_norms;

end
