function[TT,a,np] = sweep(tmax);
%
%   sweep.m  version 1.0    September 13, 2005
%   By Tom Irvine   Email:  tomirvine@aol.com
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ');
num=input(' Enter the number of coordinates ');
%
for(i=1:num)
        out1=sprintf(' Enter frequency(Hz) %d:  ',i);
        fff(i)=input(out1);
        out2=sprintf(' Enter accel(G)      %d:  ',i);
        aaa(i)=input(out2);  
        disp(' ')
end
%
f1 = fff(1);
f2 = fff(num);
%
% sr = input('  Enter the sample rate (samples/seconds) ');
%
sr=12*f2;
%
disp(' ');
disp(' Enter sweep rate type ');
sweepr = input('  1=linear  2=log  ');
%
if( sr<12*f2)
    sr=12*f2;
end    
dt=1./sr;
ns=tmax*sr;
%  
oct=log(f2/f1)/log(2.);
%
out1 = sprintf('\n Number of octaves = %8.2f  ',oct);    
disp(out1);       
%
ntimes = ns;
cycles=0.;
t2=tmax;
dur=tmax;
%
tpi=2.*pi;
%
if(sweepr==2)
    rate=oct/dur;
%
    out1 = sprintf('\n Sweep rate = %8.2f oct/sec ',rate);    
    disp(out1);  
    out1 = sprintf('\n            = %8.2f  oct/min \n',60.*rate);    
    disp(out1);     
%
else
    rate=(f2-f1)/dur;
end
%
maxn=1000000;
%
if(ntimes>maxn);
    disp(' ')
    disp(' Time history truncated to 1000000 points, ')
    ntimes=maxn;  
    out5 = sprintf('or %g sec',maxn*dt);
    disp(out5);        
end
%
disp(' Calculating time vector... ')
TT=linspace(dt,(ntimes+1)*dt,ntimes); 
%
a = zeros(1,ntimes);
arg = zeros(1,ntimes);
freq = zeros(1,ntimes);
%
if(sweepr==2)  % log
    disp(' ')
    disp(' Calculating log sweep....')
    disp(' ')    
%  
%			fspectral = f1*pow(2.,rate*t);
%
    for(i=1:ntimes)
        arg(i)=-1.+2^(rate*TT(i));
    end
%    
else           % linear
    disp(' ')
    disp(' Calculating linear sweep....')
    disp(' ')     
%
% 0.5 factor is necessary to obtain correct number of cycles for linear case.
%
%			fspectral= (     rate*t ) + f1;
    fmax=0.5*(f2-f1)+f1;
    freq=linspace(f1,fmax,ntimes); 
end
%
if(sweepr==1) % linear
   arg=tpi*freq.*TT;
else  %log
   arg=tpi*f1*arg/(rate*log(2)); 
end
%
ntimes=length(a);
%
if(sweepr==1) % linear
    spectral=linspace(f1,f2,ntimes);   
else
    f1=log10(f1);
    f2=log10(f2);
    spectral=logspace(f1,f2,ntimes);   
end
%
limit=length(fff)-1;
amplitude=zeros(1,ntimes);
%
for(i=1:ntimes)
   for(j=1:limit)
       if(fff(j)<=spectral(i) & fff(j+1)>=spectral(i))
            x=spectral(i)-fff(j);
            L=fff(j+1)-fff(j);
            c2=x/L;
            c1=1-x/L;
            amplitude(i)=c1*aaa(j)+c2*aaa(j+1);
            break;
       end
   end
end
%
a=sin(arg);
%
a(1)=0.;
a=a.*amplitude;
TT=TT-dt;
%
% disp(' ')
% disp(' End of sweep ')
% disp(' ')
%
%% sine_sweep=[TT;a]'; 
%% out5 = sprintf('\n\n The time history is renamed as "sine_sweep" ');
%% disp(out5); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
np=ntimes;
%
disp(' ')
disp(' End of sweep ')
disp(' ')