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

% select the observations to be studied (OC415_4)
n1 = 12440;
n2 = 14417;

% load the adcp velocity data
load contour_uv.mat
load contour_xy.mat
%nd=22;
u=uv(6:10,1:2:end);
v=uv(6:10,2:2:end);
%u=uv(1:22,1:2:end);
%v=uv(1:22,2:2:end);
um=nanmean(u);
um = um(n1:n2);
vm=nanmean(v);
vm = vm(n1:n2);
% get rid of outliers
find_outu = find(abs(um)>1);
um(find_outu) = NaN;
find_outv = find(abs(vm)>1);
vm(find_outv) = NaN;

lon=xyt(1,n1:n2)-360;
lat=xyt(2,n1:n2);
time=xyt(3,n1:n2);
frstday=datenum(2005,1,1,0,0,0);
time = time + frstday;

plot(lon,lat,'m-','linewidth',2);
h2=vecplot2(lon,lat,um,vm,1.0,'cm/sec');

% center determined by subjective analysis of ADCP data display on Dennis' computer
[date, tme, latc, lonc] = textread('P:\oc415_4\adcp\positions\trajectory.dat','%s %d %f %f',...
    'headerlines',1);
dtnum = datenum(date)+(tme/2400);
figure(10)
plot(dtnum,latc,'.-')
figure(11)
plot(dtnum,lonc,'.-')

% fit the center lat and lon to straight lines
plat = polyfit(dtnum,latc,1);
plon = polyfit(dtnum,lonc,1);
ylat = polyval(plat,dtnum(1):0.05:dtnum(end));
ylon = polyval(plon,dtnum(1):0.05:dtnum(end));
figure(10)
hold on
plot(dtnum(1):0.05:dtnum(end),ylat,'r.-')
figure(11)
hold on
plot(dtnum(1):0.05:dtnum(end),ylon,'r.-')

% determine a center point for each 5-min ADCP velocity sample
latc = plat(2) + plat(1)*time;
lonc = plon(2) + plon(1)*time;

% express coordinates relative to center in meters
%dlon = lon - lonc;
%dlat = lat - latc;
for bb = 1:length(lon)
    [dlon(bb),dlat(bb)] = dist2(lat(bb),lon(bb),latc(bb),lonc(bb));
end
dlon = 1000.*dlon;
dlat = 1000.*dlat;

% 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,um);
  % 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.^2 + vm.^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(20,1);
sig = zeros(20,1);
rad = [3000:6000:120000];
% convert radial distance to meters
%rad = rad.*(60*1600);
for ii = 1:20
    index = find((radius < ((ii*6000)-1)) & (radius >= (6000*(ii-1))));
    avgvel(ii) = nanmean(v_perp(index));
    sig(ii) = nanstd(v_perp(index));
end

figure
hold on
% Convert to meters
%radius_m = radius*60*1600;
radius_m = radius;
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';'OC415-4 30 Aug-5 Sep 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);

% break the results into two groups (Group1 <= 20km; Group2 >20km)
brkpoint = 20000;
group1_index = find(radius_num <= brkpoint);
group2_index = find(radius_num > brkpoint);
radius_group1 = radius_num(group1_index);
radius_group2 = radius_num(group2_index);
v_perp_group1 = v_perp_num(group1_index);
v_perp_group2 = v_perp_num(group2_index);

coeff_group1 = polyfit(radius_group1, v_perp_group1, 1)
coeff_group2 = polyfit(radius_group2, v_perp_group2, 1)
xx1 = [0:1000:brkpoint];
yy1 = coeff_group1(1)*xx1+coeff_group1(2);
plot(xx1, yy1, 'g', 'LineWidth', 2)
xx2 = [brkpoint:1000:120000];
yy2 = coeff_group2(1)*xx2+coeff_group2(2);
plot(xx2, yy2, 'g', 'LineWidth', 2)


%%%%%%%%%% Compute fraction of omega for points beyond 16 km %%%%%%%%%%%%%%%%%%%%%%%%%%
figure
xx1 = [0:1000:120000];
yy1 = coeff_group1(1)*xx1+coeff_group1(2);
xx2 = [0:1000:120000];
yy2 = coeff_group2(1)*xx2+coeff_group2(2);
plot(xx2(17:end), yy2(17:end)./yy1(17:end), 'g.-', 'LineWidth', 2)


%{
% 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)])

%}
