function[a,x1r,x2r,x3r,error_rth,syna] = sfa_engine_60(a,t,drate,num2,flow,fup,ie,nt,x1r,x2r,x3r,ijk,nlast);
%
%    sfa_engine  ver 1.7   June 12, 2007
%
rand('state',0)
clear ta;
clear error_rth
clear syna;
%
tp=2*pi;
ave=mean(a);
sd=std(a);
%
am=2.*sd;
n=num2;
pu=1.;
pl=0.;
%
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=2.*am*x1;
	x2=((fup-flow)*x2+flow)*tp;
	x3=x3*tp;
%
    if(rand>0.4 & j>100)
        x2=x2r(ie);
    end
    if(rand>0.6 & j>fix(nt/10))
%      
		x1=x1r(ie)*(0.98+0.04*rand);
  		x2=x2r(ie)*(0.999+0.002*rand);
		x3=x3r(ie)+tp*(0.05-0.1*rand);
%    
    end    
%
    fd=fup-flow;
	if(x2 > fup*tp )
                x2 = (fd*rand+flow)*tp;
    end    
%
	if(x2 < flow*tp)
                x2 = (fd*rand+flow)*tp;
    end    
%
	if(j==1)
                x1=0.;
    else
        if(x1<1.0e-12)
            x1=am*rand;
        end
    end
    x1=abs(x1);
%
    if(ijk>1 && j==1)
      	x1=x1r(ie);
  		x2=x2r(ie);
		x3=x3r(ie);  
    end
    if(ijk==nlast)
        x1=x1r(ie);
    end
%
	clear error;
    clear y;
%        
    ta=t-t(1);
%    
%%    for(i=1:n)
%   
	y=x1*sin(x2*ta-x3);
%   
    ssaa=size(a);
    if(ssaa(1)>ssaa(2))
        a=a';
    end
    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;
%
                out4 = sprintf(' %6ld  %13.4e  %10.4g %9.4f %9.4f  ',j,errormax,x1,x2/tp,x3);
                disp(out4)    
%			
     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)