function plot_eddy(sled_start,sled_end,drogue2use,ref_time)
% plot_eddy  plot sled tow position in rotating eddy coordinates
% =========================================================================
% plot_eddy  Version 1.0 02-Sep-2005
%
% Usage:
%   plot_eddy(sled_start,sled_end,drogue2use,ref_time)
%
% Description:
%   This function loads position of the sampler sled tows and plots this
%   information on a polar plot which represents the rotating eddy coordinates.
%
%
% Input:
%   sled_start - (int) number of the first sled tow you want plotted
%   sled_end - (int) number of the last sled tow you want plotted
%   drogue2use - (string) which drogue results to use ('52484','52486','all')
%   ref_time - (string) time to reference plot to (format '03-Sep-2005 13:34:34')
%
% Output:
%   Figure only
%
% Called by:
%   n/a
%
% Calls:
%   n/a
%
% Author:
%   Blair Greenan
%   Bedford Institute of Oceanography
%   Matlab 7
%   greenanb@mar.dfo-mpo.gc.ca
% =========================================================================
%

% History
% Version 1.0 02-Sep-2005

%%%%%%%%%%%%%%%% RESULTS FROM FIT OF DROGUE DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set a reference year day that represents the start of the fit to the drogue data
%ref_yrday = 238.0;
%ref_yrday = 240.1;

switch lower(drogue2use)  % convert string to lower case for the switch statement
    case '52484'
      ref_yrday = 240.0;
        % Initial latitude and longitude of the eddy center
        %        init_lon = -69.05;
        %        init_lat = 30.02;
        init_lon = -68.931;
        init_lat = 29.851;
        % translation velocities (m/s)
        %        u = -0.038; % positive to the east
        %        v = -0.033; % positive to the north
        u = -0.049; % positive to the east
        v = -0.0024; % positive to the north
        % rotation rate (rad/s)
        %        omega = -1.216e-5; % 5.98 days to complete one rotation
        omega = -1.0927e-5; % 6.66 days to complete one rotation
    case '52486'
      ref_yrday = 240.0;
        % Initial latitude and longitude of the eddy center
        %        init_lon = -69.01;
        %        init_lat = 30.04;
        init_lon = -68.970;
        init_lat = 29.846;
        % translation velocities (m/s)
        %        u = -0.028; % positive to the east
        %        v = -0.046; % positive to the north
        u = -0.068; % positive to the east
        v = -0.017; % positive to the north
        % rotation rate (rad/s)
        %omega = -1.1507e-5; % 6.32 days to complete one rotation % based on James's drifter analysis
        %rot_period = 6.32;
        %omega = -1.2120e-5; % 6.00 days to complete one rotation % based on James's drifter analysis
        %rot_period = 6.00;
        omega = -1.0388e-5; % 7.00 days to complete one rotation % based on James's drifter analysis
        rot_period = 7.00;
        %omega = -1.206e-5; % 6.03 days to complete one rotation
        %omega = -1.016e-5; % 7.16 days to complete one rotation
        %omega = -1.423e-5; % 5.11 days to complete one rotation based on ADCP fit
        %    case '52486'
        % Initial latitude and longitude of the eddy center
        %        init_lon = -68.97;
        %        init_lat = 29.846;
        % translation velocities (m/s)
        %        u = -0.0675; % positive to the east
        %        v = -0.0166; % positive to the north
        % rotation rate (rad/s)
        %        omega = -1.016e-5; % 7.16 days to complete one rotation
    case 'all'
      ref_yrday = 240.0;
        % Initial latitude and longitude of the eddy center
        init_lon = (-69.05-69.01)./2 ;
        init_lat = (30.02+30.04)./2;
        % translation velocities (m/s)
        u = (-0.038-0.028)./2; % positive to the east
        v = (-0.033-0.046)./2; % positive to the north
        % rotation rate (rad/s)
        omega = (-1.216e-5-1.206e-5)./2; % 6.01 days to complete one rotation
end


%%%%%%%%%%%%%%%

for ii = sled_start:1:sled_end,
    eval(['load ed', num2str(ii), '_pos.mat'])

    % Fix eddy position for each tow
    switch ii,
        case 1,
            lonc = -69.20 + 0*m_st; % m_st comes from the ed#_pos.mat
            latc = 29.84 + 0*m_st;
        case 2,
            lonc = -69.235 + 0*m_st;
            latc = 29.785 + 0*m_st;
        case 3,
            lonc = -69.26 + 0*m_st;
            latc = 29.79 + 0*m_st;
        case 4,
            lonc = -69.34 + 0*m_st;
            latc = 29.72 + 0*m_st;
        case 5, % interpolatede from Tow 4 and Tow 6 since we did not start the tow near center
            lonc = -69.39 + 0*m_st;
            latc = 29.71 + 0*m_st;
        case 6,
            lonc = -69.44 + 0*m_st;
            latc = 29.70 + 0*m_st;
        case 7,
            lonc = -69.47 + 0*m_st;
            latc = 29.68 + 0*m_st;
        case 8,
            lonc = -69.50 + 0*m_st;
            latc = 29.675 + 0*m_st;
        otherwise
            % Use extrapolated center from drifter model
            lonc = init_lon + (u*(m_st-ref_yrday)*86400)/1.11e5/cos(init_lat*pi/180); % longitude of eddy center
            latc = init_lat + (v*(m_st-ref_yrday)*86400)/1.11e5; % latitude of eddy center
    end

    if 0
        mlon = interp1([m_st(1) m_st(end)],[mlon(1) mlon(end)],m_st)';
        mlat = interp1([m_st(1) m_st(end)],[mlat(1) mlat(end)],m_st)';
        clf;plot(lonc,latc,'.-');grid on;
        hold on; plot(init_lon,init_lat,'rp');
        plot(mlon,mlat,'m.-');
        pause
    end

    rlon = mlon' - lonc;
    rlat = mlat' - latc;
    % convert distance from degrees to m
    rx = rlon*cos(init_lat*pi/180)*1.11e5;
    ry = rlat*1.11e5;
    radius = sqrt(rx.^2 + ry.^2);
    theta = atan2(ry,rx);
    % compute the current year day
    frstday=datenum(2005,1,1,0,0,0);
    %nowdate=now -0;
    nowdate=datenum(ref_time);
    netda=nowdate-frstday;
    % rotate the angle to the current time

    %clf; plot(netda-m_st,omega.*((netda-m_st)*86400),'.-'); grid on; pause
    % Eddy A4 has solid body rotation out to 16 km radius based on ADCP results
    all_solid = 0;
    if all_solid
        theta_rot = theta + (omega.*((netda-m_st)*86400));  % convert to seconds and then multiply by omega to get rad/s * s
    else
        solid_index = find(radius <= 20000); % 20 km breakpoint in solid body rotation
        theta_rot(solid_index) = theta(solid_index) + (omega.*((netda-m_st(solid_index))*86400));
        unsolid_index = find(radius > 20000);
        yy1 = -3.88e-3 + (1.2928e-5*radius(unsolid_index));
        yy2 = 2.81e-1 - (1.3197e-6*radius(unsolid_index));
        omega_ratio = yy2./yy1;
        theta_rot(unsolid_index) = theta(unsolid_index) + ((omega.*omega_ratio).*((netda-m_st(unsolid_index))*86400));
        theta_rot = theta_rot';
    end
    %theta_rot = theta;
    
    % load the carousel data from the sled
    if ii ~= 8
        if ii < 10
            eval(['load Car_0', num2str(ii), '.txt'])
            eval(['SF6_area = Car_0', num2str(ii), '(:,6);'])
        else
            eval(['load Car_', num2str(ii), '.txt'])
            eval(['SF6_area = Car_', num2str(ii), '(:,6);'])
        end
    end

    figure(1)
    h1 = polarTrueNorth2(theta_rot,radius);
    set(h1,'color','k')
    hold on
    if ii ~= 8
        jet_c = jet;
        for aa = 1:length(theta_rot)
            h2 = polarTrueNorth2(theta_rot(aa),radius(aa));
            if ~isnan(SF6_area(aa))
                jet_row = floor((SF6_area(aa)/1000000)*64);
                if jet_row == 0
                    jet_row = 1;
                end
                if jet_row > 64
                    jet_row = 64;
                end
                jet_color = jet_c(jet_row,:);
                set(h2,'MarkerSize',3,'Marker','o','MarkerEdgeColor',jet_color,'MarkerFaceColor',jet_color);
                clear h2
            end
        end
    end
    %h1 = polar(theta_rot,radius);
    xx = radius(end).*cos(theta_rot(end));
    yy = radius(end).*sin(theta_rot(end));
    tow = ['Tow ' num2str(ii)];
    text(xx,yy,tow);
    hold on
    figure(2)
    h2 = plot(mlon,mlat);
    tow = ['Tow ' num2str(ii)];
    text(mlon(end),mlat(end),tow);
    xlabel('Longitude')
    ylabel('Latitude')
    grid on
    hold on
    plot(lonc,latc,'rd','markerfacecolor','r')
    set(gca,'dataaspectratio',[1 cos(init_lat*pi/180) 1]);
    clear theta_rot
end

figure(1)
eval(['title(''Tow Tracks Referenced to 08-Sep-2005 2300 UTC with Eddy Rotation Period = ', ...
    num2str(rot_period), ' days'')'])


