% 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]);

load wave1_dat;
load amp1_dat;

load wave_dat;
load amp_dat;

nt1=size(wave1_dat,1);
nf1=size(amp1_dat,1);
nt2=size(wave_dat,1);
nf2=size(amp_dat,1);

for i=1:nt1,
t1(i)=wave1_dat(i,1);
x1(i)=wave1_dat(i,2);
y1(i)=wave1_dat(i,3);
z1(i)=wave1_dat(i,4);
end
xmax1=max(x1);ymax1=max(y1);zmax1=max(z1);
xmin1=min(x1);ymin1=min(y1);zmin1=min(z1);

for i=1:nt2,
t2(i)=wave_dat(i,1);
x2(i)=wave_dat(i,2);
y2(i)=wave_dat(i,3);
z2(i)=wave_dat(i,4);
end
xmax2=max(x2);ymax2=max(y2);zmax2=max(z2);
xmin2=min(x2);ymin2=min(y2);zmin2=min(z2);

xmax=xmax1;
if xmax2 > xmax
   xmax=xmax2;
end
ymax=ymax1;
if ymax2 > ymax
   ymax=ymax2;
end
zmax=zmax1;
if zmax2 > zmax
   zmax=zmax2;
end
vmax=xmax;
if ymax > vmax
   vmax=ymax;
end
if zmax > vmax
   vmax=zmax;
end

xmin=xmin1;
if xmin2 < xmin
   xmin=xmin2;
end
ymin=ymin1;
if ymin2 < ymin
   ymin=ymin2;
end
zmin=zmin1;
if zmin2 < zmin
   zmin=zmin2;
end
vmin=xmin;
if ymin < vmin
   vmin=ymin;
end
if zmin < vmin
   vmin=zmin;
end

for i=1:nf1,
f1(i)=amp1_dat(i,1)/3.141592654/2;
xr1(i)=amp1_dat(i,2);
xi1(i)=amp1_dat(i,3);
xa1(i)=sqrt(xr1(i)*xr1(i)+xi1(i)*xi1(i));
yr1(i)=amp1_dat(i,4);
yi1(i)=amp1_dat(i,5);
ya1(i)=sqrt(yr1(i)*yr1(i)+yi1(i)*yi1(i));
zr1(i)=amp1_dat(i,6);
zi1(i)=amp1_dat(i,7);
za1(i)=sqrt(zr1(i)*zr1(i)+zi1(i)*zi1(i));
end
ax1=xa1(1);ay1=ya1(1);az1=za1(1);

for i=1:nf2,
f2(i)=amp_dat(i,1)/3.141592654/2;
xr2(i)=amp_dat(i,2);
xi2(i)=amp_dat(i,3);
xa2(i)=sqrt(xr2(i)*xr2(i)+xi2(i)*xi2(i));
yr2(i)=amp_dat(i,4);
yi2(i)=amp_dat(i,5);
ya2(i)=sqrt(yr2(i)*yr2(i)+yi2(i)*yi2(i));
zr2(i)=amp_dat(i,6);
zi2(i)=amp_dat(i,7);
za2(i)=sqrt(zr2(i)*zr2(i)+zi2(i)*zi2(i));
end
ax2=xa2(1);ay2=ya2(1);az2=za2(1);

dt1=t1(2)-t1(1);
dx1(1)=0.0;
dy1(1)=0.0;
dz1(1)=0.0;
for i=2:nt1,
dx1(i)=dx1(i-1)+(x1(i)+x1(i-1))*dt1/2;
dy1(i)=dy1(i-1)+(y1(i)+y1(i-1))*dt1/2;
dz1(i)=dz1(i-1)+(z1(i)+z1(i-1))*dt1/2;
end
dxmax1=max(dx1);dymax1=max(dy1);dzmax1=max(dz1); 
dxmin1=min(dx1);dymin1=min(dy1);dzmin1=min(dz1); 

dt2=t2(2)-t2(1);
dx2(1)=0.0;
dy2(1)=0.0;
dz2(1)=0.0;
for i=2:nt2,
dx2(i)=dx2(i-1)+(x2(i)+x2(i-1))*dt2/2;
dy2(i)=dy2(i-1)+(y2(i)+y2(i-1))*dt2/2;
dz2(i)=dz2(i-1)+(z2(i)+z2(i-1))*dt2/2;
end
dxmax2=max(dx2);dymax2=max(dy2);dzmax2=max(dz2); 
dxmin2=min(dx2);dymin2=min(dy2);dzmin2=min(dz2); 

dxmax=dxmax1;
if dxmax2 > dxmax
   dxmax=dxmax2;
end
dymax=dymax1;
if dymax2 > dymax
   dymax=dymax2;
end
dzmax=dzmax1;
if dzmax2 > dzmax
   dzmax=dzmax2;
end
dmax=dxmax;
if dymax > dmax
   dmax=dymax;
end
if dzmax > dmax
   dmax=dzmax;
end

dxmin=dxmin1;
if dxmin2 < dxmin
   dxmin=dxmin2;
end
dymin=dymin1;
if dymin2 < dymin
   dymin=dymin2;
end
dzmin=dzmin1;
if dzmin2 < dzmin
   dzmin=dzmin2;
end
dmin=dxmin;
if dymin < dmin
   dmin=dymin;
end
if dzmin < dmin
   dmin=dzmin;
end

dur=max(t1);
if max(t2) > dur
   dur=max(t2);
end
fmax=max(f1);
if max(f2) > fmax
   fmax=max(f2);
end

xamax=max(xa1);
if max(xa2) > xamax
   xamax=max(xa2);
end
yamax=max(ya1);
if max(ya2) > yamax
   yamax=max(ya2);
end
zamax=max(za1);
if max(za2) > zamax
   zamax=max(za2);
end

amax=xamax;
if yamax > amax
   amax=yamax;
end
if zamax > amax
   amax=zamax;
end

subplot(3,3,1)
plot(f1,xa1,f2,xa2,'--')
title(['NS: Fourier Amplitudes'])
xlabel('frequency (Hz)')
text(fmax*0.4,amax*0.95,['a1=',num2str(ax1)])
text(fmax*0.4,amax*0.85,['a2=',num2str(ax2)])
axis([0 fmax 0 amax*1.1])

subplot(3,3,2)
plot(f1,ya1,f2,ya2,'--')
title(['EW: Fourier Amplitudes'])
text(fmax*0.4,amax*0.95,['a1=',num2str(ay1)])
text(fmax*0.4,amax*0.85,['a2=',num2str(ay2)])
axis([0 fmax 0 amax*1.1])

subplot(3,3,3)
plot(f1,za1,f2,za2,'--')
title(['UD: Fourier Amplitudes'])
xlabel('frequency (Hz)')
text(fmax*0.4,amax*0.95,['a1=',num2str(az1)])
text(fmax*0.4,amax*0.85,['a2=',num2str(az2)])
axis([0 fmax 0 amax*1.1])

subplot(3,3,4)
plot(t1,x1,t2,x2,'--')
if xmax1 > -xmin1
  title(['max1=',num2str(xmax1),', max2=',num2str(xmax2)])
else
  title(['min1=',num2str(xmin1),', min2=',num2str(xmin2)])
end
ylabel('Velocity (cm/s)')
axis([0 dur vmin*1.1 vmax*1.1])

subplot(3,3,5)
plot(t1,y1,t2,y2,'--')
if ymax1 > -ymin1
  title(['max1=',num2str(ymax1),', max2=',num2str(ymax2)])
else
  title(['min1=',num2str(ymin1),', min2=',num2str(ymin2)])
end
ylabel('Velocity (cm/s)')
axis([0 dur vmin*1.1 vmax*1.1])

subplot(3,3,6)
plot(t1,z1,t2,z2,'--')
if zmax1 > -zmin1
  title(['max1=',num2str(zmax1),', max2=',num2str(zmax2)])
else
  title(['min1=',num2str(zmin1),', min2=',num2str(zmin2)])
end
ylabel('Velocity (cm/s)')
axis([0 dur vmin*1.1 vmax*1.1])

subplot(3,3,7)
plot(t1,dx1,t2,dx2,'--')
if dxmax1 > -dxmin1
  title(['max1=',num2str(dxmax1),', max2=',num2str(dxmax2)])
else
  title(['min1=',num2str(dxmin1),', min2=',num2str(dxmin2)])
end
ylabel('Displacement (cm)')
xlabel('time (sec)')
axis([0 dur dmin*1.1 dmax*1.1])

subplot(3,3,8)
plot(t1,dy1,t2,dy2,'--')
if dymax1 > -dymin1
  title(['max1=',num2str(dymax1),', max2=',num2str(dymax2)])
else
  title(['min1=',num2str(dymin1),', min2=',num2str(dymin2)])
end
ylabel('Displacement (cm)')
xlabel('time (sec)')
axis([0 dur dmin*1.1 dmax*1.1])

subplot(3,3,9)
plot(t1,dz1,t2,dz2,'--')
if dzmax1 > -dzmin1
  title(['max1=',num2str(dzmax1),', max2=',num2str(dzmax2)])
else
  title(['min1=',num2str(dzmin1),', min2=',num2str(dzmin2)])
end
ylabel('Displacement (cm)')
xlabel('time (sec)')
axis([0 dur dmin*1.1 dmax*1.1])

print waveVD2.ps
