洛馬普里埃塔地震

這個例子展示了如何分析和可視化真實世界的地震數據。Authors: C. Denham, 1990, C. Moler, August, 1992.

東西/原料

matlab軟件
電腦

方式/步調

1
號令行鍵入：
load quake e n v
whos e n v
2
按“Enter”鍵。
如圖1所示。
3
號令窗口鍵入：
g = 0.0980;
e = g*e;
n = g*n;
v = g*v;
delt = 1/200;
t = delt*(1:length(e))';
這是加快度的曲線圖。
號令行鍵入：
yrange = [-250 250];
limits = [0 50 yrange];
subplot(3,1,1), plot(t,e,'b'), axis(limits), title('East-West acceleration')
subplot(3,1,2), plot(t,n,'g'), axis(limits), title('North-South acceleration')
subplot(3,1,3), plot(t,v,'r'), axis(limits), title('Vertical acceleration')
4
按“Enter”鍵。
如圖2所示。
5
號令行鍵入：
t1 = 8*[1;1];
t2 = 15*[1;1];
subplot(3,1,1)
hold on
plot([t1 t2],yrange,'k','LineWidth',2)
hold off
subplot(3,1,2)
hold on
plot([t1 t2],yrange,'k','LineWidth',2)
hold off
subplot(3,1,3)
hold on
plot([t1 t2],yrange,'k','LineWidth',2)
hold off
6
按“Enter”鍵。
如圖3所示。
7
號令行鍵入：
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,'b')
axis(ax)
title('East-West acceleration')
subplot(3,1,2)
plot(t,n,'g')
axis(ax)
title('North-South acceleration')
subplot(3,1,3)
plot(t,v,'r')
axis(ax)
title('Vertical acceleration')
8
按“Enter”鍵。
如圖4所示。
9
號令行鍵入：
subplot(1,1,1)
k = length(t);
k = round(max(1,k/2-100):min(k,k/2+100));
plot(e(k),n(k),'.-')
xlabel('East'), ylabel('North');
title('Acceleration During a One Second Period');
10
按“Enter”鍵。
如圖5所示。
11
將加快度積分兩次，計較三維空間中點的速度和位置。
號令行鍵入：
edot = cumsum(e)*delt;
edot = edot - mean(edot);
ndot = cumsum(n)*delt;
ndot = ndot - mean(ndot);
vdot = cumsum(v)*delt;
vdot = vdot - mean(vdot);

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,1);
plot(t,[edot+25 ndot vdot-25]);
axis([trange min(vdot-30) max(edot+30)])
xlabel('Time'), ylabel('V - N - E'), title('Velocity')

subplot(2,1,2);
plot(t,[epos+50 npos vpos-50]);
axis([trange min(vpos-55) max(epos+55)])
xlabel('Time'), ylabel('V - N - E'), title('Position')
12
按“Enter”鍵。
如圖6所示。
13
由位置數據界說的軌跡可以用三個分歧的二維投影顯示出來。下面是第一個帶有t注解的值的例子。
號令行鍵入：
subplot(1,1,1);
cla;
subplot(2,2,1)
plot(npos,vpos,'b');
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];
axis([na va]);
xlabel('North');
ylabel('Vertical');

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
14
按“Enter”鍵。
如圖7所示。
15
近似的代碼發生了別的兩個二維視圖。
號令行鍵入：
subplot(2,2,2)
plot(epos,vpos,'g');
for j = k
text(epos(j),vpos(j),['o ' int2str(t(j))]);
end
axis([ea va]);
xlabel('East');
ylabel('Vertical');

subplot(2,2,3)
plot(npos,epos,'r');
for j = k
text(npos(j),epos(j),['o ' int2str(t(j))]);
end
axis([na ea]);
xlabel('North');
ylabel('East');
16
按“Enter”鍵。
如圖8所示。
17
第四個子情節是軌跡的三維視圖。
號令行鍵入：
subplot(2,2,4)
plot3(npos,epos,vpos,'k')
for j = k
text(npos(j),epos(j),vpos(j),['o ' int2str(t(j))]);
end
axis([na ea va]);
xlabel('North');
ylabel('East')
zlabel('Vertical');
box on
18
按“Enter”鍵。
如圖9所示。
19
最后，在每10個點上畫一個點，點之間的間距暗示速度。
號令行鍵入：
subplot(1,1,1)
plot3(npos,epos,vpos,'r')
hold on
step = 10;
plot3(npos(1:step:end),epos(1:step:end),vpos(1:step:end),'.')
hold off
box on
axis tight
xlabel('North-South')
ylabel('East-West')
zlabel('Vertical')
title('Position (cms)')
20
按“Enter”鍵。
如圖10所示。
END