%
%    sfa_engine.m  ver 1.9   April 26, 2013
%
function[a,x1r,x2r,x3r,error_rth,syna,cnew] = ...
              sfa_engine(dur,a,t,~,~,flow,fup,ie,nt,x1r,x2r,x3r,...
                                            original,running_sum,istr)
%
clear ta;
clear error_rth
clear syna;
%
tp=2*pi;
%
yr=zeros(nt,1);
%
zcf=0;
%
if(istr==2)
    amp=a-mean(a);
    nzc=0;
    pzc=0;
    for i=1:length(a)-1
%
       if(amp(i) <=0 && amp(i+1) > 0)
           pzc=pzc+1;
       end
%       
       if(amp(i) >=0 && amp(i+1) < 0)
           nzc=nzc+1;      
       end 
    end
%
    zcf=((pzc+nzc)/2)/dur;  
%
    out1=sprintf(' zcf = %8.4g Hz ',zcf);
    disp(out1);
%
end
%
sd=std(a);
%
am=sd;
%
errormax=1.0e+53;
%
disp(' ');
disp('  Trial     Error      Amplitude   Freq(Hz)   Phase(rad)  ');
%
jk=0;
%
for j=1:nt
%
	jk=jk+1;
 %           
    if(jk==10000)
		out4=sprintf('\n %ld ',j);
        disp(out4);
		jk=0;
    end
%
	x1=rand;
	x2=rand;
	x3=rand;
%      
	x1= am*x1^2;
	x2=((fup-flow)*x2+flow)*tp;
	x3=x3*tp;
%
    if(istr==2)
       if(rand<0.8)
          x2=zcf*tp*(0.98+0.04*rand);
       end
    end
%
    if(rand<0.5 && j>fix(nt/10))
%      
		x1=x1r(ie);
		x2=x2r(ie);
		x3=tp*rand;
%    
    end
%
    if(rand<0.4 && j>fix(nt/10))
%      
		x1=x1r(ie);
		x2=x2r(ie);
		x3=x3r(ie)+tp*(0.05-0.1*rand);
%    
    end
%
    if(rand<0.3 && j>fix(nt/10))
%      
		x1=x1r(ie)*(0.98+0.04*rand);
		x2=x2r(ie);
		x3=x3r(ie);
%    
    end
    if(rand<0.2 && j>fix(nt/10))
%      
		x1=x1r(ie);
		x2=x2r(ie)*(0.98+0.04*rand);
		x3=x3r(ie);
%    
    end    
    if(rand<0.1 && j>fix(nt/10))
%      
		x1=x1r(ie)*(0.98+0.04*rand);
		x2=x2r(ie)*(0.98+0.04*rand);
		x3=x3r(ie)+tp*(0.05-0.1*rand);
%    
    end    
%
	if(x2 > fup*tp )
       x2 = fup*tp;
    end    
%
	if(x2 < flow*tp)
       x2 = flow*tp;
    end    
%
	if(j==1)
        x1=0.;
    else
        if(x1<1.0e-12)
            x1=am*rand;
        end
    end
    x1=abs(x1);
%
	clear error;
    clear y;
%   
    ta=t-t(1);
	y=x1*sin(x2*ta-x3);
%   
    error=a-y;
%
	if(std(error)<errormax)
%
                syna=-(error-a);
                error_rth=error;
	            errormax=std(error);
				x1r(ie)=x1;
				x2r(ie)=x2;
				x3r(ie)=x3;
                yr=y;
%
                out4 = sprintf(' %6ld  %13.4e  %10.4g %9.4f %9.4f  ',j,errormax,x1,x2/tp,x3);
                disp(out4)    
%		
       figure(2)
       cnew=running_sum+yr;
       plot(t,original,t,cnew);
       legend ('Input Data','Synthesis',2);
       xlabel('Time(sec)');
%
    end
%				
end
%   
a=error_rth;
%
ave=mean(a);  
sd=std(a);
%	
out4=sprintf('\n  ave=%12.4g  sd=%12.4g \n',ave,sd);
disp(out4)    
%