function [accel, veloc, displ] = pulseVx ( Vp, Tp, Nc, t, t0, fig ) 
% [accel, veloc displ] = pulseVx ( Vp, Tp, Nc, t, t0, fig )  
% Computes an earthquake-like acceleration, velocity, and displacement pulse 
%     INPUT                                                  DEFAULT
%   =========                                                =======
%       Vp   - max velocity of pulse                           1.0
%       Tp   - time period of pulse                            1.0
%       Nc   - number of cycles in pulse ... approximate       1.0
%       t    - time record ... ( [1:points] * delta_t )        [1:1000]*0.01
%       t0   - time at which the pulse starts                  0.25
%       fig  - figure number for plotting, 0: no plots         1       
%
%    OUTPUT
%   =========
%      accel - earthquake ground acceleration
%      veloc - earthquake ground velocity
%      displ - earthquake ground displacement


 if nargin < 1,	Vp  =  1.0;	end
 if nargin < 2,	Tp  =  1.0;	end
 if nargin < 3,	Nc  =  1.0;	end
 if nargin < 4,	t    = (0:1000)*0.01; end
 if nargin < 5,	t0   = 0.25 ; end
 if nargin < 6, fig  = 1; end
 
 points  = length(t);			% number of points
 delta_t = t(2) - t(1);			% time step interval
 p0 = floor((t0-t(1))/delta_t);		% number of inital points w/o motion
 pF = points-p0;			% number of points w/ motion
 tau = Nc*Tp/4.0;			% decay time constant of pulse
 tT = tau/Tp;				% Nc / 4
 n  = 2.0;				% rise-time value for pulse

% set phase of pulse such that terminal velocity is zero, or close to it.
 phi = atan2( (3*(2*pi*tT)^2-1) , (3*2*pi*tT - (2*pi*tT)^3) );
 accel = ((t-t0)/tau).^n .* exp(-(t-t0)/tau) .* cos(2*pi*(t-t0)/Tp - phi);
 accel(1:p0) = 0;			% shift the offset

% baseline correction for zero terminal velocity and displacement 
  z = 1/4;				% correction over duration z*b
  ts =  (t-n*tau-t0)/(z*tau);		% scaled time
  exp_ts = exp(-ts);
  dc = 1./(1+exp_ts);			% sigmoidal displacement correction
% vc = exp_ts .* (1+exp_ts).^(-2.0)/(z*tau);                    % d(dc)/dt
  ac = exp_ts .* (1+exp_ts).^(-3.0).*(exp_ts-1)/(z*tau)^2;      % d(vc)/dt

  maxIter = 15;
  dspT = zeros(maxIter,1);
  for iter = 1:maxIter;
	dspT(iter)  = [points:-1:1]*accel'*delta_t^2;
    if ( abs(dspT(iter)) < 1e-5*Tp ); break; end
    accel = accel - dspT(iter)*ac;
  end

 veloc = cumtrapz(accel)*delta_t; % trapezoidal rule for the pulse veloc.
 displ = cumtrapz(veloc)*delta_t; % trapezoidal rule for the pulse displ.

%  if Nc > 5,
% 	maxV = 0.086 * Tp;
%  else
%  	maxV = max( [ 4.063*Nc^-2.165*exp(-4.403/Nc) ;
%                       2.329*Nc^-1.336*exp(-5.693/Nc) ] ) * Tp;
%  end

 maxV = max(abs(veloc));
 scale = Vp / maxV;
 accel = accel * scale;
 veloc = veloc * scale;
 displ = displ * scale;

if fig > 0 % plot the pulse data

   figure(fig) 
    subplot(3,1,1)
     plot (t,accel) %, t, scale*dspT*ac );
      axis tight
      ylabel('accel, m/s/s')
    subplot(3,1,2)
     plot (t,veloc)
      axis tight
      ylabel('veloc, m/s')
    subplot(3,1,3)
     plot (t,displ) %, t,scale*dspT*dc);
      axis tight
      ylabel('displ, m')
      xlabel('time, s')

   drawnow

end
%  pulseVx  2 Oct 2010; 18 Nov 2010; 24 May 2011; 19 Dec 2011
