vrs_engine(fi,ai,damp,df);

function[fn,a_vrs,rd_vrs] = vrs_engine(fi,ai,damp,df);
%
    disp(' Calculating VRS.... ')
    disp(' ')
%
	fn(1)=5.;
    oct=1./24.;
    tpi=2.*pi;
    last=length(fi);
%
	for i=1:1000   % natural frequency loop
%
%   absolute acceleration
%    fid = fopen('vrs.dat','w')
%
        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;
%
%       fprintf(fid,'%9.3g %9.3g %9.3g %9.3g %9.3g %9.3g t=%9.3g ai(j)=%9.3g df=%9.3g sum=%9.3g \n',fi(j),fn(i),rho,tdr,c1,c2,t,ai(j),df,sum);
%       
        end
    a_vrs(i)=sqrt(sum);
%    fclose(fid);           
%
%   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 );
%
			if(j==1)
               t=t/2.;
            end   
%
            sum = sum + t*ai(j)*df;
        end
         rd_vrs(i)=sqrt(sum)*386./(tpi^2.);
%
        if(fn(i) > 10000. )
            break;
        end
        if(fn(i) > 2.*fi(last) )
            break;
        end    
     	fn(i+1)=fn(i)*(2.^oct);   
     end