% **************************************************************************
%                              
%                            ***** TMA1 *****
%
%                   Prepared by D.Bernal - February, 1999.
%
% **************************************************************************

% Program computes the response of linear MDOF systems using the transition
% matrix approach. The solution is exact for the assumed variation of the
% load within the step. Two load variations are possible
% 1) Linear (default)
% 2) Constant (Requires setting flag='cst')


% INPUT
% m,c,k = mass dampig and stiffness matrices for the system.
% dt = time step
% b2 = load influence matrix. Column j contains the spatial distribution of input j.
% p = inputs. Row j contains the time history of input j.
% d0,v0 = initial conditions
% nsteps = number of time steps in the solution.  
% lvs -- if lvs ='' the solution is carried out assuming a linear variation of
% the load within the step. Set lv='cst' for constant load within the step.

% OUTPUT
% u = matrix of displacements
% ud = matrix of velocities
% udd = matrix of accelerations (if flag=7);
% (each column contains one DOF);

% LIMITATIONS
% Mass matrix must be positive definite
% ***************************************************************************

function [u,ud,udd]= tma1(m,c,k,dt,b2,p,d0,v0,nsteps,lvs,dtf);

% Make sure that the length of specified loading is adequate
[dum,ns]=size(p);
if nsteps*dt>ns*dtf;
   disp('nsteps*dt is larger than the time over which the loading is defined');
   break;
else
end;

% Interpolate or extrapolate the applied load to match the integration time step.
if dtf/dt~=1;
x=[0:dtf:(ns-1)*dtf]';
xi=[0:dt:nsteps*dt]';
F=interp1(x,p',xi,'spline');
F=F';
else
   F=p;
   end;
p=F;
% To eliminate the computation of accelerations flag~=7;
flag=7;

% Basic matrices.
[dof,dum]=size(m);
A=[eye(dof)*0 eye(dof);-inv(m)*k -inv(m)*c];
B=[eye(dof)*0 eye(dof)*0;eye(dof)*0 inv(m)];

% Place d0,v0 in column form in case they where input as a row vectors
[r,cx]=size(d0);
if r==1;
d0=d0';
end;
[r,cx]=size(v0);
if r==1;
v0=v0';
end;

% Compute the matrices used in the marching algorithm.
AA=A*dt;
F=expm(AA);
X=inv(A);
P1=X*(F-eye(2*dof))*B;
P2=X*(B-P1/dt);

% Set P2 to zero if the load is constant within the step
if lvs=='cst'
P2=P2*0;
else
   end;

% Initialize.
yo=[d0;v0];

% Perform the integration.
for i=1:nsteps;
y1(:,i)=(P1+P2)*[d0*0;b2*p(:,i)]-P2*[d0*0;b2*p(:,(i+1))]+F*yo;
yo=y1(:,i);
end;

u=y1(1:dof,:)';
u=[d0';u];
ud=y1(dof+1:2*dof,:)';
ud=[v0';ud];

% Calculate the acceleration vector if requested.
if flag==7;
X=inv(m);
for i=1:nsteps+1;
uddi=X*(b2*p(:,i)-c*ud(i,:)'-k*u(i,:)');
udd=[udd uddi];
end;
udd=udd';
end;