%Arrays contain: [Depth; 2/20/98 v(cm/s); 2/21/98 v(cm/s)]
v_northSouth = csvread('1998_Feb_20_NorthSouth.csv',8,1);
v_northSouth = v_northSouth(:,1:51);
v_northSouth(1,:) = [5 v_northSouth(1, 1:50)];
v_northSouth(v_northSouth == -1e+34) = NaN;

v_eastWest = csvread('1998_Feb_20_EastWest.csv',8,1);
v_eastWest = v_eastWest(:,1:51);
v_eastWest(1,:) = [5 v_eastWest(1, 1:50)];
v_eastWest(v_eastWest == -1e+34) = NaN;

% Calculate velocity magnitudes and bearings
v_mag(1,:) = v_northSouth(1,:);
v_mag(2:3,:) = hypot(v_northSouth(2:3,:),v_eastWest(2:3,:));
v_bearing(1,:) = v_northSouth(1,:);
v_bearing(2,:) = radtodeg(pi/2 - atan2(v_northSouth(2,:), v_eastWest(2,:)));
v_bearing(3,:) = radtodeg(pi/2 - atan2(v_northSouth(3,:), v_eastWest(3,:)));
v_mag(4,:) = v_mag(1,:) / speed;
figure
plot(v_mag(4,:), v_mag(2,:))

%Plot velocity profiles
figure('Units', 'inches', 'Position', [2 4 10.5 6]);
% subplot(1,4,1);
% plot(v_northSouth(2,:), v_northSouth(1,:));
% xlabel('North-South Velocity (cm/s)');
% ylabel('Depth (m)');
% title('North-South Velocity Profile');
% subplot(1,4,2);
% plot(v_eastWest(2,:), v_eastWest(1,:));
% xlabel('East-West Velocity (cm/s)');
% %ylabel('Depth (m)');
% title('East-West Velocity Profile');
subplot(1,2,1)
plot(v_mag(2,:), v_mag(1,:));
xlabel('Velocity Magnitude (cm/s)');
ylabel('Depth (m)');
title('Velocity Magnitude Profile');
%set(gca, 'XAxisLocation', 'top');
set(gca, 'YDir', 'reverse');
subplot(1,2,2);
plot(v_bearing(2,:), v_bearing(1,:));
xlabel('Bearing (Degrees from north)');
ylabel('Depth (m)');
title('Current velocity bearing');
set(gca, 'YDir', 'reverse');

%Plot velocity vector field
figure(2)
q = quiver3(zeros(1,51), zeros(1,51), v_northSouth(1,:), v_eastWest(2,:), v_northSouth(2,:), zeros(1,51));
q.MaxHeadSize = 0.01;
title('Velocity Profile Quiver');

% Differentiate to get acceleration
speed = 0.1; %m/s
dist = 5; %m (Measurements taken every 5m)
dt = dist/speed;

accel_mag_2_20 = [diff(v_mag(2,:))/dt, NaN];
accel_mag_2_21 = [diff(v_mag(3,:))/dt, NaN];

accel_mag_2_20 = accel_mag_2_20 / 100 / 9.8; % Converting to g's
accel_mag_2_21 = accel_mag_2_21 / 100 / 9.8;

%Plot acceleration for both days
figure(3)
subplot(1,2,1);
plot(accel_mag_2_20, v_northSouth(1,:));
xlabel('Acceleration (g)');
ylabel('Depth (m)');
title('Acceleration for 2/20/1998');
subplot(1,2,2);
plot(accel_mag_2_21, v_northSouth(1,:));
xlabel('Acceleration (g)');
ylabel('Depth (m)');
title('Acceleration for 2/21/1998');