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

clear
%close

%fig=figure('Units','centimeters','Position',[0.0 0.0 21.0 29.5]...
%,'PaperUnits','centimeters','PaperPosition',[0.0 0.0 21.0 29.5]);

fig=figure('Units','centimeters','Position',[0.0 0.0 21.0 26.0]...
,'PaperUnits','centimeters','PaperPosition',[0.0 0.0 21.0 26.0]);

dur=10
fmax=5

load wave1_d;
size(wave1_d)
load amp1_dat;
size(amp1_dat)
nt=size(wave1_d,1);
nf=size(amp1_dat,1);

%x0=wave1_d(1,2);
%y0=wave1_d(1,3);
%z0=wave1_d(1,4);
%for i=2:nt,
%t(i)=wave1_d(i,1);
%x(i)=wave1_d(i,2)-x0;
%y(i)=wave1_d(i,3)-y0;
%z(i)=wave1_d(i,4)-z0;
%end

for i=1:nt,
t(i)=wave1_d(i,1);
x(i)=wave1_d(i,2);
y(i)=wave1_d(i,3);
z(i)=wave1_d(i,4);
end
xmax=max(x);ymax=max(y);zmax=max(z);
xmin=min(x);ymin=min(y);zmin=min(z);

for i=1:nf,
f(i)=amp1_dat(i,1)/3.141592654/2;
xr(i)=amp1_dat(i,2);
xi(i)=amp1_dat(i,3);
xa(i)=sqrt(xr(i)*xr(i)+xi(i)*xi(i))*100;
yr(i)=amp1_dat(i,4);
yi(i)=amp1_dat(i,5);
ya(i)=sqrt(yr(i)*yr(i)+yi(i)*yi(i))*100;
zr(i)=amp1_dat(i,6);
zi(i)=amp1_dat(i,7);
za(i)=sqrt(zr(i)*zr(i)+zi(i)*zi(i))*100;
end


dt=t(2)-t(1);
dx(1)=0.0;
dy(1)=0.0;
dz(1)=0.0;
for i=2:nt-1,
dx(i)=dx(i-1)+(x(i-1)+4*x(i)+x(i+1))*dt/6;
dy(i)=dy(i-1)+(y(i-1)+4*y(i)+y(i+1))*dt/6;
dz(i)=dz(i-1)+(z(i-1)+4*z(i)+z(i+1))*dt/6;
end
dx(nt)=dx(nt-1);
dy(nt)=dy(nt-1);
dz(nt)=dz(nt-1);
dx0=dx(nt);
dy0=dy(nt);
dz0=dz(nt);
dmax=max(dx); 
dmay=max(dy); 
dmaz=max(dz); 
dmix=min(dx); 
dmiy=min(dy); 
dmiz=min(dz); 
if -dmix > dmax
  dmax=-dmix
end
if -dmiy > dmay
  dmay=-dmiy
end
if -dmiz > dmaz
  dmaz=-dmiz
end

amax=dmax; 
if dmay > amax
  amax=dmay
end
if dmaz > amax
  amax=dmaz
end
amin=dmix; 
if dmiy < amin
  amin=dmiy
end
if dmiz < amin
  amin=dmiz
end
ad1=0.1*amax
ad2=-0.1*amin
if ad2 > ad1
  ad1=ad2
end
amax=amax+ad1
amin=amin-ad1

subplot(3,3,1)
plot(t,dx,'w-')
%title(['FN:Amax=',num2str(dmax),', d0=',num2str(dx0)])
%title(['NS:Amax=',num2str(dmax),', d0=',num2str(dx0)])
title(['FP:Amax=',num2str(dmax),', d0=',num2str(dx0)])
%title(['NS (max=',num2str(dmax),')'])
ylabel('Displacement (cm)')
xlabel('time (sec)')
axis([0 dur amin amax])

subplot(3,3,2)
plot(t,dy,'w-')
%title(['FP:Amax=',num2str(dmay),', d0=',num2str(dy0)])
%title(['EW:Amax=',num2str(dmay),', d0=',num2str(dy0)])
title(['FN:Amax=',num2str(dmay),', d0=',num2str(dy0)])
%title(['EW (max=',num2str(dmay),')'])
%ylabel('Displacement (cm)')
xlabel('time (sec)')
axis([0 dur amin amax])

subplot(3,3,3)
plot(t,dz,'w-')
title(['UD:Amax=',num2str(dmaz),', d0=',num2str(dz0)])
%title(['UD (max=',num2str(dmaz),')'])
%ylabel('Displacement (cm)')
xlabel('time (sec)')
axis([0 dur amin amax])

amax=xmax; 
if ymax > amax
  amax=ymax
end
if zmax > amax
  amax=zmax
end
amin=xmin; 
if ymin < amin
  amin=ymin
end
if zmin < amin
  amin=zmin
end
ad1=0.1*amax
ad2=-0.1*amin
if ad2 > ad1
  ad1=ad2
end
amax=amax+ad1
amin=amin-ad1

subplot(3,3,4)
plot(t,x,'w-')
%title(['FN:max=',num2str(xmax),', min=',num2str(xmin)])
%title(['NS:max=',num2str(xmax),', min=',num2str(xmin)])
title(['FP:max=',num2str(xmax),', min=',num2str(xmin)])
xlabel('time (sec)')
ylabel('Velocity (cm/s)')
axis([0 dur amin amax])

subplot(3,3,5)
plot(t,y,'w-')
%title(['FP:max=',num2str(ymax),', min=',num2str(ymin)])
%title(['EW:max=',num2str(ymax),', min=',num2str(ymin)])
title(['FN:max=',num2str(ymax),', min=',num2str(ymin)])
xlabel('time (sec)')
%ylabel('Velocity (cm/s)')
axis([0 dur amin amax])

subplot(3,3,6)
plot(t,z,'w-')
title(['UD:max=',num2str(zmax),', min=',num2str(zmin)])
xlabel('time (sec)')
%ylabel('Velocity (cm/s)')
axis([0 dur amin amax])

xmax=max(xa);
ymax=max(ya);
zmax=max(za);
amax=xmax; 
if ymax > amax
  amax=ymax
end
if zmax > amax
  amax=zmax
end
amax=amax*1.1

subplot(3,3,7)
plot(f,xa,'w-')
%title(['FN; Amp(0 Hz)=',num2str(xa(1))])
%title(['NS; Amp(0 Hz)=',num2str(xa(1))])
title(['FP; Amp(0 Hz)=',num2str(xa(1))])
xlabel('frequency (Hz)')
ylabel('Fourier Amplitude (cm)')
axis([0 fmax 0 amax])

subplot(3,3,8)
plot(f,ya,'w-')
%title(['FP; Amp(0 Hz)=',num2str(ya(1))])
%title(['EW; Amp(0 Hz)=',num2str(ya(1))])
title(['FN; Amp(0 Hz)=',num2str(ya(1))])
xlabel('frequency (Hz)')
ylabel('Fourier Amplitude (cm)')
axis([0 fmax 0 amax])

subplot(3,3,9)
plot(f,za,'w-')
title(['UD; Amp(0 Hz)=',num2str(za(1))])
xlabel('frequency (Hz)')
ylabel('Fourier Amplitude (cm)')
axis([0 fmax 0 amax])
