function shaft_orbit_animation(caseName, nRev, fps)
% SHAFT_ORBIT_ANIMATION  Animate a rotating-shaft orbit measured by two
% proximity probes mounted 90 deg apart, with the time-domain X(t) and Y(t)
% waveforms displayed alongside the orbit (ADRE / DSA style layout).
%
%   Two eddy-current proximity probes observe the shaft surface 90 deg apart.
%   Probe X -> horizontal displacement, Probe Y -> vertical displacement.
%   Plotting Y(t) vs X(t) traces the shaft centerline path (the "orbit").
%   A once-per-rev keyphasor pulse provides the timing reference; counting
%   keyphasor marks on the orbit reveals the dominant order (1X, 2X, ...).
%
% USAGE
%   shaft_orbit_animation                      % default: 'unbalance'
%   shaft_orbit_animation('misalign')
%   shaft_orbit_animation('oilwhirl', 8, 30)   % case, # revs, frames/rev
%
% CASES
%   'balanced'  - light unbalance, near-circular 1X
%   'unbalance' - strong 1X, elliptical (anisotropic support stiffness)
%   'misalign'  - 1X + 2X, flattened / banana orbit
%   'figure8'   - 2X dominant, figure-eight (severe misalignment)
%   'oilwhirl'  - sub-synchronous ~0.45X, precessing loop
%   'rub'       - rotor rub, truncated/looped orbit with harmonics
%
% Tom Irvine style: direct, practical, self-contained. No toolboxes required.

    if nargin < 1 || isempty(caseName), caseName = 'unbalance'; end
    if nargin < 2 || isempty(nRev),     nRev     = 6;           end
    if nargin < 3 || isempty(fps),      fps      = 60;          end

    % ----------------------------------------------------------------
    % Signal model.  fr = rotating speed (Hz).  The probe signals are
    % built from synchronous (1X) and other order components per case.
    % ----------------------------------------------------------------
    fr  = 50;                 % rotating speed, Hz (3000 rpm)
    Trev = 1/fr;              % one revolution, s
    npr  = 200;               % samples per revolution for the smooth path
    [xfun, yfun, info] = get_case(caseName);

    % High-resolution time base for the full record (used for the static
    % waveform traces and the closed orbit path).
    tEnd = nRev*Trev;
    t    = linspace(0, tEnd, nRev*npr + 1);
    th   = 2*pi*fr*t;                 % shaft angle, rad
    X    = xfun(th);
    Y    = yfun(th);

    % Normalize to a common displacement scale (mils-like, arbitrary units).
    s = max([max(abs(X)), max(abs(Y)), eps]);
    X = X/s;  Y = Y/s;

    % Keyphasor fires once per rev at theta = 0, 2pi, 4pi, ...
    kpTh = (0:nRev)*2*pi;
    kpX  = xfun(kpTh)/s;
    kpY  = yfun(kpTh)/s;

    % ----------------------------------------------------------------
    % Figure / axes layout
    %   Left  : orbit (Y vs X)
    %   Right : X(t) on top, Y(t) on bottom
    % ----------------------------------------------------------------
    fig = figure('Color','w','Name','Shaft Orbit + Probe Waveforms',...
                 'Position',[100 100 1100 560]);

    % --- Orbit axis ---
    axO = axes('Parent',fig,'Position',[0.06 0.12 0.42 0.78]); hold(axO,'on');
    axis(axO,'equal');
    lim = 1.25;
    set(axO,'XLim',[-lim lim],'YLim',[-lim lim],'Box','on','FontSize',10);
    grid(axO,'on');
    plot(axO,[-lim lim],[0 0],'Color',[.7 .7 .7],'LineWidth',0.8);
    plot(axO,[0 0],[-lim lim],'Color',[.7 .7 .7],'LineWidth',0.8);
    % bearing clearance reference circle
    ang = linspace(0,2*pi,200);
    plot(axO,lim*0.95*cos(ang),lim*0.95*sin(ang),':','Color',[.8 .8 .8]);
    xlabel(axO,'Probe X  (horizontal), normalized');
    ylabel(axO,'Probe Y  (vertical), normalized');
    title(axO,sprintf('Shaft Orbit  -  %s', info.name),'FontWeight','bold');

    % full orbit path (faint, for context)
    plot(axO, X, Y, '-', 'Color',[0.62 0.78 0.95],'LineWidth',1.0);
    % keyphasor marks on orbit
    plot(axO, kpX, kpY, 'o','MarkerEdgeColor',[0.85 0.35 0.19],...
         'MarkerFaceColor','none','MarkerSize',9,'LineWidth',1.6);

    % animated handles
    hTrail = plot(axO, NaN, NaN, '-','Color',[0.10 0.37 0.65],'LineWidth',2.0);
    hDot   = plot(axO, NaN, NaN, 'o','MarkerFaceColor',[0.10 0.37 0.65],...
                  'MarkerEdgeColor',[0.10 0.37 0.65],'MarkerSize',9);

    % --- X(t) axis ---
    axX = axes('Parent',fig,'Position',[0.56 0.56 0.40 0.34]); hold(axX,'on');
    set(axX,'XLim',[0 tEnd*1e3],'YLim',[-lim lim],'Box','on','FontSize',10);
    grid(axX,'on');
    plot(axX, t*1e3, X, '-','Color',[0.62 0.78 0.95],'LineWidth',1.0);
    for k=1:numel(kpTh)
        tk = kpTh(k)/(2*pi*fr)*1e3;
        plot(axX,[tk tk],[-lim lim],':','Color',[0.85 0.35 0.19]);
    end
    ylabel(axX,'Probe X');
    title(axX,'Time-domain probe signals','FontWeight','bold');
    hXtrail = plot(axX, NaN, NaN, '-','Color',[0.10 0.37 0.65],'LineWidth',1.8);
    hXdot   = plot(axX, NaN, NaN, 'o','MarkerFaceColor',[0.10 0.37 0.65],...
                   'MarkerEdgeColor',[0.10 0.37 0.65],'MarkerSize',7);

    % --- Y(t) axis ---
    axY = axes('Parent',fig,'Position',[0.56 0.12 0.40 0.34]); hold(axY,'on');
    set(axY,'XLim',[0 tEnd*1e3],'YLim',[-lim lim],'Box','on','FontSize',10);
    grid(axY,'on');
    plot(axY, t*1e3, Y, '-','Color',[0.62 0.78 0.95],'LineWidth',1.0);
    for k=1:numel(kpTh)
        tk = kpTh(k)/(2*pi*fr)*1e3;
        plot(axY,[tk tk],[-lim lim],':','Color',[0.85 0.35 0.19]);
    end
    xlabel(axY,'Time, ms');
    ylabel(axY,'Probe Y');
    hYtrail = plot(axY, NaN, NaN, '-','Color',[0.10 0.37 0.65],'LineWidth',1.8);
    hYdot   = plot(axY, NaN, NaN, 'o','MarkerFaceColor',[0.10 0.37 0.65],...
                   'MarkerEdgeColor',[0.10 0.37 0.65],'MarkerSize',7);

    % annotation with diagnostic summary
    annotation(fig,'textbox',[0.06 0.90 0.42 0.07],'String',...
        sprintf('Dominant: %s   |   Shape: %s', info.freq, info.shape),...
        'EdgeColor','none','FontSize',10,'Color',[0.25 0.25 0.25]);

    % ----------------------------------------------------------------
    % Animation loop.  A sliding window of points forms the moving trail.
    % ----------------------------------------------------------------
    nPts = numel(t);
    win  = round(npr*1.0);            % trail length ~ one revolution
    dt   = 1/fps;

    for i = 1:nPts
        i0 = max(1, i-win);
        set(hTrail,'XData',X(i0:i),'YData',Y(i0:i));
        set(hDot,  'XData',X(i),   'YData',Y(i));

        set(hXtrail,'XData',t(i0:i)*1e3,'YData',X(i0:i));
        set(hXdot,  'XData',t(i)*1e3,   'YData',X(i));

        set(hYtrail,'XData',t(i0:i)*1e3,'YData',Y(i0:i));
        set(hYdot,  'XData',t(i)*1e3,   'YData',Y(i));

        drawnow limitrate;
        pause(dt);
        if ~ishghandle(fig), return; end   % allow closing mid-run
    end
end

% ====================================================================
function [xfun, yfun, info] = get_case(name)
% Return anonymous functions x(theta), y(theta) and a descriptor struct.
% theta is shaft angle in radians.  Amplitudes are arbitrary (normalized
% later).  Anisotropic support stiffness is modeled by unequal X/Y gains.

    switch lower(name)
        case 'balanced'
            xfun = @(th) cos(th);
            yfun = @(th) 0.92*sin(th);
            info = mk('Balanced (light unbalance)','1X','Near-circular');

        case 'unbalance'
            xfun = @(th) 1.00*cos(th);
            yfun = @(th) 0.55*sin(th);
            info = mk('Mass unbalance','1X','Ellipse');

        case 'misalign'
            xfun = @(th) cos(th);
            yfun = @(th) 0.60*sin(th) + 0.35*sin(2*th);
            info = mk('Misalignment','1X + 2X','Banana / flattened');

        case 'figure8'
            xfun = @(th) sin(2*th);
            yfun = @(th) 0.85*sin(th);
            info = mk('Severe misalignment','2X','Figure-eight');

        case 'oilwhirl'
            % sub-synchronous component at ~0.45X superimposed on 1X
            w = 0.45;
            xfun = @(th) cos(th) + 0.55*cos(w*th);
            yfun = @(th) sin(th) + 0.55*sin(w*th);
            info = mk('Oil whirl','~0.45X','Precessing loop');

        case 'rub'
            % harmonic-rich, partially truncated arc (contact flattening)
            xfun = @(th) min(0.78, 1+0.18*sin(3*th)).*cos(th);
            yfun = @(th) 0.80*(1+0.18*sin(3*th)).*sin(th);
            info = mk('Rotor rub','1X + harmonics','Truncated / looped');

        otherwise
            error('Unknown case "%s". Try: balanced, unbalance, misalign, figure8, oilwhirl, rub.', name);
    end
end

function s = mk(n,f,sh)
    s.name = n;  s.freq = f;  s.shape = sh;
end
