function [t, acc, vel, disp, accRaw] = sine_burst_acc(f0, nCycles, fs, amp, r)
%SINE_BURST_ACC
% Generate an acceleration sine burst with stable velocity and displacement.
%
% Velocity and displacement are corrected so they begin and end at zero.
%
%   [t, acc, vel, disp, accRaw] = sine_burst_acc(f0, nCycles, fs, amp, r)

    if nargin < 4
        amp = 1;
    end

    if nargin < 5
        r = 0.2;       % Tukey taper ratio
    end

    if f0 <= 0 || nCycles <= 0 || fs <= 0
        error('f0, nCycles, and fs must be positive.');
    end

    if r < 0 || r > 1
        error('Tukey ratio r must be between 0 and 1.');
    end

    dt = 1/fs;
    T  = nCycles/f0;

    N = round(T*fs) + 1;
    t = (0:N-1)'*dt;

    % ------------------------------------------------------------
    % Raw acceleration sine burst
    % ------------------------------------------------------------
    accRaw = amp*sin(2*pi*f0*t) .* tukeywin(N,r);

    % ------------------------------------------------------------
    % First integration: acceleration -> velocity
    % ------------------------------------------------------------
    velRaw = cumtrapz(t, accRaw);

    % Force velocity to start and end at zero
    velDrift = velRaw(1) + (velRaw(end)-velRaw(1))*(t/t(end));
    vel1 = velRaw - velDrift;

    % ------------------------------------------------------------
    % Second integration: velocity -> displacement
    % ------------------------------------------------------------
    dispRaw = cumtrapz(t, vel1);

    % Force displacement to start and end at zero
    phi = (t/t(end)).^2;
    disp = dispRaw - dispRaw(end)*phi;

    % Enforce exact endpoint values
    disp(1)   = 0;
    disp(end) = 0;

    % ------------------------------------------------------------
    % Recompute velocity and acceleration consistently
    % ------------------------------------------------------------
    vel = gradient(disp, dt);
    acc = gradient(vel, dt);

    % Enforce exact velocity endpoint values
    vel(1)   = 0;
    vel(end) = 0;

end