%
%   quake_srs_synth_phase_33_flow.m   ver 1.9
%   Fully vectorized SRS evaluation (Smallwood recursion)
%   Dimensions-safe wavelet injection
%
function [raccel,ramp,last_wavelet,f,amp,NHS,td,...
          wavelet_low,wavelet_up,alpha,beta] = ...
          quake_srs_synth_phase_33_flow(f,dt,amp,NHS,td,t,last_wavelet,...
              srs_spec_12,dlimit,fn,raccel,a1,a2,b1,b2,b3,nsec,...
                   error_limit_low,error_limit_high,wavelet_low,wavelet_up,alpha,beta,flow)

num  = length(fn);
tpi  = 2*pi;
num2 = length(t);

duration = t(end) - t(1);
sr = 1/dt;

fl = 3/duration;
if fl < flow
    fl = flow;
end

fu = sr/10;
if fu > max(fn)
    fu = max(fn);
end

fh = (fu+fl)/2;
omegah = fh*tpi;
omegau = fu*tpi;

fprintf(' fl = %10.4g Hz  fh = %10.4g Hz  fu = %10.4g Hz \n\n', fl, fh, fu);

x1r = zeros(num2,1);
x2r = zeros(num2,1);
x3r = zeros(num2,1);
x4r = zeros(num2,1);

A = max(amp)*0.3;

fm1 = (fu-fl)*rand() + fl;
fm2 = (fu-fl)*rand() + fl;

error_max   = 1.0e+90;
error_peaks = 1.0e+90;

progressbar;

for ie = 1:nsec

    iar = 1;
    ibr = 2;

    ferr1 = fm1;
    ferr2 = fm2;

    progressbar(ie/nsec);
    fprintf(' frequency case %d \n', ie);

    for ij = 1:1200

        x1 = A*(-0.5 + rand());

        if (ij < 10 || rand() < 0.5)
            x2 = (fu-fl)*rand() + fl;
        end

        x3 = 3 + 2*round(50.*rand()*rand());
        x4 = duration*0.5*rand();

        x2 = x2*tpi;

        try
            if ij>20 && rand()>0.7
                x1 = x1r(ie)*(0.98 + 0.04*rand());
                x2 = x2r(ie)*(0.98 + 0.04*rand());
                x3 = x3r(ie);
                x4 = x4r(ie)*(0.98 + 0.04*rand());
            end
        catch
            disp('qw1');
            ij, ie, x1r(ie)
        end

        if rand()<0.85
            if rand()<0.5
                x2 = tpi*ferr1*(0.98+0.04*rand());
            else
                x2 = tpi*ferr2*(0.98+0.04*rand());
            end
        end

        while x2 > omegau
            x2 = tpi*((fu-fl)*rand() + fl);
        end
        if x2 < tpi*flow
            x2 = tpi*((fu-fl)*rand() + fl);
        end

        % --- ensure duration constraint ---
        iflag = 0;
        while 1
            for iv = 1:100
                if tpi*x3/(2*x2) + x4 < duration
                    iflag = 1;
                    break;
                else
                    iflag = 2;
                    x2 = x2*1.01;
                    x4 = x4*0.95;
                end
            end
            if iflag == 1
                break;
            else
                x3 = x3 - 2;
            end
        end

        if x3 < 3
            x2 = omegau;
            x3 = 3;
            x4 = 0;
        end

        if x2 > omegau
            x2 = tpi*((fu-fl)*rand() + fl);
        end

        if ij==1
            x1 = 1e-10;
        end
        if abs(x1)<1e-10
            x1 = 1e-10;
        end

        % ---------------------------------------------------------
        %                 WAVELET INJECTION
        % ---------------------------------------------------------
        aa = raccel;

        t1 = x4 + t(1);
        t2 = t1 + tpi*x3/(2*x2);

        ia = ceil(t1/dt);
        ib = ia + floor((t2-t1)/dt);

        if ia < 1, ia = 1; end
        if ib > num2, ib = num2; end

        iq  = ia:ib;
        arg = x2 * (t(iq) - t1);

        % safe vector multiply
        aa(iq) = aa(iq) + x1*( sin(arg/double(x3)) .* sin(arg) );

        % ---------------------------------------------------------
        %        ** NEW: BATCHED SRS (NO more filter() loop!) **
        % ---------------------------------------------------------
        %
        [xmax, xmin] = ws_srs_recursion(aa, a1, a2, b1, b2, b3);

        % ---------------------------------------------------------
        %          SAME ERROR METRICS AS ORIGINAL CODE
        % ---------------------------------------------------------
        error = 0;
        lerr  = 0;
        merr  = 0;
        ierr  = 0;

        for j = 1:num
            if fn(j) >= flow

                e1 = 20*log10(xmax(j)/srs_spec_12(j));
                e2 = 20*log10(xmin(j)/srs_spec_12(j));

                e1a = abs(e1);
                e2a = abs(e2);

                if e1a > merr
                    merr = e1a;
                    fm1  = fn(j);
                end
                if e2a > merr
                    merr = e2a;
                    fm2  = fn(j);
                end

                if e1 < lerr
                    lerr = e1;
                end
                if e2 < lerr
                    lerr = e2;
                end

                if e1a > ierr, ierr = e1a; end
                if e2a > ierr, ierr = e2a; end

                error = error + e1a + e2a;

            end
        end

        error   = error   + num*(abs(lerr)+abs(merr));
        error_p = abs(lerr)+abs(merr);

        if ( (error < error_max || error_p < error_peaks) && ...
             x2 < omegau && x2 >= tpi*flow)

            x1r(ie) = x1;
            x2r(ie) = x2;
            x3r(ie) = x3;
            x4r(ie) = x4;

            iar = ia;
            ibr = ib;

            if error < error_max
                error_max = error;
                fprintf('%d %9.5g %6.4g %6.4g %6.3g %6.3g %6.3g %6.3g\n',...
                       ij,error_max,merr,abs(lerr),x1,x2/tpi,x3,x4);
            end

            if error_p < error_peaks
                error_peaks = error_p;
                fprintf('%d %9.5g %6.4g %6.4g %6.3g %6.3g %6.3g %6.3g *\n',...
                        ij,error_max,merr,abs(lerr),x1,x2/tpi,x3,x4);
            end

            if merr < error_limit_high && lerr > error_limit_low
                break;
            end
        end
    end

    % ============================================================
    %            UPDATE WAVELET PARAMETERS (same as original)
    % ============================================================
    last_wavelet = last_wavelet + 1;

    f(last_wavelet)   = x2r(ie)/tpi;
    amp(last_wavelet) = x1r(ie);

    if x3r(ie)==0, x3r(ie)=3; end
    NHS(last_wavelet) = x3r(ie);

    td(last_wavelet) = x4r(ie);

    beta(last_wavelet)  = tpi*f(last_wavelet);
    alpha(last_wavelet) = beta(last_wavelet)/double(NHS(last_wavelet));

    wavelet_low(last_wavelet) = iar;
    wavelet_up(last_wavelet)  = ibr;

    % ---------------- FINAL ACCEL ADDITION -----------------
    L = last_wavelet;
    accel = zeros(num2,1);

    ia = wavelet_low(L);
    ib = wavelet_up(L);

    sa = sin(alpha(L)*(t(ia:ib)-td(L)));
    sb = sin(beta(L)*(t(ia:ib)-td(L)));

    % SAFE VERSION (prevents dimension mismatch)
    accel(ia:ib) = accel(ia:ib) + amp(L)*(sa(:).*sb(:));

    raccel = raccel + accel;

    fprintf('**  %8.4g %8.4g %8.4g \n', merr, abs(lerr), error_limit_high);

    if merr < error_limit_high && lerr > error_limit_low
        progressbar(1);
        pause(0.5);
        break;
    end

end

progressbar(1);
pause(0.1);

ramp = amp;

end
