function out = misalignment_timehistories(varargin)
% Synthesize realistic proximity probe relative-displacement time histories
% for a shaft with combined 1X (imbalance) and 2X (misalignment) content.
%
% Two orthogonal probes (Probe 1 = X axis, Probe 2 = Y axis, mounted 90 deg
% apart). EACH probe measures the full 1X + 2X motion along its own axis.
%
% Name/value parameters (defaults in parentheses):
%   rpm     (1800)      shaft speed, RPM            -> 1X = rpm/60 Hz
%   D1      (1.0)       1X displacement amplitude, mils 0-pk
%   D2      (0.7)       2X displacement amplitude, mils 0-pk
%   phi1    (0)         1X phase, rad (Probe 1 reference)
%   phi2    (-pi/6)     2X phase offset, rad
%   probe_sep (pi/2)    angular separation between probes, rad (90 deg typ.)
%   fs      (50000)     sample rate, Hz
%   T       (1)         record length, sec
%   verbose (true)      print summary to command window

    p = inputParser;
    p.addParameter('rpm',        1800);
    p.addParameter('D1',         1.0);
    p.addParameter('D2',         0.7);
    p.addParameter('phi1',       0);
    p.addParameter('phi2',      -pi/6);
    p.addParameter('probe_sep',  pi/2);
    p.addParameter('fs',         50000);
    p.addParameter('T',          1);
    p.addParameter('verbose',    true);
    p.parse(varargin{:});
    r = p.Results;

    f1 = r.rpm/60;        % 1X frequency, Hz
    f2 = 2*f1;            % 2X frequency, Hz
    w1 = 2*pi*f1;
    w2 = 2*pi*f2;

    dt = 1/r.fs;
    t  = (0:dt:r.T-dt).';

    % --- Probe 1 (X axis): full 1X + 2X content ---
    p1 = r.D1*sin(w1*t + r.phi1) ...
       + r.D2*sin(w2*t + r.phi1 + r.phi2);

    % --- Probe 2 (Y axis): same physical motion, sensed 90 deg around ---
    % Each harmonic is shifted by (harmonic order) x probe_sep, because a
    % geometric rotation of the sensing axis advances the n-th harmonic by
    % n times the angular separation.
    p2 = r.D1*sin(w1*t + r.phi1 + 1*r.probe_sep) ...
       + r.D2*sin(w2*t + r.phi1 + r.phi2 + 2*r.probe_sep);

    out.probe1    = [t, p1];   % Probe 1, full content
    out.probe2    = [t, p2];   % Probe 2, full content
    out.t         = t;
    out.f1        = f1;
    out.f2        = f2;
    out.disp_unit = 'mils';

    if r.verbose
        fprintf('\n--- Proximity Probe Synthesis ---\n');
        fprintf('  Shaft speed     : %.1f RPM\n', r.rpm);
        fprintf('  1X / 2X freq    : %.2f / %.2f Hz\n', f1, f2);
        fprintf('  D1 (1X) / D2(2X): %.3f / %.3f mils 0-pk\n', r.D1, r.D2);
        fprintf('  phi1 / phi2     : %.3f / %.3f rad\n', r.phi1, r.phi2);
        fprintf('  Probe sep       : %.1f deg\n', r.probe_sep*180/pi);
        fprintf('  fs / T          : %.0f Hz / %.3f s\n', r.fs, r.T);
        fprintf('  Both probes contain 1X + 2X content.\n');
    end
end