洛马普里埃塔地震

1、命令行键入：load quake e n vwhos 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 onplot([t1 t2],yrange,'k','LineWidth',2)hold offsubplot(3,1,2)hold onplot([t1 t2],yrange,'k','LineWidth',2)hold offsubplot(3,1,3)hold onplot([t1 t2],yrange,'k','LineWidth',2)hold off

6、按“Enter”键。如图3所示。

7、命令行键入：trange = sort([t1(1) t2(1)]);k = find((trang髫潋啜缅e(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))]);endaxis([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))]);endaxis([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))]);endaxis([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 onstep = 10;plot3(npos(1:step:end),epos(1:step:end),vpos(1:step:end),'.')hold offbox onaxis tightxlabel('North-South')ylabel('East-West')zlabel('Vertical')title('Position (cms)')

20、按“Enter”键。如图10所示。