%
%   quake_srs_synth_phase_333.m  ver 1.7  by Tom Irvine
%
function[dB_error,raccel,ramp,last_wavelet,f,amp,NHS,td,...
           wavelet_low,wavelet_up,alpha,beta]=...
           quake_srs_synth_phase_333(f,dt,amp,NHS,td,t,last_wavelet,...
              srs_spec_12,dlimit,fn,raccel,a1,a2,b1,b2,b3,nsec,...
                   error_limit,wavelet_low,wavelet_up,alpha,beta)
%

num=length(fn);     

tpi=2*pi;

num2=length(t);

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

fl=3/duration;
fu=sr/10;

fu = min(fu, max(fn));

fh=(fu+fl)/2;

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;

delta_f=fu-fl;

fm1=delta_f*rand()+fl;
fm2=delta_f*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);
%
        out1=sprintf(' frequency case %d ',ie);
        disp(out1);
%
        for ij=1:2000
%            
            x1=A*(-0.5+1*rand());
            
            if(ij<10 || rand()<0.5)
                x2=delta_f*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*(delta_f*rand()+fl);                
            end
            
            
            iflag=0;
            while(1)
                for iv=1:100
%
                    if( pi*x3/x2 + x4 < duration )
                        iflag=1;
                        break;
                    else
                        iflag=2;
                        x2=x2*1.01;
                        x4=0.95*x4;
                    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)
                disp(out1);                    
                x2=tpi*(delta_f*rand()+fl);                 
            end    
                
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
            
            if(ij==1)
                x1=1e-10;
            end

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

            t1=x4 + t(1);
            t2=t1 + pi*x3/x2;
                
            ia=ceil(t1/dt);
            ib=ia+floor((t2-t1)/dt);
            
            ia = max(ia, 1);
            ib = min(ib, num2);
                       
            % for iq=ia:ib
            %
            %    arg=x2*(t(iq)-t1);  
            %    aa(iq)=aa(iq)+x1*sin(arg/double(x3))*sin(arg);    
            % end

            iq = ia:ib;
            arg = x2 * (t(iq) - t1);
            aa(iq) = aa(iq) + x1 * sin(arg / x3) .* sin(arg);

            error=0;            
            lerr=0;
            merr=0;
            ierr=0;


            [xmax, xmin] = ws_srs_recursion(aa, a1, a2, b1, b2, b3);
                               
            for j=1:num
%    
                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                
                
%%%                
                lerr = min([lerr, e1, e2]);               
%%%                

                if(abs(e1)>ierr)
                    ierr=e1a;
                end
                if(abs(e2)>ierr)
                    ierr=e2a;
                end                
                
                error=error+e1a+e2a;          
                
            end
            
%            error
%            lerr
%            merr
            
            error=error+num*(abs(lerr)+abs(merr));
            error_p=abs(lerr)+abs(merr);
            

            if( (error<error_max || error_p<error_peaks) && x2<omegau)
                
                x1r(ie)=x1;
                x2r(ie)=x2;
                x3r(ie)=x3;
                x4r(ie)=x4;
                iar=ia;
                ibr=ib;
                               
                if(error<error_max)
                    error_max=error;
                    out1=sprintf('%d %9.5g %6.4g %6.4g %6.3g %6.3g %6.3g %6.3g',...
                              ij,error_max,merr,abs(lerr),x1,x2/tpi,x3,x4);                    
                end
                if(error_p<error_peaks)
                    error_peaks=error_p;
                    out1=sprintf('%d %9.5g %6.4g %6.4g %6.3g %6.3g %6.3g %6.3g *',...
                               ij,error_max,merr,abs(lerr),x1,x2/tpi,x3,x4);
                end             
               
                disp(out1);
                
                if(merr<error_limit && abs(lerr)<error_limit)
                    break;
                end                
            end          
        end

        
        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;

        L=last_wavelet;

        nt=length(t);

        i=L;
        
        %
        ia=wavelet_low(i);
        ib=wavelet_up(i);   
        %
        tx=t(ia:ib)-td(i);
        sa=sin( alpha(i)*tx );        
        sb=sin(  beta(i)*tx );
        %

        accel=zeros(nt,1);
        accel(ia:ib)=accel(ia:ib)+(amp(i)*sa.*sb);
              
        raccel=raccel+accel;

        fprintf('**  %8.4g %8.4g %8.4g \n',merr,abs(lerr),error_limit);
        
        if(merr<error_limit && abs(lerr)<error_limit)
            progressbar(1);
            pause(0.1);
            break;
        end
end

dB_error=max([ merr,abs(lerr) ]);

progressbar(1)
pause(0.1);

ramp=amp;