function [accel,veloc,displ,scale,Ag,Bg,Cg] = eqgm_1D(PGA,fg,zg,aa,tt,t,figNo,seed)
% [accel, veloc, displ] = eqgm_1D ( PGA, fg, zg, aa, tt, t, figNo, seed )
%
% Generate an artificial (synthetic) earthquake ground motion record
%
%  Input Variable        Description                            Default
%  --------------        -----------                            -------
%      PGA       peak ground acceleration                       3.5
%      fg        ground frequency, Hz				1.5
%      zg        ground damping ratio				0.9
%      aa        envelope rise-time parameter			4.0
%      tt        envelope decay parameter, s			2.0
%      t         the time axis, sec    			0.01 * [ 1 : 3000] 	
%      figNo     figure number; 0: don't draw plots		9
%      seed      a seed for Gaussean noise generation
%
% suggested ground motion parameters
%
%median peak  spectrum ground motion parameters
%
%median peak  spectrum ground motion parameters
%
% peak     peak                                      peak 
% ground   ground    ground  ground  rise    decay   ----     ATC-63
% accel    veloc     freq.   damping expont  time     RMS     data 
% m/s/s     m/s      Hz      -       s       s       dt=0.01  set   
% -------  --------  ------  ------  ------  ------  -------  -------------
% PGA=3.5; PGV=0.33; fg=1.5; zg=0.9; aa=4.0; tt=2.0; P2R=2.7; % ATC-63-FF 
% PGA=5.0; PGV=0.52; fg=1.3; zg=1.1; aa=3.0; tt=2.0; P2R=2.7; % ATC-63-NFNP
% PGA=4.7; PGV=0.80; fg=0.5; zg=1.8; aa=1.0; tt=2.0; P2R=2.5; % ATC-63-NFP
%
%  Output Variable        Description
%  --------------         -----------
%     accel               simulated earthquake ground acceleration  (1xn)
%     veloc               simulated earthquake ground velocity      (1xn)
%     displ               simulated earthquake ground displacement  (1xn)

% Henri Gavin, Civil and Environmental Engineering, Duke Univsersity, 
% Jul 2007, Nov 2007, Mar 2008, Apr 2008, Nov 2009, Jan 2010, Mar 2010, Apr 2010
% Aug 2010, Dec 2010, Aug 2011, Aug 2019
%G = 9.81633;		% gravitational acceleration, g, m/s/s 
% --- default values ---

 if nargin < 1, PGA = 3.5; end  % peak ground acceleration, value
 if nargin < 2, fg  = 1.5; end  % ground frequency value, Hz
 if nargin < 3, zg  = 0.9; end  % ground damping ratio
 if nargin < 4, aa  = 4.0; end  % envelope rise time constant
 if nargin < 5, tt  = 2.0; end  % envelope decay time constant, sec
 if nargin < 6, t   = 0.01 * [ 1:3000 ]; end  % time step constant
 if nargin < 7, figNo = 0; end  % do not draw plots
 if nargin > 7,	% random number generator seed
	rand ('state', seed);
	randn('state', seed);
 end

 dt = t(2)-t(1);		% time step, s

 t0 = 1.0;			% initial time of no motion, s
 T  = t0+5.74*aa^0.42*tt;	% total duration of earthquake record
 if ( max(t) < T && figNo > 0 )
 	P  = floor(T/dt);
	printf('eqgm_1D: time series may not be long enough, T=%3.0f, P=%d\n', T,P);
 end
 P = length(t);

 Pz = floor(t0/dt);
 Pr = round(Pz / 2);

% envelope of time series
  envlp = [zeros(1,Pz-1), ((t(Pz:P)-t0)/(aa*tt)).^aa .* exp(aa-(t(Pz:P)-t0)/tt)];
  envlp = envlp + 0.005*PGA;
  envlp = envlp .* [ 1-cos([1:Pr]*pi/Pr), 2*ones(1,P-2*Pr),  1-cos([Pr:-1:1]*pi/Pr)] * 0.5;

  % model P2R based on fg*tt, z, and aa
  Tpk   = aa*tt;
  ttfg  = tt*fg;
  ttfg  = Tpk*fg;
  lzg   = log10(zg);
  c = [   2.0327e+00
          2.2535e-01
          1.9215e-02
          2.7183e-01
          2.3593e-02
         -3.1858e-02
         -4.4708e-03
         -2.4680e-02
         -6.3575e-04
         -2.8971e-01 ];
  
  B = [ 1, ttfg, aa, lzg, ttfg*aa, ttfg*lzg, aa*lzg, ttfg^2, aa^2, lzg^2 ];
  
  P2R = B*c; % peak to RMS ratio    
  
RMSa  = PGA / P2R; 
a_bar = PGA / P2R  / sqrt(2*pi*zg*fg);

% linear state variable model
  Ag = [ 0.0 1.0 ; -4*pi^2*fg^2  -4*pi*zg*fg ]; 
  Bg = [ 0.0 ; a_bar ];
  Cg = [ 0.0 , 4*pi*zg*fg ]; 

RMSa = sqrt(diag(real(Cg*liap(Ag,Bg*Bg')*Cg')));    % RMS input acceleration

% damp(Ag)					% check dynamics matrix
 
% Unit Gaussian white noise has a power spectral density value of 1.0
% and a mean square value of 1/dt.
 u = envlp .* randn(1,P) / sqrt(dt);

 accel = lsym( Ag, Bg, Cg, 0, u, t);

 figure(1);
 clf;
 hold on;
 [maxA,idx] = max(abs(accel));
 plot(t,accel);
 plot(t, envlp*RMSa);
 plot([0 P]*dt, [1 1]*RMSa, '--k');
 plot([0 P]*dt,-[1 1]*RMSa, '--k');
 plot([0 P]*dt, [1 1]*PGA, '--', 'color', [0,0.8,0] );
 plot([0 P]*dt,-[1 1]*PGA, '--', 'color', [0,0.8,0] );
 plot(t(idx), accel(idx), '*', 'MarkerSize', 10, 'color', [0,0.8,0] )
 hold off

%rmsA = std(accel);

% remove some bias and drift from the acceleration
% with aa = 1.0*dt , the response time of the filter is  5 s
% with aa = 0.5*dt , the response time of the filter is 10 s
% with aa = 0.2*dt , the response time of the filter is 25 s
% with aa = 0.1*dt , the response time of the filter is 50 s
 t0 = 1.2*max(find(abs(accel)/max(abs(accel)) > 0.98 ))*dt;
 tt = 0.3;
 accel = accel2displ(accel, t, 'SRA', 0.5*dt, t0, tt);

 veloc = cumtrapz(accel)*dt;    % veloc from detrended accel 
 displ = cumtrapz(veloc)*dt;    % displ from detrended accel 

% normalize to specified median PGA
%{
  Amax = max(abs(accel));
  scale = PGA / Amax;
  accel = accel * scale;
  veloc = veloc * scale;
  displ = displ * scale;
  Bg    = Bg    * scale;
%}

% normalize to specified median PGV
%{
  Vmax = max(abs(veloc));
  scale = PGA / Vmax;
  accel = accel * scale;
  veloc = veloc * scale;
  displ = displ * scale;
  Bg    = Bg    * scale;
%}

% if ( max(abs(veloc)) < 0.99*PGV || max(abs(veloc)) > 1.01*PGV ) 
%fprintf(' eqgm_1D: warning - velocity not scaled to PGV\n');
% end

  Amax = max(abs(accel));
  Vmax = max(abs(veloc));
  Dmax = max(abs(displ));

% eqgm_data = [ accel ; veloc ; displ ];

 if figNo > 0	% --- plot the earthquake ground motion
   
   T = 30;      % maximum time displayed in the plots

   figure(figNo) 
    clf
    subplot(311)
     hold on
     plot( t,Amax*envlp,'-r') 
     plot( t,-Amax*envlp,'-r')
     plot( t,accel,'-b' )
     hold off
     ylabel('accel')
     axis('tight')
    subplot(312)
     plot( t, veloc,'-b' )
     ylabel('veloc')
     axis('tight')
    subplot(313)
     plot(  t, displ,'-b' )
     ylabel('displ')
     axis('tight')


%{
  figure(figNo+2)	 	% double-check ground motion model dynamics
   f = [0.05:0.15:14];
   [mag,pha] = mimoBode(Ag,Bg,Cg,0,2*pi*f);	% check transfer function
   clf
   f1 = fg * sqrt(1+2*zg^2 - 2*zg*sqrt(1+zg^2));
   f2 = fg * sqrt(1+2*zg^2 + 2*zg*sqrt(1+zg^2));
   plot( [fg fg],[0,0.95], '-r', [f1,f2],[.707,.707],'-r', f,mag/max(mag), '-b')
   text(1.05*fg,0.05,'f_g', 'color', 'red', 'FontSize',14)
   text((f1+f2)/2,0.76,'2 f_g \zeta_g', 'color', 'red', 'FontSize',14)
   ylabel('acceleration amplitude spectrum')
   xlabel('frequency, Hz')
   axis( [ 0.00 10 0 1.05] )
% zero_2zgfg = f2 - f1 - 2*fg*zg
%}
 drawnow

 end 	        % --- plot the earthquake ground motion

% ---------------------------------------------------------------- EQGM_1D
% Nov  7 2009, 
% Jan 24 2010, Jan 27 2010, Jan 28 2010, Mar 15 2010, Apr 6 2010, Aug 13 2010,
% Jun 16 2011, Aug 15 2011, Dec 2 2014, Aug 2019, Aug 2020
