function plot_eddy(sled_start,sled_end,ref_time,rot_period)
% 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,ref_time,rot_rate)
%
% 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
%   ref_time - (string) time to reference plot to (format '03-Sep-2005 13:34:34')
%   rot_period - (float) rotation period in days
%
% 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

%%%%%%%%%%%%%%%% Load Eddy Center Trajectory Data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% center determined by subjective analysis of ADCP data display on Dennis' computer
[date, tme, latc1, lonc1] = textread('..\adcp\positions\trajectory.dat','%s %d %f %f',...
    'headerlines',1);
dtnum = datenum(date)+(tme/2400);


%%%%%%%%%%%%%%%%%% Load SF6 Concentration Data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load cex.mat

%%%%%%%%%%%%%%%%%%% Rotation rate (rad/s) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the rotation rate based on the user input of period

omega = (-2*pi)./(rot_period*86400); % 86,400 sec/day converts rot_period to seconds
                                     % negative sign indicates anit-cyclonic (clockwise) rotation 


%%%%%%%%%%%%%%% Loop to process each tow %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set an intial latitude, just used to scale rx below
init_lat = 29.846;

for ii = sled_start:1:sled_end,
    clear lonc latc
    % load the tow lat/lon data for Medusa
    eval(['load ed', num2str(ii), '_pos.mat'])
    % center lat/lon for the tow (the 0*m_st just makes the lonc variable the right length)
    lonc = lonc1(ii) + 0*m_st; % m_st comes from the ed#_pos.mat
    latc = latc1(ii) + 0*m_st;
    % relative distance from center in degrees
    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

    % Eddy A4 has solid body rotation out to 20 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;
    
    % Produce polar plot of tow results
    figure(1)
    h1 = polarTrueNorth2(theta_rot,radius);
    set(h1,'color','k')
    hold on
    jet_c = jet;
    for aa = 1:length(theta_rot)
        h2 = polarTrueNorth2(theta_rot(aa),radius(aa));
        if ~isnan(cmat(aa,ii))
            jet_row = floor((cmat(aa,ii)/500)*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
    %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
    
    % Plot of tows in lat/lon coordinates
    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

% add title and colorbar to polar plot
figure(1)
eval(['title(''Tow Tracks Referenced to ', ref_time, ' UTC with Eddy Rotation Period = ', ...
    num2str(rot_period), ' days'')'])
cbar_axes = colorbar;
set(cbar_axes,'YTick',[0 0.2 0.4 0.6 0.8 1],'YTickLabel',{'0','100','200','300','400','500'})


