%  ws_th_from_wavelet_table.m  ver 1.3  (allocation-lean + NHS odd≥3)
%  Based on your v1.1 signature; uses all inputs.
%
%  [acceleration, velocity, displacement] = ws_th_from_wavelet_table( ...
%      iwin, store_amp, store_NHS, store_td, dur, dt, nt, iunit, f)
%
function [acceleration, velocity, displacement] = ws_th_from_wavelet_table( ...
    iwin, store_amp, store_NHS, store_td, dur, dt, nt, iunit, f)

    % ---------- Units scale ----------
    isu = 1.0;
    if iunit == 1
        isu = 386.0;        % in/s^2 per g
    elseif iunit == 2
        isu = 9.81;         % m/s^2 per g
    end

    % ---------- Pull active row ----------
    amp = store_amp(iwin,:).';
    NHS = double(store_NHS(iwin,:)).';
    td  = store_td(iwin,:).';
    f   = f(:);
    last_wavelet = length(f);

    % ---------- Enforce NHS rule (odd and >=3) ----------
    NHS(NHS < 3) = 3;
    evenMask = mod(NHS,2)==0;
    NHS(evenMask) = NHS(evenMask) + 1;

    % ---------- Time vector ----------
    t = (0:nt-1).' * dt;

    % ---------- Precompute params ----------
    beta  = 2*pi*f(1:last_wavelet);
    alpha = beta ./ NHS(1:last_wavelet);
    upper = td(1:last_wavelet) + (NHS(1:last_wavelet) ./ (2 .* f(1:last_wavelet)));

    APB   = alpha + beta;
    AMB   = alpha - beta;

    den_guard = 1e-15;
    APB_abs = max(abs(APB), den_guard);
    AMB_abs = max(abs(AMB), den_guard);
    APB2    = APB_abs.^2;
    AMB2    = AMB_abs.^2;

    % ---------- Index windows ----------
    wavelet_low = zeros(last_wavelet,1);
    wavelet_up  = zeros(last_wavelet,1);
    for i = 1:last_wavelet
        ia = round( 0.5 +  (td(i)/dur)*nt );
        ib = round(-0.5 + (upper(i)/dur)*nt );
        if ia < 1, ia = 1; end
        if ib > nt, ib = nt; end
        if ib < ia, ib = ia; end
        wavelet_low(i) = ia;
        wavelet_up(i)  = ib;
    end

    % ---------- Outputs ----------
    accel = zeros(nt,1);
    velox = zeros(nt,1);
    dis   = zeros(nt,1);

    % ---------- Acceleration accumulation ----------
    for i = 1:last_wavelet
        ia = wavelet_low(i); ib = wavelet_up(i);
        if ia > ib, continue; end
        tt = t(ia:ib) - td(i);
        s1 = sin(alpha(i) * tt);
        s2 = sin(beta(i)  * tt);
        accel(ia:ib) = accel(ia:ib) + amp(i) * (s1 .* s2);
    end

    % ---------- Velocity (analytic integral) ----------
    for i = 1:last_wavelet
        ia = wavelet_low(i); ib = wavelet_up(i);
        if ia > ib, continue; end
        tt = t(ia:ib) - td(i);
        v_add = 0.5 * amp(i) * ( -sin(APB(i)*tt)./APB_abs(i) + sin(AMB(i)*tt)./AMB_abs(i) );
        velox(ia:ib) = velox(ia:ib) + v_add;
    end

    % ---------- Displacement (analytic integral) ----------
    for i = 1:last_wavelet
        ia = wavelet_low(i); ib = wavelet_up(i);
        if ia > ib, continue; end
        tt = t(ia:ib) - td(i);
        x_add = 0.5 * amp(i) * ( (1 - cos(APB(i)*tt))./APB2(i) - (1 - cos(AMB(i)*tt))./AMB2(i) );
        dis(ia:ib) = dis(ia:ib) + x_add;
    end

    % ---------- Units (as in your v1.1) ----------
    velox = velox * isu;
    dis   = dis   * isu;

    % ---------- Package outputs ----------
    acceleration = [t, accel];
    velocity     = [t, velox];
    displacement = [t, dis];
end
