%Tadd Truscott modified 3-31-05
% Jason Dahl
% Find Values and Make plots for 13.42 Lab 3
% 5/5/03
% Use this file to anlyze the .asc runs from the 13.42 lab here you will
% get heave amplitude, wave frequency, wave amplitude, (heave amplitude)/(wave
% amplitude), pitch response amplitude, and the phase angle.

%%$$$$$$ PLEASE DON"T FORGET TO INPUT THE DISTANCE FROM THE SHIP TO THE WAVEPROBE!

%Enter the calibration coefficients supplied by you friendly TA for the
%heave and pitch sensors, as well as the wave probe calibration that you
%obtained, in the appropriate place below:
clear all;
%----------------------------------------------------------------------
%----------------------------------------------------------------------
%Enter calibration data here

heave_coefficient = .7094;    %in/volt
heave_offset =     1.06937;  %volt
pitch_coefficient =  15.135;   %rad/volt
pitch_offset =  0.009599311;   %volt
wave_probe_coefficient = -1.2719; %in/volt
wave_probe_offset =  -0.113094;    %volt

%----------------------------------------------------------------------
%----------------------------------------------------------------------

name = input('What is .asc filename in directory:  ','s');

%input asci file into a matrix. rename columns. For this lab columns are
%arranged as follows.
%time, pitch, heave, upstream waveprobe, midship waveprobe.
%and apply calibration coefficient and offset for all sensors

M=dlmread(name,'\t',8,0);
time=M(:,1);
heave_position=(M(:,3)+heave_offset)*heave_coefficient;
surface_height=(M(:,5)+wave_probe_offset)*wave_probe_coefficient; %midships wave probe.
pitch_angle=(M(:,2)+pitch_offset)*pitch_coefficient;

figure(1);clf;hold off;
plot(time,heave_position);
xlabel('time, s');
ylabel('heave, in');
title('Pick peak to peak amplitude of heave response: y1-y2');
zoom on; pause;
[junk,y] = ginput(2);
heave_amp = y(1)-y(2);
disp(sprintf('Heave amplitude: %g in',heave_amp));

figure(2);clf;hold off;
plot(time,surface_height);
xlabel('time, s');
ylabel('surface height, in');
title('Pick points (one wavelength) to determine frequency: 1/(t2-t1)');
zoom on; pause;
[x, junk2] = ginput(2);
omega = 1/(x(2)-x(1));
disp(sprintf('Wave Frequency: %g hz',omega));

figure(3);clf;hold off;
plot(time,surface_height);
xlabel('time, s');
ylabel('surface height, in');
title('Pick peak to peak amplitude of wave surface height: y1-y2');
zoom on; pause;
[junk,y] = ginput(2);
eta = y(1)-y(2);
disp(sprintf('Wave amplitude: %g in',eta));
disp(sprintf('Heave/Amp: %g',heave_amp./eta));

figure(4);clf;hold off;
plot(time,pitch_angle);
xlabel('time, s');
ylabel('pitch, rad');
title('Pick peak to peak amplitude of pitch response: y1-y2');
zoom on; pause;
[junk,y] = ginput(2);
pitch = y(1)-y(2);
disp(sprintf('Pitch Response: %g rad',pitch));

%This plot will help find the phase shift
%figure(5);clf;hold off;
%plot(time,heave_position, time, surface_height,'--');
%xlabel('time, s');
%ylabel('heave, eta (in)');
%title('Phase Shift - 0.85 Hz');
%legend('heave','wave surface');
%zoom on; pause;
%[x,junk] = ginput(2);
%t = x(2)-x(1);
%phi0 = 2*pi*omega*t;
%k = (2*pi*omega)^2/9.81;
%xm = 0;  % m, distance between ship and wave probe. if using midship this is zero
%phi1 = phi0-k*xm;
%disp(sprintf('Phase shift: %g rad', phi1));
%disp(sprintf('Phase shift: %g deg',phi1*180/pi));=