
clear
clf
pause off
hold off

% the seismometer equation


% some initialisations
f0=1;  % eigenfrequency of seismometer
w=2*pi*f0; % ... in radians
dt=.01; % time increment for numerical scheme
fu0=input('Give frequency of forcing (e.g. f=1 Hz) : '); % central frequency of forcing
h=input('Give damping parameter (e.g. h=0.5) : '); % damping 
% initialize ground motion
p=1./fu0;
nts=2*p/dt;
uii=1:nts;
uii=uii*0;
t0=p/dt;
a_ricker=4/p;
for it=1:nts
            t=(it-t0)*dt;uii(it) = -2*a_ricker*t*exp(-(a_ricker*t)^2);
end
% end initialize ground motion

nt=round(5*1/fu0/dt); % total number of time steps

src=(1:nt)*0;
src(1:length(uii))=uii;

x=(1:nt)*0;

% initial condition

eps=h*w;
x=0;
xnow=0;
xold=0;

for i=1:nt,
if i==1,xold=xnow;end

xnew=1./(1./dt^2+2*eps/(2*dt)) * (-src(i) -w^2*xnow + 2*eps/(2*dt) * xold - (xold-2*xnow)/dt^2 );
xold=xnow;
xnow=xnew;
x(i)=xnow;  

plot((1:i)*dt,x(1:i)/max(abs(x(1:i))),'r-','LineWidth',1.5),hold on,
plot((1:i)*dt,src(1:i)/max(src(1:i)),'b-','LineWidth',1.5),hold off
title(sprintf('F0=1Hz, SRC=%g Hz, h=%g ',fu0,h));
axis([0 nt*dt -1 1]), pause(0.5),drawnow
xlabel('Time (s)')
ylabel('Displacement')

end

