function cmp = plate_compare(A, B)
% PLATE_COMPARE  Compare two solutions of the same plate problem.
%
%   cmp = plate_compare(A, B) takes two result structs (for example from
%   plate_modal and plate_fd) and compares them over their common time window.
%   B is interpolated onto the time axis of A, so A should be the solution
%   with the coarser sampling in time if they differ.
%
%     cmp.t, cmp.w_a, cmp.w_b, cmp.vm_a, cmp.vm_b   aligned probe series
%     cmp.w_rel_l2, cmp.vm_rel_l2                    ||a - b|| / ||a||
%     cmp.w_peak_a, cmp.w_peak_b, cmp.vm_peak_a, cmp.vm_peak_b
%     cmp.f_peak_a, cmp.f_peak_b                     frequency of the largest spectral peak
%     cmp.frames_rel_l2                              per-frame relative field error of w, when
%                                                    the two results share a grid and frame times
%     cmp.frames_t_mismatch                          largest difference between the frame times [s]
%
%   Relative errors in stress are expected to exceed those in displacement,
%   since stress is a second derivative of the displacement.

T = min(A.t(end), B.t(end));
keep = A.t <= T;
t = A.t(keep);
w_a = A.probe.w(keep);
vm_a = A.probe.vm(keep);
w_b = interp1(B.t, B.probe.w, t);
vm_b = interp1(B.t, B.probe.vm, t);

cmp.t = t;
cmp.w_a = w_a; cmp.w_b = w_b;
cmp.vm_a = vm_a; cmp.vm_b = vm_b;
cmp.w_rel_l2 = norm(w_a - w_b) / max(norm(w_a), eps);
cmp.vm_rel_l2 = norm(vm_a - vm_b) / max(norm(vm_a), eps);
cmp.w_peak_a = max(abs(w_a)); cmp.w_peak_b = max(abs(w_b));
cmp.vm_peak_a = max(vm_a); cmp.vm_peak_b = max(vm_b);

fmax = min(A.spectrum.f(end), B.spectrum.f(end));
ia = A.spectrum.f <= fmax & A.spectrum.f > 0; ib = B.spectrum.f <= fmax & B.spectrum.f > 0;   % skip the DC bin
[~, ka] = max(A.spectrum.w(ia)); fa = A.spectrum.f(ia); cmp.f_peak_a = fa(ka);
[~, kb] = max(B.spectrum.w(ib)); fb = B.spectrum.f(ib); cmp.f_peak_b = fb(kb);

% Frames are compared directly when the grids match and the frame times agree to
% within half a sample of each solver (each snaps its frames to its own time
% samples, so they cannot coincide exactly).
tol = 0.5 * ((A.t(2) - A.t(1)) + (B.t(2) - B.t(1))) + 1e-12;
same_grid = isequal(size(A.frames.w), size(B.frames.w)) && ...
    max(abs(A.grid.x - B.grid.x)) < 1e-9 && max(abs(A.grid.y - B.grid.y)) < 1e-9 && ...
    max(abs(A.frames.t - B.frames.t)) <= tol;
cmp.frames_t_mismatch = max(abs(A.frames.t(1:min(end, numel(B.frames.t))) - B.frames.t(1:min(end, numel(A.frames.t)))));
if same_grid
    nf = numel(A.frames.t);
    cmp.frames_t = A.frames.t;
    cmp.frames_rel_l2 = zeros(1, nf);
    for k = 1:nf
        da = A.frames.w(:, :, k); db = B.frames.w(:, :, k);
        fin = isfinite(da) & isfinite(db);       % nodes off an elliptical plate are NaN
        cmp.frames_rel_l2(k) = norm(da(fin) - db(fin)) / max(norm(da(fin)), eps);
    end
    cmp.frames_diff = A.frames.w - B.frames.w;
else
    cmp.frames_t = [];
    cmp.frames_rel_l2 = [];
    cmp.frames_diff = [];
end
end