function [dep,vit,acc,qs] = NewmarkMDDL(m,k,c,p,DT,Nstep,a,d,depin,vitin)
% [dep,vit,acc,qs] = NewmarkMDDL(m,k,c,p,DT,Nstep,a,d)

nddl = length(k);

dep = zeros(nddl,Nstep);
vit = zeros(nddl,Nstep);
acc = zeros(nddl,Nstep);
qs  = zeros(nddl,Nstep);

if nargin>8
    dep(:,1)=depin;
    vit(:,1)=vitin;
end

acc(:,1) = m\(p(:,1)-c*vit(:,1)-k*dep(:,1));

akf = m/a/DT^2 + d/a/DT*c + k;
invakf = inv(akf);
invk = inv(k);

for i=2:Nstep
    dep(:,i) =invakf*(p(:,i)+m*(1/a/DT^2*dep(:,i-1)+1/a/DT*vit(:,i-1)+(1/2/a-1)*acc(:,i-1)) ...
        + c*(d/a/DT*dep(:,i-1) + (d/a-1)*vit(:,i-1) + DT/2*(d/a-2)*acc(:,i-1)));
    vit(:,i) = d/a/DT*(dep(:,i)-dep(:,i-1))+(1-d/a)*vit(:,i-1)+DT*(1-d/2/a)*acc(:,i-1);
    acc(:,i) = 1/a/DT^2*(dep(:,i)-dep(:,i-1))-1/a/DT*vit(:,i-1)-(1/2/a-1)*acc(:,i-1);
    qs(:,i) = invk*p(:,i);
end
