function [accel, veloc, displ,N,t] = chirp( ao,af, fo,ff, cycles, t, p,n, figNo,units)
% [accel,veloc,displ] = chirp ( Ao, Af, Fo, Ff, CYCLES, t, P, n, figNo, units )  
% Compute a sine-sweep time record of linearly changing frequency and amplitude
%	Ao - starting velocity amplitude
%	Af - stopping velocity amplitude 
%	Fo - starting frequency , Hz 
%	Ff - stopping frequency , Hz 
%   CYCLES - number of cycles of data in the sweep
%       t  - time values
%	P  - power of geometric increase in frequency, default = 1
%       n  - characteristic time for amplitude variation, default = 1
%    figNo - figure number of plotted responses
%    units - units of Ao and Af
%
% The  first output of SSWEEP contains an 'acceleration' record.
% The second output of SSWEEP contains a    'velocity'   record.
% The third  output of SSWEEP contains a  'displacement' record.

% Henri Gavin, Dept. Civil Engineering, Duke Univ., 6-99 henri.gavin@duke.edu

 if ( nargin < 8 ) p = 2; n = 1; end
 if ( nargin < 9 ) figNo = 0; end
 if ( nargin < 10 ) units = 'm'; end

dt = t(2)-t(1);

% ---  more cycles at lower frequency ...
 p = abs(p);
 if ( fo < ff && p < 1 )  
	p = 1/p;	
 end
 if ( fo > ff && p > 1 )  
	p = 1/p;	
 end

 T = cycles / ( fo + (ff-fo)/(p+1));	% length of time in the record

 N = round(T/dt);
 t = [ 0:N-1]*dt;

% ------- amplitude linear in time
% amp = ao +  [1:N] * (af-ao) * dt/T; 	% linear amplitude

% ------- amplitude a power law in time
% q = log(ao/af) / log(1+T);
% amp = ao ./ ( 1 + [1:N]*dt ).^q;       % power-law amplitude

% ------- amplitude exponential in time
% n = 1;
 r = (1/T^n) * log(ao/af);
 amp = ao *  exp(-r*([1:N]*dt).^n);	% exponential amp.

% ------- phase envelope

 phase = 2*pi * [ t*fo + t.^(p+1) * (ff-fo) / ((p+1)*T^p) ];

 veloc = taper ( amp .* sin(phase), floor(N/10), floor(N/10) );

 accel = cdiff(veloc) / dt;
 displ = cumsum( veloc ) * dt;

 fprintf('     %d cycles in %f seconds with %d points\n',  ...
			cycles, length(veloc)*dt, length(veloc) );

if  (figNo > 0) % plot chirp record

  figure(figNo); 
   subplot(311)
    plot(t,accel )
     ylabel(sprintf('acceleration, %s/s^2',units))
     axis('tight')
   subplot(312)
    plot(t,veloc)
     ylabel(sprintf('velocity , %s/s',units))
     axis('tight')
   subplot(313)
    plot(t,displ )
     xlabel('time, s')
     ylabel(sprintf('displacement, %s',units))
     axis('tight')

end

% endfunction # -------------------------------------------------------- CHIRP
