function [ramp,error_i,wavelet_low,wavelet_up,alpha,beta,t] = ...
    quake_srs_synth_phase_1_flow_fast(f,dt,damp,amp,NHS,td,last_wavelet,iter, ...
        fn,srs_spec,dlimit,error_limit,iu,uewf,dur,flow,tdel)

tpi = 2*pi;
nt  = round(dur/dt);
t   = (0:nt)*dt;

% --- precompute ---
alpha = (tpi*f(:))./double(NHS(:));
beta  = tpi*f(:);
ud    = NHS(:)./(2*f(:));
wavelet_low = max(1,floor(td(:)/dt));
wavelet_up  = min(nt, wavelet_low + floor(ud/dt));

[a1,a2,b1,b2,b3,~,~,~,~,~] = srs_coefficients(fn,damp,dt);

amp(f > max(fn)) = 0;
ramp  = amp;
error_max = inf;
drec = inf;
iflag = 0;
progressbar;

% --- preallocate once ---
nfn = numel(fn);
xmax = zeros(nfn,1); xmin = zeros(nfn,1);

for i = 1:iter
    progressbar(i/iter);

    if iflag == 0
        accel = generate_time_history_wavelet_table_q_fast(f,amp,NHS,td,t,wavelet_low,wavelet_up,alpha,beta);
    end

    % --- vectorized SRS filter loop ---

    try
        [~, xmax, xmin] = srs_iir_peaks_mex(accel, a1,a2,b1,b2,b3);
    catch    
        % --- compute SRS filter responses efficiently ---
        xmax(:) = 0; xmin(:) = 0;  % reuse preallocated vectors
        for j = 1:nfn
            % Filter one band only (this is unavoidable; coefficients differ)
            resp = filter([b1(j) b2(j) b3(j)], [1 -a1(j) -a2(j)], accel);
            xmax(j) = max(resp);
            xmin(j) = abs(min(resp));
        end
    end





    % --- vectorized error calculation ---
    valid = f(:) >= flow & srs_spec(:) > 0;
    xratio = srs_spec(valid)./max(xmin(valid),1e-16);
    yratio = srs_spec(valid)./max(xmax(valid),1e-16);
    errA = abs(20*log10(xratio));
    errB = abs(20*log10(yratio));
    errA(xmin(valid)>srs_spec(valid)) = errA(xmin(valid)>srs_spec(valid))*uewf;
    errB(xmax(valid)>srs_spec(valid)) = errB(xmax(valid)>srs_spec(valid))*uewf;
    err_comb = max([errA;errB]);
    error_i = max(err_comb);
    e2 = sum(errA+errB);
    error_val = error_i*e2;

    % --- displacement check ---
    v = cumtrapz(t,accel);
    d = cumtrapz(t,v);
    dratio = max(abs(d));
    if iu==1, dratio = dratio*386; else, dratio = dratio*9.81*1000; end

    % --- accept improvement ---
    if (error_val<error_max) && (dratio<dlimit || dratio<drec)
        drec = min(drec,dratio);
        error_max = error_val;
        ramp = amp;
        ramp_srs = (xmax+xmin)/2;
        raccel = accel;
        fprintf('%d  %8.4g  %8.4g  %8.4g\n',i,error_val,error_i,dratio);
        if error_i < error_limit, break; end
    end

    % --- adaptive mutation ---
    if i < iter
        jf = find(f>=flow,1);
        if i>max(20,round(0.03*iter))
            iflag=1;
            jk = randi(last_wavelet);
            amp = ramp;
            if rand()<0.5 && f(jk)>=flow
                [~,ii] = min(abs(fn-f(jk)));
                sc = (srs_spec(ii)/max(ramp_srs(ii),1e-16))^0.6;
                amp(jk) = ramp(jk)*sc;
            else
                amp(jk) = ramp(jk)*(0.95+0.1*rand());
            end
            if f(jk)>flow && td(jk)>tdel
                delta = generate_time_history_wavelet_table_q_fast(f(jk),amp(jk)-ramp(jk),NHS(jk),td(jk),t,...
                    wavelet_low(jk),wavelet_up(jk),alpha(jk),beta(jk));
            else
                delta = 0;
            end
            accel = raccel + delta;
        else
            iflag=0;
            scale = 0.5+rand();
            amp(jf+1:end)=ramp(jf+1:end)*scale;
        end
    end
end
progressbar(1);
pause(0.05);
end
