% Algorithm to locate the center of the eddy based on 
% a track of current velocity measurements and a little 
% pointing and clicking


%load data (u,v,x,y)
load Z:\oc415_2\adcp\contour_uv
load Z:\oc415_2\adcp\contour_xy

% Depth-average all of the speeds over the 22 depth bins
% Each bin is 10 m deep
nbins=22;
u=uv(1:nbins,1:2:end);
v=uv(1:nbins,2:2:end);
% take the mean of all the depth bins, ignoring NaNs
um=nanmean(u);
vm=nanmean(v);

% Convert x,y,t to latitude and longitude
% time is in Julian Days 
lon=xyt(1,:)-360;
lat=xyt(2,:);
time=xyt(3,:);

% indicate the beginning and end points of the track in question
n1 = 231;
n2 = 245;

% Make array containing the beginning and end of each straight track
limit = [136 155
         158 175
         176 188
         189 205
         231 245
         260 372
         386 412
         543 560
         645 660
         718 738
         739 744
         745 787
         788 797
         798 817
         819 841];

% Calculate slope perpendicular to the local velocity & the angle to the
% horizontal
m_perp = -um./vm;
theta = atan(m_perp);
intercept = lat - m_perp.*lon;

% % plot all the perpendicular lines
% hold on
% for i=n1:1:n2
%     x = [-0.5:.01:0.5];
%     x = x + lon(i);
%     y = m_perp(i)*x+intercept(i);
%     plot(x,y,'--r');
% end
% hold off
    
% calculate magnitude of each velocity vector
speed = sqrt(um.^2 + vm.^2);

% Associated with each speed is a range of probable distances to the 
% center of the eddy, assuming the eddy rotates once every 5-7 days
T_min = 5; %Minimum rotational period in days
T_max = 7; 
rmin_m = speed*2*pi*T_min*86400;
rmax_m = speed*2*pi*T_max*86400;
% NB: distances in units of meters
% %rmin_m  = speed/(3.6*10^(-5));
% %rmax_m = speed/(1.5*10^(-5));
% % Use radii such that the eddy rotates every 5-7 days
% rmin_m = speed/(1.5*10^(-5));
% rmax_m = speed/(1.0*10^(-5));


% Convert into units of latitude and longitude
rmin  = rmin_m/(1600*60);
rmax = rmax_m/(1600*60);

% Find x-coordinates that correspond to the desired radii
% Calculate x-projection of a vector of length rmin along 
% the perpendicular to the velocity
xin = rmin.*cos(theta);
xout = rmax.*cos(theta);

xdiff = abs(xout - xin);

hold on
for ii=n1:n2
    x_neg(ii,:) = [lon(ii)-xout(ii):xdiff(ii)/50:lon(ii)-xin(ii)];
    y_neg(ii,:) = m_perp(ii)*x_neg(ii,:)+intercept(ii);
    plot(x_neg(ii,:),y_neg(ii,:), 'r', 'LineWidth', 1)

    x_pos(ii,:) = [lon(ii)+xin(ii):xdiff(ii)/50:lon(ii)+xout(ii)];
    y_pos(ii,:) = m_perp(ii)*x_pos(ii,:)+intercept(ii);
    plot(x_pos(ii,:),y_pos(ii,:), 'g', 'LineWidth', 1)        
end
hold off

% Take the center point by clicking on the figure where all the 
% little red line segments intersect:
[x,y] = ginput(1)


