% azimuthalv.m
% Find the azimuthal velocity as a fxn of distance from center
%
hold off

% select the observations to be studied
n1 = 136;
n2 = 360;

% locate the center
% for the purposes of this exercise, we will assume it to be:
lonc = -66.97;
latc = 30.52;

% express coordinates relative to center
dlon = lon(n1:n2) - lonc;
dlat = lat(n1:n2) - latc;

% express all velocities in polar coordinates wrt center
radius = sqrt((dlat).^2 + (dlon).^2);
theta = atan2(dlat,dlon);

% calculate the azimuthal projection of the velocities
% find the angle to horizontal (E-W) of the velocity
alpha = atan2(vm(n1:n2),um(n1:n2));
% From alpha & theta, find the angle of the projection
% onto azimuthal
phi = -pi/2 + theta - alpha;
%phi = pi/2 - (theta + alpha);
%phi = mod(phi, 2*pi)-pi;
% find the total magnitude of the velocity
speed = sqrt( um(n1:n2).^2 + vm(n1:n2).^2 );
v_perp = speed.*cos(phi);

% Take out bad points (from when ADCP was down)
ind = find(speed > 1);
v_perp(ind) = NaN;

% bin the velocities into histogram of radius
% declare a vector to store the average velocities
avgvel = zeros(9,1);
sig = zeros(9,1);
rad = [0.025:0.05:0.425];
% convert radial distance to meters
rad = rad.*(60*1600);
for ii = 1:9
index = find(radius < ii*.05 & radius >= (ii-1)*.05);
avgvel(ii) = nanmean(v_perp(index));
sig(ii) = nanstd(v_perp(index));
end

hold on
% Convert to meters
radius_m = radius*60*1600;
plot(radius_m,v_perp, '.')
xlabel('Distance from Center of Eddy (m)')
ylabel('Azimuthal Velocity (m/s)')
plot(rad, avgvel,'r')
errorbar(rad, avgvel, sig,'r')
axis([0 35000 -0.1 0.5])
title({'Rotational Velocity of Eddy A4 at 100m Depth';'OC415-2 29-31 Jul 2005'})

% Do a linear best fit of the velocities
% Remove all the NaN values
no_good1 = isnan(radius_m);
no_good2 = isnan(v_perp);
all_no_good = no_good1 | no_good2;
real_num = find(all_no_good == 0);
radius_num = radius_m(real_num);
v_perp_num = v_perp(real_num);
coeff = polyfit(radius_num, v_perp_num, 1);
xx = [0:1:40000];
yy = coeff(1)*xx+coeff(2);
plot(xx, yy, 'g', 'LineWidth', 3)
% Annotate figure with rotation rate from best fit
% Convert slope to units of days/cycle
period = 2*pi*(coeff(1)*86400)^(-1);
text(25000,0,['Period of Eddy Rotation: ',num2str(period)])

hold off