%wavedata.m 
%Program to determine wave characteristics 
%author - Joshua Davis 
%MIT Towing Tank 
%Nov 2000 
 
hold off; zoom off; clear all; close all; 
 
%cal3;

%Load calibration data 
load cal; 
 
%open window to get  data file: 
[filen, pathn] = uigetfile('*.i32','OPEN Wave Data File'); 
filestub = [pathn '' filen]; 
 
%Load i32 format data and create array with channel 1,2,3 in columns: 
[dataall]= parsedl2(filestub,3)'; 
 
%Apply calibration coefficient and offset for wave probe 3 
data=(dataall(:,3)-zero_off(3))*ccoefs(3); 
 
%Develop time scale 
dt=.01; %Sampling at 100Hz 
time=0:dt:(length(data)-1)*dt; 
 
%Get region of interest and shorten data 
plot(time,data); 
title('Pick region of wave action (two horizontal points)'); 
[a,b]=ginput(2); 
[a]=round(a)/dt; 
data=data(a(1):a(2)); 
time=time(a(1):a(2)); 
 
 
%Plot data and calculate maximum wave displacement 
plot(time,data); 
title('Pick peak-to-peak wave amplitude (two vertical points)'); 
xlabel('Time (s)'); 
ylabel('Wave Elevation (cm)'); 
[x,y]=ginput(2); 
eta=abs(diff(y))/2; 
 
%Perform FFT to get peak wave frequency 
dataf=fft(data); 
n = a(2) - a(1) + 1 ; 
w_nyq = pi/dt ;    
w_scale = ([0:n/2-1])/(n/2) * w_nyq ;    % freq axis spacing 
 
figure(1);clf ;hold off ;clf; 
 
lw=2; 
upper_freq = 30 ; 
plot(w_scale,abs([dataf(1:round(n/2-1))])*2/n,'b','LineWidth',lw) ;
a = axis ; 
axis([0 upper_freq 0 a(4)]) ; 
grid; 
ylabel('spectral density'); 
title('wave frequency spectra: pick largest peak'); 
xlabel('frequency, rad/s'); 
 
[xt,yt] = ginput(1) ; 
 
a = axis ;
% (aspect ratio is about 0.66)
[junk,index] = min( ( (abs(dataf(1:round(n/2-1)))*2/n-yt) / (a(4)-a(3)) * .666).^2 +... 
		   ((w_scale'-xt)/upper_freq).^2           ) ; 
hold on; 
plot(w_scale(index),abs(dataf(index))*2/n,'ro','LineWidth',lw); 
hold off; 
pause(1) ; 

f_peak = w_scale(index)/(2*pi)  ; 
period=1/f_peak; 
g=9.81; 
 
dwwl=(g*period^2)/(2*pi); 
 
 
%All wave probes 
wave1=(dataall(:,1)-zero_off(1))*ccoefs(1); 
wave2=(dataall(:,2)-zero_off(2))*ccoefs(2); 
wave3=(dataall(:,3)-zero_off(3))*ccoefs(3); 
 
%Get phase velocity by plotting closely spaced probes 
time=0:dt:(length(wave1)-1)*dt; 
plot(time,wave2); 
hold on; 
plot(time,wave1,'g'); 
title('Zoom in on a few crests. Then press a key to continue.'); 
zoom on; 
pause; 
title('Pick green crest [Probe 1] then pick following blue crest [Probe 2]'); 
xlabel('Time (s)'); 
ylabel('Wave Elevation (cm)'); 
[x,y]=ginput(2); 
delay=diff(x); 
smalldist=12*2.54; %13.021 Wave Lab Nov 2000 
vp=smalldist/delay; 
hold off; 
 
%Get group velocity 
plot(time,wave1,'g'); 
hold on; 
plot(time,wave3,'b'); 
title('Zoom to show beginning of wave packet at both probes. Then press a key to continue.'); 
zoom on; 
pause; 
title('Pick a green crest [Probe 1] then pick the same blue crest [Probe 3]'); 
[x,y]=ginput(2); 
delay=diff(x); 
bigdist=30*12*2.54; 
vg=bigdist/delay; 
 
%Display data 
disp(sprintf('------------RESULTS-----------------')); 
disp(sprintf(' Wave Amplitude (eta):  %8.2f cm', eta)); 
disp(sprintf('            Frequency:  %8.2f Hz   (%3.2f rad/s)', f_peak,f_peak*2*pi)); 
disp(sprintf('Deep-Water Wavelength:  %8.2f m', dwwl)); 
disp(sprintf('       Phase Velocity:  %8.2f cm/s', vp)); 
disp(sprintf('       Group Velocity:  %8.2f cm/s', vg));