%   rainflow_bar3.m  ver 1.2  by Tom Irvine
function [Z, ampCenters, meanCenters, ampEdges, meanEdges] = rainflow_bar3(fig_num, YS, data, nAmpBins, nMeanBins, ampLimits, meanLimits)
% RAINFLOW_BAR3  Bin [amplitude, mean, cycles] into a 2D grid and plot with bar3.
%
% Inputs
%   fig_num     : figure number handle/id
%   YS          : label string for title (e.g., test name)
%   data        : Nx3 [amplitude, mean, cycles]
%   nAmpBins    : # bins along Amplitude (default 40)
%   nMeanBins   : # bins along Mean (default 40)
%   ampLimits   : [amin amax] (optional; default from data)
%   meanLimits  : [mmin mmax] (optional; default from data)
%
% Outputs
%   Z           : nMeanBins x nAmpBins matrix of summed cycles per bin
%   ampCenters  : 1 x nAmpBins
%   meanCenters : 1 x nMeanBins
%   ampEdges    : 1 x (nAmpBins+1)
%   meanEdges   : 1 x (nMeanBins+1)

    if nargin < 4 || isempty(nAmpBins),  nAmpBins  = 40; end
    if nargin < 5 || isempty(nMeanBins), nMeanBins = 40; end

    % Validate & sanitize
    if size(data,2) ~= 3
        error('data must be Nx3: [amplitude, mean, cycles].');
    end
    amp = data(:,1);
    mu  = data(:,2);
    cyc = data(:,3);

    good = isfinite(amp) & isfinite(mu) & isfinite(cyc);
    amp = amp(good);  mu = mu(good);  cyc = cyc(good);

    % ROBUST: handle empty after filtering
    if isempty(amp)
        ampEdges = linspace(0,1,nAmpBins+1);
        meanEdges = linspace(0,1,nMeanBins+1);
        ampCenters  = (ampEdges(1:end-1)  + ampEdges(2:end))  / 2;
        meanCenters = (meanEdges(1:end-1) + meanEdges(2:end)) / 2;
        Z = zeros(nMeanBins, nAmpBins);
        figure(fig_num); clf
        h = bar3(Z.'); %#ok<NASGU>
        title(sprintf('Rainflow Spectra  %s', YS));
        xlabel('Amplitude'); ylabel('Mean'); zlabel('Cycles');
        colormap(jet); colorbar; view(-100,56); grid on
        return
    end

    % Limits
    if nargin < 6 || isempty(ampLimits)
        amin = min(amp); amax = max(amp);
    else
        amin = ampLimits(1); amax = ampLimits(2);
    end
    % FIX: meanLimits is the 7th input
    if nargin < 7 || isempty(meanLimits)
        mmin = min(mu); mmax = max(mu);
    else
        mmin = meanLimits(1); mmax = meanLimits(2);
    end

    % Guard for degenerate ranges
    if amax == amin, amax = amin + eps(max(1,abs(amin))); end
    if mmax == mmin, mmax = mmin + eps(max(1,abs(mmin))); end

    % Edges & centers
    ampEdges    = linspace(amin, amax, nAmpBins+1);
    meanEdges   = linspace(mmin, mmax, nMeanBins+1);
    ampCenters  = (ampEdges(1:end-1)  + ampEdges(2:end))  / 2;
    meanCenters = (meanEdges(1:end-1) + meanEdges(2:end)) / 2;

    % -------- Bin weighted by cycles --------
    Z = zeros(nMeanBins, nAmpBins);
    try
        % Newer MATLABs: weights supported
        Z = histcounts2(mu, amp, meanEdges, ampEdges, 'Weights', cyc);
    catch
        % Fallback: discretize -> accumarray
        iM = discretize(mu,  meanEdges);   % row (mean) bin index
        iA = discretize(amp, ampEdges);    % col (amp) bin index
        mask = ~isnan(iM) & ~isnan(iA) & isfinite(cyc);
        subs = [iM(mask), iA(mask)];
        Z = accumarray(subs, cyc(mask), [nMeanBins, nAmpBins], @sum, 0);
    end

    % -------- Plot with bar3 --------
    figure(fig_num); clf
    h = bar3(Z.');   %#ok<NASGU> % transpose so x-group = Amplitude, y = Mean
    title(sprintf('Rainflow Spectra  %s', YS));
    xlabel('Amplitude');  ylabel('Mean');  zlabel('Cycles');
    colormap(jet); colorbar
    view(-100,56); grid on

    % Axis ticks at physical bin centers (thinned for readability)
    ax = gca;
    ax.XTick = 1:nAmpBins; ax.YTick = 1:nMeanBins;
    if nAmpBins > 24
        ax.XTick = unique([1, round(linspace(1,nAmpBins,8)), nAmpBins]);
    end
    if nMeanBins > 24
        ax.YTick = unique([1, round(linspace(1,nMeanBins,8)), nMeanBins]);
    end
    ax.XTickLabel = compose('%.2g', ampCenters(ax.XTick));
    ax.YTickLabel = compose('%.2g', meanCenters(ax.YTick));
    ax.TickLabelInterpreter = 'none';

    % Map face colors to heights
    kids = get(gca,'Children');
    for k = 1:numel(kids)
        if isa(kids(k),'matlab.graphics.chart.primitive.Bar')
            zdata = get(kids(k),'ZData');
            set(kids(k),'CData', zdata, 'FaceColor','interp');
        end
    end
end
