%QUAKE  Loma Prieta Earthquake.
echo off

%   C. Denham, 1990, C. Moler, August, 1992.
%   Copyright 1984-2000 The MathWorks, Inc. 
%   $Revision: 5.13 $  $Date: 2000/08/16 21:00:24 $

clf reset
clc
echo on

% The file QUAKE.MAT contains 200Hz data from the October 17, 1989
% Loma Prieta earthquake in the Santa Cruz Mountains.
% The data are courtesy of Joel Yellin at the Charles F. Richter
% Seismological Laboratory, University of California, Santa Cruz.
% The data can be loaded with the command:

load quake

% Press any key to continue after pauses.
pause

clc

% In the workspace now are three variables containing time
% traces from an accelerometer in the Natural Sciences'
% building at UC Santa Cruz.  The accelerometer recorded the
% main shock of the earthquake.  The variables n, e, v refer to
% the three directional components measured by the instrument,
% which was aligned parallel to the fault, with its N direction
% pointing in the direction of Sacramento.
% The data is uncorrected for the response of the instrument.
% The data needs to be scaled by the gravitational acceleration.
% A fourth variable, t, containing a time base, is also created.

g = 0.0980;
delt = 1/200;
e = g*e;
n = g*n;
v = g*v;
t = delt*(1:length(e))';
pause

clc

% Here are plots of the accelerations.

yrange = [-250 250];
ax = [0 50 yrange];
subplot(3,1,1), plot(t,e,'y-'), title('East-West acceleration'), axis(ax)
set(gca, 'Color', 'black');
subplot(3,1,2), plot(t,n,'m-'), title('North-South acceleration'), axis(ax)
set(gca, 'Color', 'black');
subplot(3,1,3), plot(t,v,'c-'), title('Vertical acceleration'), axis(ax)
set(gca, 'Color', 'black');
pause

clc

% Let's look at the interval from t=8 sec to t=15 sec
% Draw white lines at the selected spots.
% All subsequent calculations will involve the selected interval.

t1=8*[1 1];
subplot(3,1,1), hold on, plot(t1,yrange,'w-','erase','none', 'LineWidth', 2); hold off
subplot(3,1,2), hold on, plot(t1,yrange,'w-','erase','none', 'LineWidth', 2); hold off
subplot(3,1,3), hold on, plot(t1,yrange,'w-','erase','none', 'LineWidth', 2); hold off

%t2 = ginput(1)
%t2(2) = t2(1);
t2=15*[1 1];
subplot(3,1,1), hold on, plot(t2,yrange,'w-','erase','none', 'LineWidth', 2); hold off
subplot(3,1,2), hold on, plot(t2,yrange,'w-','erase','none', 'LineWidth', 2); hold off
subplot(3,1,3), hold on, plot(t2,yrange,'w-','erase','none', 'LineWidth', 2); hold off
pause

clc

% Zoom in on the selected time interval.

trange = sort([t1(1) t2(1)])
k = find((trange(1)<=t) & (t<=trange(2)));
e = e(k);
n = n(k);
v = v(k);
t = t(k);
ax = [trange yrange]

subplot(3,1,1), plot(t,e,'y-'), title('East-West acceleration'), axis(ax)
set(gca, 'Color', 'black');
subplot(3,1,2), plot(t,n,'m-'), title('North-South acceleration'), axis(ax)
set(gca, 'Color', 'black');
subplot(3,1,3), plot(t,v,'c-'), title('Vertical acceleration'), axis(ax)
set(gca, 'Color', 'black');
pause

clf
clc

% The COMET function produces a dynamic plot of one second's worth
% of horizontal acceleration.

k = length(t);
k = round(max(1,k/2-100):min(k,k/2+100));
xlabel('East'), ylabel('North');
title('Acceleration During a One Second Period');
comet(e(k),n(k))
pause

clc

% The velocity and position of a point in 3-D space can be obtained
% by integrating the accelerations twice.

edot = cumsum(e)*delt;  edot = edot - mean(edot);
ndot = cumsum(n)*delt;  ndot = ndot - mean(ndot);
vdot = cumsum(v)*delt;  vdot = vdot - mean(vdot);
subplot(2,1,1);
plot(t,[edot+25 ndot vdot-25])
axis([trange min(vdot-30) max(edot+30)])
xlabel('Time'), ylabel('V - N - E'), title('Velocity')
drawnow

epos = cumsum(edot)*delt;  epos = epos - mean(epos);
npos = cumsum(ndot)*delt;  npos = npos - mean(npos);
vpos = cumsum(vdot)*delt;  vpos = vpos - mean(vpos);
subplot(2,1,2)
plot(t,[epos+50 npos vpos-50])
xlabel('Time'), ylabel('V - N - E'), title('Position')
axis([trange min(vpos-55) max(epos+55)])
pause

clf
na = max(abs(npos)); na = 1.05*[-na na];
ea = max(abs(epos)); ea = 1.05*[-ea ea];
va = max(abs(vpos)); va = 1.05*[-va va];
clc

% The trajectory defined by the position data can be displayed
% with three different 2-dimensional projections.  The first is

subplot(2,2,1)
plot(npos,vpos, 'r-');
axis([na va])
xlabel('North'); ylabel('Vertical');
pause

% The trajectory can be annotated with a few values of t. 

nt = ceil((max(t)-min(t))/6);
k = find(fix(t/nt)==(t/nt))';
for j = k, text(npos(j),vpos(j),['o ' int2str(t(j))]); end
pause

clc

% Similar code produces two more 2-D views.

subplot(2,2,2)
plot(epos,vpos, 'color',[1 0 1],'linestyle','-');
for j = k; text(epos(j),vpos(j),['o ' int2str(t(j))]); end
axis([ea va])
xlabel('East'); ylabel('Vertical');
drawnow

subplot(2,2,3)
plot(npos,epos, 'color',[0 .5 0],'linestyle','-');
for j = k; text(npos(j),epos(j),['o ' int2str(t(j))]); end
axis([na ea])
xlabel('North'); ylabel('East');
pause

clc

% The fourth subplot is a 3-D view of the trajectory.
%quktest=-100
subplot(2,2,4)
plot3(npos,epos,vpos,'b-')
for j = k;text(npos(j),epos(j),vpos(j),['o ' int2str(t(j))]);end
xlabel('North'); ylabel('East'), zlabel('Vertical');
axis([na ea va])
set(gca,'box','on')
pause

clf
clc

% Another 3-D view of position is obtained with the MESH function.

n = 50;
p = fix((n-1)*(epos-min(epos))/(max(epos)-min(epos))+1);
q = fix((n-1)*(npos-min(npos))/(max(npos)-min(npos))+1);
r = vpos - min(vpos);
z = zeros(n,n);
echo off;
for k = 1:length(r);z(p(k),q(k)) = r(k);end
echo on;
mesh(z)
axis([0 n 0 n 0 max(r)])
set(gca, 'Color', 'black')
title('The movement of a point in 3-dimensions')
zlabel('Vertical');
xlabel('North-South')
ylabel('East-West')
colormap([0 0 1]);
pause

clf
clc

% Finally, Handle Graphics allows us to produce a dancing dot which traces
% the actual movement of a point on the earth during this portion of the quake.

echo off
plt = plot3(npos,epos,vpos,'b',0,0,0,'.');
plt = plt(2);
set(plt,'erasemode','xor','markersize',24);
axis([min(npos) max(npos) min(epos) max(epos) min(vpos) max(vpos)])
xlabel('North-South')
ylabel('East-West')
zlabel('Vertical')
title('Position (cms)')
dk = ceil(length(epos)/1500);
s = [];
s = [s 'for k = 1:dk:length(epos), '];
s = [s 'set(plt,''xdata'',npos(k),''ydata'',epos(k),''zdata'',vpos(k)), '];
s = [s 'drawnow, end'];
eval(s)
uicontrol('Units','normal','Position',[.4 .01 .1 .06], ...
          'String','Again','callback',s)
