% 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 wave_dat;
load amp_dat;

load wave1_dat;
load amp1_dat;

nt=size(wave_dat,1);
nf=size(amp_dat,1);
nt1=size(wave1_dat,1);
nf1=size(amp1_dat,1);

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

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

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:nf,
f(i)=amp_dat(i,1)/3.141592654/2;
xr(i)=amp_dat(i,2);
xi(i)=amp_dat(i,3);
xa(i)=sqrt(xr(i)*xr(i)+xi(i)*xi(i));
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));
end

dt=t(2)-t(1);
dx(1)=0.0;
for i=2:nt,
dx(i)=dx(i-1)+(x(i)+x(i-1))*dt/2;
end
dmax=max(dx); 
dmin=min(dx); 

dt1=t1(2)-t1(1);
dx1(1)=0.0;
for i=2:nt1,
dx1(i)=dx1(i-1)+(x1(i)+x1(i-1))*dt1/2;
end
dmax1=max(dx1); 

subplot(3,1,1)
plot(t1,dx1,t,dx,'--')
title(['NS (max1=',num2str(dmax1),', max2=',num2str(dmax),')'])
ylabel('Displacement (cm)')
xlabel('time (sec)')
%axis([0 dur -amax amax])

subplot(3,1,2)
plot(t1,x1,t,x,'--')
title(['NS (max=',num2str(xmax1),', min=',num2str(xmin1),')'])
ylabel('Velocity (cm/s)')
%axis([0 dur -amax amax])

subplot(3,1,3)
plot(f1,xa1,f,xa,'--')
title(['Fourier Amplitudes (NS)'])
xlabel('frequency (Hz)')
%axis([0 dur -amax amax])

print waveVD2.ps
