function [dep,vit,acc,qs] = Newmark(m,ksi,f,p,DT,Nstep,a,d,dep0,vit0)
%[dep,vit,acc,qs] = Newmark(m,ksi,f,p,DT,Nstep,a,d,dep0,vit0)

if nargin<9; dep0=0; vit0=0; end
if nargin<7; a=0.25; d=0.50; end

omega = 2*pi*f;
k = m * omega^2;
c = 2*m*omega*ksi;

dep=zeros(Nstep+1,1); dep(1)=dep0;
vit=zeros(Nstep+1,1); vit(1)=vit0;
acc=zeros(Nstep+1,1);
akf = m/a/DT^2 + d/a/DT*c + k;

for i=2:Nstep+1
    dep(i) = (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))) / akf;
    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) = p(i)/k;
end
