home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 3.ddi / DEMOS.DI$ / QUAKE.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  6.1 KB  |  232 lines
  1. %QUAKE    Loma Prieta Earthquake.
  2. echo off
  3.  
  4. %    C. Denham, 1990, C. Moler, August, 1992.
  5. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  6.  
  7. clf
  8. clc
  9. echo on
  10.  
  11. % The file QUAKE.MAT contains 200Hz data from the October 17, 1989
  12. % Loma Prieta earthquake in the Santa Cruz Mountains.
  13. % The data are courtesy Joel Yellin at the Charles F. Richter
  14. % Seismological Laboratory, University of California, Santa Cruz.
  15. % The data can be loaded with the command:
  16.  
  17. load quake
  18.  
  19. % Press any key to continue after pauses.
  20. pause
  21.  
  22. clc
  23.  
  24. % In the workspace now are three variables containing time
  25. % traces from an accelerometer in the Natural Sciences'
  26. % building at UC Santa Cruz.  The accelerometer recorded the
  27. % main shock of the earthquake.  The variables n, e, v refer to
  28. % the three directional components measured by the instrument,
  29. % which was aligned parallel to the fault, with its N direction
  30. % pointing in the direction of Sacramento.
  31. % The data is uncorrected for the response of the instrument.
  32. % The data needs to be scaled by the gravitational acceleration.
  33. % A fourth variable, t, containing a time base, is also created.
  34.  
  35. g = 0.0980;
  36. delt = 1/200;
  37. e = g*e;
  38. n = g*n;
  39. v = g*v;
  40. t = delt*(1:length(e))';
  41. pause
  42.  
  43. clc
  44.  
  45. % Here are plots of the accelerations.
  46.  
  47. yrange = [-250 250];
  48. ax = [0 50 yrange];
  49. subplot(3,1,1), plot(t,e,'y-'), title('East-West acceleration'), axis(ax)
  50. subplot(3,1,2), plot(t,n,'m-'), title('North-South acceleration'), axis(ax)
  51. subplot(3,1,3), plot(t,v,'c-'), title('Vertical acceleration'), axis(ax)
  52. pause
  53.  
  54. clc
  55.  
  56. % Now, click twice with the mouse along one of the time axes
  57. % to select an interesting time interval.
  58. % Draw white lines at the selected spots.
  59. % All subsequent calculations will involve the selected interval.
  60.  
  61. t1 = ginput(1);
  62. t1(2) = t1(1);
  63. subplot(3,1,1), hold on, plot(t1,yrange,'w-','erase','none'); hold off
  64. subplot(3,1,2), hold on, plot(t1,yrange,'w-','erase','none'); hold off
  65. subplot(3,1,3), hold on, plot(t1,yrange,'w-','erase','none'); hold off
  66.  
  67. t2 = ginput(1);
  68. t2(2) = t2(1);
  69. subplot(3,1,1), hold on, plot(t2,yrange,'w-','erase','none'); hold off
  70. subplot(3,1,2), hold on, plot(t2,yrange,'w-','erase','none'); hold off
  71. subplot(3,1,3), hold on, plot(t2,yrange,'w-','erase','none'); hold off
  72. pause
  73.  
  74. clc
  75.  
  76. % Zoom in on the selected time interval.
  77.  
  78. trange = sort([t1(1) t2(1)]);
  79. k = find((trange(1)<=t) & (t<=trange(2)));
  80. e = e(k);
  81. n = n(k);
  82. v = v(k);
  83. t = t(k);
  84. ax = [trange yrange];
  85. subplot(3,1,1), plot(t,e,'y-'), title('East-West acceleration'), axis(ax)
  86. subplot(3,1,2), plot(t,n,'m-'), title('North-South acceleration'), axis(ax)
  87. subplot(3,1,3), plot(t,v,'c-'), title('Vertical acceleration'), axis(ax)
  88. pause
  89.  
  90. clf
  91. clc
  92.  
  93. % The COMET function produces a dynamic plot of one second's worth
  94. % of horizontal acceleration.
  95.  
  96. k = length(t);
  97. k = max(1,k/2-100):min(k,k/2+100);
  98. comet(e(k),n(k))
  99. xlabel('East'), ylabel('North');
  100. title('Acceleration During a One Second Period');
  101. pause
  102.  
  103. clc
  104.  
  105. % The velocity and position of a point in 3-D space can be obtained
  106. % by integrating the accelerations twice.
  107.  
  108. edot = cumsum(e)*delt;  edot = edot - mean(edot);
  109. ndot = cumsum(n)*delt;  ndot = ndot - mean(ndot);
  110. vdot = cumsum(v)*delt;  vdot = vdot - mean(vdot);
  111. subplot(2,1,1);
  112. plot(t,[edot+25 ndot vdot-25])
  113. axis([trange min(vdot-30) max(edot+30)])
  114. xlabel('Time'), ylabel('V - N - E'), title('Velocity')
  115. drawnow
  116.  
  117. epos = cumsum(edot)*delt;  epos = epos - mean(epos);
  118. npos = cumsum(ndot)*delt;  npos = npos - mean(npos);
  119. vpos = cumsum(vdot)*delt;  vpos = vpos - mean(vpos);
  120. subplot(2,1,2)
  121. plot(t,[epos+50 npos vpos-50])
  122. xlabel('Time'), ylabel('V - N - E'), title('Position')
  123. axis([trange min(vpos-55) max(epos+55)])
  124. pause
  125.  
  126. clf
  127. na = max(abs(npos)); na = 1.05*[-na na];
  128. ea = max(abs(epos)); ea = 1.05*[-ea ea];
  129. va = max(abs(vpos)); va = 1.05*[-va va];
  130. clc
  131.  
  132. % The trajectory defined by the position data can be displayed
  133. % with three different 2-dimensional projections.  The first is
  134.  
  135. subplot(2,2,1)
  136. plot(npos,vpos, 'y-');
  137. axis([na va])
  138. xlabel('North'); ylabel('Vertical');
  139. pause
  140.  
  141. % The trajectory can be annotated with a few values of t. 
  142.  
  143. nt = ceil((max(t)-min(t))/6);
  144. k = find(fix(t/nt)==(t/nt))';
  145. for j = k
  146.    text(npos(j),vpos(j),['o ' int2str(t(j))])
  147. end
  148. pause
  149.  
  150. clc
  151.  
  152. % Similar code produces two more 3-D views.
  153.  
  154. subplot(2,2,2)
  155. plot(epos,vpos, 'm-');
  156. for j = k;
  157.    text(epos(j),vpos(j),['o ' int2str(t(j))])
  158. end
  159. axis([ea va])
  160. xlabel('East'); ylabel('Vertical');
  161. drawnow
  162.  
  163. subplot(2,2,3)
  164. plot(npos,epos, 'c-');
  165. for j = k;
  166.    text(npos(j),epos(j),['o ' int2str(t(j))])
  167. end
  168. axis([na ea])
  169. xlabel('North'); ylabel('East');
  170. pause
  171.  
  172. clc
  173.  
  174. % The fourth subplot is a 3-D view of the trajectory.
  175.  
  176. subplot(2,2,4)
  177. plot3(npos,epos,vpos,'r-')
  178. for j = k
  179.    text(npos(j),epos(j),vpos(j),['o ' int2str(t(j))])
  180. end
  181. xlabel('North'); ylabel('East'), zlabel('Vertical');
  182. axis([na ea va])
  183. set(gca,'box','on')
  184. pause
  185.  
  186. clf
  187. clc
  188.  
  189. % Another 3-D view of position is obtained with the MESH function.
  190.  
  191. n = 50;
  192. p = fix((n-1)*(epos-min(epos))/(max(epos)-min(epos))+1);
  193. q = fix((n-1)*(npos-min(npos))/(max(npos)-min(npos))+1);
  194. r = vpos - min(vpos);
  195. z = zeros(n,n);
  196. for k = 1:length(r)
  197.    z(p(k),q(k)) = r(k);
  198. end
  199. mesh(z)
  200. axis([0 n 0 n 0 max(r)])
  201. title('The movement of a point in 3-dimensions')
  202. C = pink(64);
  203. colormap(C(33:64,:));
  204. pause
  205.  
  206. clf
  207. clc
  208.  
  209. % Finally, Handle Graphics allows us to produce a dancing dot which traces
  210. % the actual movement of a point on the earth during this portion of the quake.
  211.  
  212. plot3(npos,epos,vpos,'b')
  213. axis([min(npos) max(npos) min(epos) max(epos) min(vpos) max(vpos)])
  214. set(gca,'box','on')
  215. xlabel('North-South')
  216. ylabel('East-West')
  217. zlabel('Vertical')
  218. title('Position (cms)')
  219. hold on
  220. plt = plot3(0,0,0,'.','erasemode','xor','markersize',24);
  221. dk = ceil(length(epos)/1500);
  222. for k = 1:dk:length(epos)
  223.    set(plt,'xdata',npos(k),'ydata',epos(k),'zdata',vpos(k))
  224.    drawnow
  225. end
  226. hold off
  227. s = 'for k = 1:dk:length(epos)';
  228. s = [s 'set(plt,''xdata'',npos(k),''ydata'',epos(k),''zdata'',vpos(k)),'];
  229. s = [s 'drawnow,end'];
  230. echo off
  231. uicontrol('Units','normal','Position',[.4 .01 .1 .06],'String','Again','callback',s)
  232.