% A "m" file for Matlab. After making a "x_dat" file from "grfftm.f", you
% can plot waveforms using Matlab. A post script data is stored in the 
% "x.ps" file.

clear
close

load xyz1_dat;
size(xyz1_dat)
nt=xyz1_dat(1,1)
nob=xyz1_dat(1,2)

load x1_dat;
size(x1_dat)
load y1_dat;
size(y1_dat)
load z1_dat;
size(z1_dat)
load time_dat;
size(time_dat)
f=max(time_dat)
amaxx=max(x1_dat)*2.5
amaxy=max(y1_dat)*2.5
amaxz=max(z1_dat)*2.5

load wave11_dat; 
load wave12_dat; 
load wave13_dat; 
load wave14_dat; 
load wave15_dat; 

for j=1:nt,
  w_x(j)=x1_dat(j);
  w_y(j)=y1_dat(j);
  w_z(j)=z1_dat(j);
  wx(1,j)=wave11_dat(j,1);
  wx(2,j)=wave12_dat(j,1);
  wx(3,j)=wave13_dat(j,1);
  wx(4,j)=wave14_dat(j,1);
  wx(5,j)=wave15_dat(j,1);
  wy(1,j)=wave11_dat(j,2);
  wy(2,j)=wave12_dat(j,2);
  wy(3,j)=wave13_dat(j,2);
  wy(4,j)=wave14_dat(j,2);
  wy(5,j)=wave15_dat(j,2);
  wz(1,j)=wave11_dat(j,3);
  wz(2,j)=wave12_dat(j,3);
  wz(3,j)=wave13_dat(j,3);
  wz(4,j)=wave14_dat(j,3);
  wz(5,j)=wave15_dat(j,3);
end
ax=-1.2*min(w_x)
ay=-1.2*min(w_y)
az=-1.2*min(w_z)
for j=1:nt,
  w_x(j)=x1_dat(nt*(nob-1)+j);
  w_y(j)=y1_dat(nt*(nob-1)+j);
  w_z(j)=z1_dat(nt*(nob-1)+j);
end
bx=+1.2*max(w_x)
by=+1.2*max(w_y)
bz=+1.2*max(w_z)

subplot(1,3,1)
hold on
for i=1:nob,
  for j=1:nt,
    w_x(j)=x1_dat(nt*(i-1)+j)+ax+amaxx*(i-1)/2;
    w(j)=wx(i,j)+ax+amaxx*(i-1)/2;
  end
  plot(time_dat,w_x,'-',time_dat,w,'--')
end
title('X')
set(gca, 'box','on')
axis([0 f 0 ax+bx+amaxx*(nob-1)/2])
ylabel('velocity')
xlabel('time (sec)')
hold off

subplot(1,3,2)
hold on
for i=1:nob,
  for j=1:nt,
    w_y(j)=y1_dat(nt*(i-1)+j)+ay+amaxy*(i-1)/2;
    w(j)=wy(i,j)+ay+amaxy*(i-1)/2;
  end
  plot(time_dat,w_y,'-',time_dat,w,'--')
end
title('Y')
set(gca, 'box','on')
axis([0 f 0 ay+by+amaxy*(nob-1)/2])
ylabel('velocity')
xlabel('time (sec)')
hold off

subplot(1,3,3)
hold on
for i=1:nob,
  for j=1:nt,
    w_z(j)=z1_dat(nt*(i-1)+j)+az+amaxz*(i-1)/2;
    w(j)=wz(i,j)+az+amaxz*(i-1)/2;
  end
  plot(time_dat,w_z,'-',time_dat,w,'--')
end
title('Z')
set(gca, 'box','on')
axis([0 f 0 az+bz+amaxz*(nob-1)/2])
ylabel('velocity')
xlabel('time (sec)')
hold off

print wave3.ps


