function [X,Xd,Xdd,t] = newmark(M,C,K,f,dt,gamma,beta,Xi,Xdi) 
%================================================================ 
%function [X,t] = newmark(M,C,K,f,dt,gamma,beta,Xi,Xdi) 
% 
%  (c) Ajax, (all rights reserved) 
%  Central South University 
%  October 13, 2006 
% 
%  This function is intended to perform the numerical 
%  integration of a structural system subjected to an 
%  external dynamic excitation such as a wind or earthquake. 
%  The structural model is assumed to be a lumped mass shear 
%  model.The integration scheme utilized for this analysis 
%  is the newmark alpha-beta method. The newmark alpha-beta 
%  method is an implicit time steping scheme so stability of 
%  the system need not be considered.
% 
%  Input Variables: 
%  [M] = Mass Matrix (nxn) 
%  [C] = Damping Matrix (nxn) 
%  [K] = Stiffness Matrix (nxn) 
%  {f} = Excitation Vector (mx1) 
%  dt = Time Stepping Increment 
%  beta= Newmark Const (1/6 or 1/4 usually) 
%  gamma = Newmark Const (1/2) 
%  Xi = Initial Displacement Vector (nx1) 
%  Xdi = Initial Velocity Vector (nx1) 
% 
%  Output Variables: 
%  {t} = Time Vector (mx1) 
%  [X] = Response Matrix (mxn) 
%================================================================ 

  n = size(M,1); 
  fdimc = size(f,2); 
  fdimr = size(f,1); 

%  Check Input Excitation 
%  ====================== 
   if(fdimc==n)
       f=f';
   end
   m=size(f,2);
   
%  Coefficients 
%  ============ 
        c0 = 1/(beta*dt*dt) ; 
        c1 = gamma/(beta*dt) ; 
        c2 = 1/(beta*dt) ; 
        c3 = 1/(beta*2) - 1 ; 
        c4 = gamma/beta - 1 ; 
        c5 = 0.5*dt*(gamma/beta - 2 ) ; 
        c6 = dt*(1 - gamma ) ; 
        c7 = dt* gamma ; 

%Initialize Stiffness Matricies 
% ============================== 
        Keff = c0*M + c1*C + K ; 
        Kinv = inv(Keff) ;

         
% Initial Acceleration
% ==================== 
        R0=f(:,1);
        Xddi= inv(M)*(f - K*Xi - C*Xdi)

% Perform First Step 
% ==================
        f(:,1) = f(:,1) + M*(c0*Xi+c2*Xdi+c3*Xddi) ... 
                 +C*(c1*Xi+c4*Xdi+c5*Xddi) ; 
        X(:,1) = Kinv*f(:,1) ; 
        Xdd(:,1)= c0*(X(:,1)-Xi) - c2*Xdi - c3*Xddi ; 
        Xd(:,1) = Xdi + c6*Xddi + c7*Xdd(:,1) ; 

%Perform Subsequent Steps 
% ======================== 
        for i=1:size(f,2)-1; 
          f(:,i+1) = f(:,i+1) + M * ...
        (c0*X(:,i)+c2*Xd(:,i)+c3*Xdd(:,i)) ... 
                              + C*(c1*X(:,i)+c4*Xd(:,i)+c5*Xdd(:,i)) ; 
          X(:,i+1) = Kinv*f(:,i+1) ; 
          Xdd(:,i+1)= c0*(X(:,i+1)-X(:,i))-c2*Xd(:,i)-c3*Xdd(:,i) ; 
          Xd(:,i+1) = Xd(:,i)+c6*Xdd(:,i)+c7*Xdd(:,i+1) ; 
        end; 

         
%  Strip Off Padded Response Zeros 
%  =============================== 
  X = X'
  Xd=Xd'
  Xdd=Xdd'
  
%  Generate the Time Vector 
%  ======================== 
  for i=0:1:size(X,1)-1 
    t(i+1,1) = i*dt; 
  end; 
  t