%
%    ws_gen_time.m  ver 2.6  (vectorized + BLAS + NHS odd≥3)  by Tom Irvine
%    Signature unchanged from your v2.3
%
%    [wavelet, th] = ws_gen_time(nhs, amp, omegaf, upper, nt, dt, td, ...
%                                igen, nv, wavelet, nspec)
%
function [wavelet, th] = ws_gen_time(nhs, amp, omegaf, upper, nt, dt, td, igen, nv, wavelet, nspec)

    % ---- Enforce valid half-sine counts: odd and >= 3 (applies to active individual) ----
    nhs_row = nhs(igen,:);                      % 1×nspec
    nhs_row(nhs_row < 3) = 3;
    evenMask = mod(nhs_row,2)==0;
    nhs_row(evenMask) = nhs_row(evenMask) + 1;
    nhs(igen,:) = nhs_row;                      % write-back

    % ---- Lightweight validation ----
    if ~isscalar(nt) || nt < 2 || ~isfinite(nt), error('ws_gen_time: nt must be >=2.'); end
    if ~isscalar(dt) || dt <= 0 || ~isfinite(dt), error('ws_gen_time: dt>0 required.'); end
    if ~isscalar(nspec) || nspec < 1 || ~isfinite(nspec), error('ws_gen_time: nspec>=1.'); end

    omegaf = omegaf(:); upper = upper(:);
    if numel(omegaf) ~= nspec || numel(upper) ~= nspec
        error('ws_gen_time: omegaf & upper must be length nspec.');
    end
    if any(~isfinite(omegaf) | omegaf<=0 | omegaf>1e20)
        error('ws_gen_time: invalid omegaf detected.');
    end
    if any(~isfinite(upper) | abs(upper)<1e-20)
        error('ws_gen_time: invalid upper window length detected.');
    end

    % ---- Build / reuse wavelets ----
    if nv == 1
        % Time (row)
        t = (0:nt-1)*dt;                         % 1×nt

        % Parameters (rows)
        td_row  = td(igen,:);                    % 1×nspec
        nhs_row = double(nhs(igen,:));           % 1×nspec (already odd≥3)

        % Local time τ per band (implicit expansion → nspec×nt)
        TT = t - td_row.';                       % nspec×nt

        % Active window mask: 0 ≤ τ ≤ upper_i
        M  = (TT >= 0) & (TT <= upper);          % nspec×nt

        % Phase argument ωτ
        TA = omegaf .* TT;                       % nspec×nt

        % Wavelet: sin(ωτ/NHS) * sin(ωτ) in active window
        W1 = sin( TA ./ nhs_row(:) );           % nspec×nt
        W2 = sin( TA );                         % nspec×nt
        W  = W1 .* W2;

        wavelet = zeros(nspec, nt);
        wavelet(M) = W(M);

        % (Optional) diagnostic
        % big = abs(wavelet) > 1e10;
        % if any(big(:)), [ii,jj]=find(big,1,'first'); fprintf('wavelet large @(%d,%d)=%.3g\n',ii,jj,wavelet(ii,jj)); end
    else
        % Reuse incoming wavelet
        if ~isequal(size(wavelet), [nspec, nt])
            error('ws_gen_time: when nv~=1, wavelet must be nspec×nt.');
        end
    end

    % ---- Assemble time history via BLAS (fast & stable) ----
    if size(amp,1) < igen || size(amp,2) < nspec
        error('ws_gen_time: amp must have at least (igen,1:nspec).');
    end
    th = amp(igen,1:nspec) * wavelet;           % 1×nt

    % ---- Final numeric sanity ----
    n2 = norm(th,2);
    if ~isfinite(n2) || n2 < 1e-20
        error('ws_gen_time: th invalid (NaN/Inf) or near-zero.');
    end
end
