function results = jeffcott_main(iu, m, k, zeta, e, rpm_run)
% JEFFCOTT_MAIN
% Jeffcott rotor whirl using closed-form Laplace-transform solution.
%
% Inputs:
%   iu      = unit system flag
%               1 -> English: m in lbm, k in lbf/in, e in in
%               2 -> SI:      m in kg,  k in N/m,    e in m
%   m       = rotor disk mass
%   k       = lateral shaft stiffness
%   zeta    = damping ratio (dimensionless)
%   e       = mass eccentricity
%   rpm_run = rotor speed for Laplace time-response plot (RPM)
%             if omitted, defaults to 0.8 * critical speed

%% Unit handling
switch iu
    case 1
        m_si = m * 0.45359237;
        k_si = k * 175.126835;
        e_si = e * 0.0254;
        disp_unit  = 'in';
        disp_scale = 1/0.0254;
        force_lbl  = 'lbf-sec/in';
        c_scale    = 1 / 175.126835;
    case 2
        m_si = m;  k_si = k;  e_si = e;
        disp_unit  = 'mm';
        disp_scale = 1000;
        force_lbl  = 'N-sec/m';
        c_scale    = 1;
    otherwise
        error('jeffcott_main:iu','iu must be 1 (English) or 2 (SI)');
end

%% Rotor properties
wn = sqrt(k_si/m_si);
fn = wn/(2*pi);
rpm_crit = 60*fn;
c = 2*zeta*m_si*wn;

fprintf('Critical speed = %.2f Hz = %.0f RPM\n', fn, rpm_crit);
fprintf('Damping ratio  = %.4f\n', zeta);
fprintf('Equivalent c   = %.4g %s\n', c*c_scale, force_lbl);

results.iu       = iu;
results.wn       = wn;
results.fn       = fn;
results.rpm_crit = rpm_crit;
results.zeta     = zeta;
results.c        = c;

%% Default run speed
if nargin < 6 || isempty(rpm_run)
    rpm_run = 0.8 * rpm_crit;
end
Omega_run = rpm_run * 2*pi/60;

%% Frequency response
rpm   = linspace(100, 12000, 2000);
Omega = rpm * 2*pi/60;
r     = Omega / wn;
X     = e_si * r.^2 ./ sqrt((1-r.^2).^2 + (2*zeta*r).^2);
phi   = atan2(2*zeta*r, 1-r.^2);

results.rpm = rpm;
results.X   = X;
results.phi = phi;

% --- Figure 1: Phase (top) + Magnitude (bottom, 2x height) ---
figure;
t2 = tiledlayout(3, 1, 'TileSpacing','compact', 'Padding','compact');
ss = sprintf('Jeffcott Rotor Unbalance Response\nCritical speed = %.2f Hz = %.0f RPM', fn, rpm_crit);
title(t2, ss);

% Phase (top)
nexttile(1);
plot(rpm, rad2deg(phi), 'LineWidth', 1.8);
grid on;
xlabel('Rotor Speed (RPM)');
ylabel('Phase (deg)');
xline(rpm_crit, '--', '');
ylim([0 180]);
yticks(0:45:180);

% Magnitude (bottom, 2x height)
nexttile(2, [2 1]);
plot(rpm, X*disp_scale, 'LineWidth', 1.8);
grid on;
xlabel('Rotor Speed (RPM)');
ylabel(sprintf('Whirl Amplitude (%s)', disp_unit));
xl = xline(rpm_crit, '--', 'Critical Speed');
xl.LabelHorizontalAlignment = 'right';
xl.LabelVerticalAlignment   = 'bottom';

%% Orbit plots at 0.5x, 1x, 1.5x critical
rpm_cases = [0.5 1.0 1.5] * rpm_crit;
t = linspace(0, 0.1, 3000);

figure;
for i = 1:length(rpm_cases)
    Om = rpm_cases(i) * 2*pi/60;
    rr = Om / wn;
    A  = e_si * rr^2 / sqrt((1-rr^2)^2 + (2*zeta*rr)^2);
    ph = atan2(2*zeta*rr, 1-rr^2);
    x  = A * cos(Om*t - ph);
    y  = A * sin(Om*t - ph);

    subplot(1, 3, i);
    plot(x*disp_scale, y*disp_scale, 'LineWidth', 1.5);
    axis equal; grid on;
    xlabel(sprintf('x (%s)', disp_unit));
    ylabel(sprintf('y (%s)', disp_unit));
    title(sprintf('Orbit at %.0f RPM', rpm_cases(i)), ...
          'Units','normalized', 'Position',[0.5, 1.14, 0]);

    [lim, ticks_orb] = neat_axes(A*disp_scale);
    xlim([-lim lim]); ylim([-lim lim]);
    xticks(ticks_orb);
    yticks(ticks_orb);
end
% Scale figure 1.3x larger
pos = get(gcf, 'Position');
set(gcf, 'Position', [pos(1), pos(2), pos(3)*1.3, pos(4)*1.3]);

%% Laplace time response at rpm_run
t = linspace(0, 1, 5000);

[x, y] = jeffcott_laplace_solution(t, m_si, k_si, zeta, e_si, Omega_run);

figure;
plot(t, x*disp_scale); hold on;
plot(t, y*disp_scale);
grid on;
xlabel('Time (sec)');
ylabel(sprintf('Displacement (%s)', disp_unit));
title(sprintf('Laplace Solution Time Response at %.0f RPM', rpm_run));
legend('x', 'y');

[ylim_val, ytick_val] = neat_axes(max(abs([x*disp_scale, y*disp_scale])));
ylim([-ylim_val ylim_val]);
yticks(ytick_val);

%% Laplace transient orbit
figure;
plot(x*disp_scale, y*disp_scale, 'LineWidth', 1.1);
axis equal; grid on;
xlabel(sprintf('x (%s)', disp_unit));
ylabel(sprintf('y (%s)', disp_unit));
title(sprintf('Laplace Solution Transient Orbit at %.0f RPM', rpm_run));

[lim_tr, ticks_tr] = neat_axes(max(abs([x*disp_scale, y*disp_scale])));
xlim([-lim_tr lim_tr]); ylim([-lim_tr lim_tr]);
xticks(ticks_tr);
yticks(ticks_tr);

results.t       = t;
results.x       = x;
results.y       = y;
results.rpm_run = rpm_run;
end

%% ------------------------------------------------------------------------
function [x, y] = jeffcott_laplace_solution(t, m, k, zeta, e, Omega)
wn  = sqrt(k/m);
wd  = wn * sqrt(1 - zeta^2);
r   = Omega / wn;
A   = e * r^2 / sqrt((1-r^2)^2 + (2*zeta*r)^2);
phi = atan2(2*zeta*r, 1-r^2);

xss = A * cos(Omega*t - phi);
yss = A * sin(Omega*t - phi);

x0 = xss(1);
v0 = -A * Omega * sin(-phi);
y0 = yss(1);
w0 =  A * Omega * cos(-phi);

xh = exp(-zeta*wn*t) .* ( ...
     -x0 * cos(wd*t) + ...
     ((-v0 - zeta*wn*x0) / wd) .* sin(wd*t));
yh = exp(-zeta*wn*t) .* ( ...
     -y0 * cos(wd*t) + ...
     ((-w0 - zeta*wn*y0) / wd) .* sin(wd*t));

x = xss + xh;
y = yss + yh;
end

%% ------------------------------------------------------------------------
function [lim, ticks] = neat_axes(v)
% NEAT_AXES  Clean symmetric axis limit and tick positions for value v.
% Adapts ytick_linear_double_sided logic by Tom Irvine.

lim   = 0;
ticks = 0;

for i = 1:30
    n = i - 20;
    X = 10^n;

    if v <= 0.13*X
        lim   = 0.15*X;
        ticks = lim * [-1  -2/3  -1/3  0  1/3  2/3  1];
        return;
    end
    if v <= 0.17*X
        lim   = 0.2*X;
        ticks = lim * [-1  -0.5  0  0.5  1];
        return;
    end
    if v <= 0.26*X
        lim   = 0.3*X;
        ticks = lim * [-1  -2/3  -1/3  0  1/3  2/3  1];
        return;
    end
    if v <= 0.375*X
        lim   = 0.4*X;
        ticks = lim * [-1  -0.5  0  0.5  1];
        return;
    end
    if v <= 0.48*X
        lim   = 0.5*X;
        ticks = lim * [-1  -0.8  -0.6  -0.4  -0.2  0  0.2  0.4  0.6  0.8  1];
        return;
    end
    if v <= 0.96*X
        lim   = 1.0*X;
        ticks = lim * [-1  -0.5  0  0.5  1];
        return;
    end
end
end