%
%   vrs_engine.m  ver 1.1  March 1, 2013
%
function[fn,a_vrs,rd_vrs] = vrs_engine(fi,ai,damp,df,iu)
%
    disp(' Calculating VRS.... ')
    disp(' ')
%
	fn(1)=5.;
    oct=1./24.;
    tpi=2.*pi;
    tpi_sq=tpi^2;
    last=length(fi);
%
	for i=1:1000   % natural frequency loop
%
%   absolute acceleration
%
        sum=0.; 
%
		for j=1:last 
%
		    rho = fi(j)/fn(i);
			tdr=2.*damp*rho;
%
            c1= tdr^ 2.;
			c2= (1.- (rho^2.))^ 2.;
%
			t= (1.+ c1 ) / ( c2 + c1 );
            sum = sum + t*ai(j)*df;
%       
        end
%        
        a_vrs(i)=sqrt(sum);
%    
%
%   relative displacement
%
		sum=0.; 
%
		for j=1:last 
%
            c1= ((fn(i)^2.)-(fi(j)^2.) )^2.;
			c2= ( 2.*damp*fn(i)*fi(j))^2.;
%
			t= 1. / ( c2 + c1 );   
%
            sum = sum + t*ai(j)*df;
        end
%
        if(iu==1)
            rd_vrs(i)=sqrt(sum)*386./tpi_sq;
        else
            rd_vrs(i)=sqrt(sum)*9.81*1000/tpi_sq;           
        end
%
        if(fn(i) > 10000. )
            break;
        end
        if(fn(i) > 2.*fi(last) )
            break;
        end    
     	fn(i+1)=fn(i)*(2.^oct);   
    end