% Plot the injection and sampling tracks in eddy (rotating) coordinates
%
hold off
clf

load Z:\oc415_2\adcp\Rebecca\tow_tracks.mat
load_med
load Z:\oc415_2\adcp\Rebecca\med_ind.mat

% Define eddy center as a function of time (given in year day)

% For the purposes of this exercise, we will assume the center
% be stationary during the 18 hour period of injection, at:
% (30.28, -67.29)
% We further assume that it is moving along the following track
% during sampling:
% y_c = -0.02*time + 34.285;
% x_c = -0.0325*time - 60.7619;

% express coordinates relative to center
lonc_i3 = -67.29;
latc_i3 = 30.28;
rlon_i3 = lon_i3 - lonc_i3;
rlat_i3 = lat_i3 - latc_i3;
rx_i3 = rlon_i3*0.866*1.11*10^2;
ry_i3 = rlat_i3*1.11*10^2;

lonc_i4 = -67.29;
latc_i4 = 30.28;
rlon_i4 = lon_i4 - lonc_i4;
rlat_i4 = lat_i4 - latc_i4;
rx_i4 = rlon_i4*0.866*1.11*10^2;
ry_i4 = rlat_i4*1.11*10^2;

lonc_s1 = -0.0325*day_s1 - 60.7619;
latc_s1 = -0.02*day_s1 + 34.285;
rlon_s1 = lon_s1 - lonc_s1;
rlat_s1 = lat_s1 - latc_s1;
rx_s1 = rlon_s1*0.866*1.11*10^2;
ry_s1 = rlat_s1*1.11*10^2;

lonc_s2 = -0.0325*day_s2 - 60.7619;
latc_s2 = -0.02*day_s2 + 34.285;
rlon_s2 = lon_s2 - lonc_s2;
rlat_s2 = lat_s2 - latc_s2;
rx_s2 = rlon_s2*0.866*1.11*10^2;
ry_s2 = rlat_s2*1.11*10^2;

lonc_s3 = -0.0325*day_s3 - 60.7619;
latc_s3 = -0.02*day_s3 + 34.285;
rlon_s3 = lon_s3 - lonc_s3;
rlat_s3 = lat_s3 - latc_s3;
rx_s3 = rlon_s3*0.866*1.11*10^2;
ry_s3 = rlat_s3*1.11*10^2;

lonc_s4 = -0.0325*day_s4 - 60.7619;
latc_s4 = -0.02*day_s4 + 34.285;
rlon_s4 = lon_s4 - lonc_s4;
rlat_s4 = lat_s4 - latc_s4;
rx_s4 = rlon_s4*0.866*1.11*10^2;
ry_s4 = rlat_s4*1.11*10^2;

lonc_s5 = -0.0325*day_s5 - 60.7619;
latc_s5 = -0.02*day_s5 + 34.285;
rlon_s5 = lon_s5 - lonc_s5;
rlat_s5 = lat_s5 - latc_s5;
rx_s5 = rlon_s5*0.866*1.11*10^2;
ry_s5 = rlat_s5*1.11*10^2;

lonc_s6 = -0.0325*day_s6 - 60.7619;
latc_s6 = -0.02*day_s6 + 34.285;
rlon_s6 = lon_s6 - lonc_s6;
rlat_s6 = lat_s6 - latc_s6;
rx_s6 = rlon_s6*0.866*1.11*10^2;
ry_s6 = rlat_s6*1.11*10^2;

lonc_s7 = -0.0325*day_s7 - 60.7619;
latc_s7 = -0.02*day_s7 + 34.285;
rlon_s7 = lon_s7 - lonc_s7;
rlat_s7 = lat_s7 - latc_s7;
rx_s7 = rlon_s7*0.866*1.11*10^2;
ry_s7 = rlat_s7*1.11*10^2;

% express all velocities in polar coordinates wrt center
% then convert radius from lat/lon units to meters
%radius_i3 = sqrt((rlat_i3).^2 + (rlon_i3*0.0866).^2);
radius_i3 = sqrt((rx_i3).^2 + (ry_i3).^2);
theta_i3 = atan2(rx_i3,ry_i3);

%radius_i4 = sqrt((rlat_i4).^2 + (rlon_i4*0.0866).^2);
radius_i4 = sqrt((rx_i4).^2 + (ry_i4).^2);
theta_i4 = atan2(rx_i4,ry_i4);

%radius_s1 = sqrt((rlat_s1).^2 + (rlon_s1*0.0866).^2);
radius_s1 = sqrt((rx_s1).^2 + (ry_s1).^2);
theta_s1 = atan2(rx_s1,ry_s1);

%radius_s2 = sqrt((rlat_s2).^2 + (rlon_s2*0.0866).^2);
radius_s2 = sqrt((rx_s2).^2 + (ry_s2).^2);
theta_s2 = atan2(rx_s2,ry_s2);

%radius_s3 = sqrt((rlat_s3).^2 + (rlon_s3*0.0866).^2);
radius_s3 = sqrt((rx_s3).^2 + (ry_s3).^2);
theta_s3 = atan2(rx_s3,ry_s3);
% mradius_s3 = sqrt((mx_s3).^2 + (my_s3).^2);
% mtheta_s3 = atan2(mx_s3, my_s3);

%radius_s4 = sqrt((rlat_s4).^2 + (rlon_s4*0.0866).^2);
radius_s4 = sqrt((rx_s4).^2 + (ry_s4).^2);
theta_s4 = atan2(rx_s4,ry_s4);

%radius_s5 = sqrt((rlat_s5).^2 + (rlon_s5*0.0866).^2);
radius_s5 = sqrt((rx_s5).^2 + (ry_s5).^2);
theta_s5 = atan2(rx_s5,ry_s5);

%radius_s6 = sqrt((rlat_s6).^2 + (rlon_s6*0.0866).^2);
radius_s6 = sqrt((rx_s6).^2 + (ry_s6).^2);
theta_s6 = atan2(rx_s6,ry_s6);

radius_s7 = sqrt((rx_s7).^2 + (ry_s7).^2);
theta_s7 = atan2(rx_s7,ry_s7);

% Correct theta for the rotation of the eddy, where we assume
% that the eddy started rotating on yearday 0, and has been
% rotating at a speed of 6.5 days per rotation ever since
omega = 1.12*10^(-5)*86400;
theta_rot_i3 = theta_i3 + omega*day_i3;
theta_rot_i4 = theta_i4 + omega*day_i4;
theta_rot_s1 = theta_s1 + omega*day_s1;
theta_rot_s2 = theta_s2 + omega*day_s2;
theta_rot_s3 = theta_s3 + omega*day_s3;
theta_rot_s4 = theta_s4 + omega*day_s4;
theta_rot_s5 = theta_s5 + omega*day_s5;
theta_rot_s6 = theta_s6 + omega*day_s6;
theta_rot_s7 = theta_s7 + omega*day_s7;
% theta_rot2_i3 = theta_i3 + omega*day_i3 + 0.2506;
% theta_rot2_i4 = theta_i4 + omega*day_i4 + 0.2506;
% theta_rot2_s1 = theta_s1 + omega*day_s1 + 0.2506;
% theta_rot2_s2 = theta_s2 + omega*day_s2 + 0.2506;
% theta_rot2_s3 = theta_s3 + omega*day_s3 + 0.2506;
% theta_rot2_s4 = theta_s4 + omega*day_s4 + 0.2506;
% theta_rot2_s5 = theta_s5 + omega*day_s5 + 0.2506;
% theta_rot2_s6 = theta_s6 + omega*day_s6 + 0.2506;
% theta_rot2_s7 = theta_s7 + omega*day_s7 + 0.2506;

% Select only the good sampling positions
index_i3 = find(radius_i3 < 50);
index_s6 = find(radius_s6 < 50);
index_s7 = find(radius_s7 < 50);

% Select the medusa samples that contained SF6
index_ms3 = index_m3([36:41 48:50]);
index_ms4 = index_m4([4:12]);
index_ms6 = index_m6([21:26]);

h1=polar(theta_rot_s2, radius_s2,'r');
set(h1, 'LineWidth', 2)
hold on

h7=polar(theta_rot_i4, radius_i4);
set(h7, 'LineWidth', 3)

polar(theta_rot_s3(index_ms3), radius_s3(index_ms3),'+k');
h2=polar(theta_rot_s3, radius_s3,'r');
set(h2, 'LineWidth', 2)
polar(theta_rot_s3(index_ms3), radius_s3(index_ms3),'+k');

h3=polar(theta_rot_s4, radius_s4,'r');
polar(theta_rot_s4(index_ms4), radius_s4(index_ms4),'+k');
set(h3, 'LineWidth', 2)

h4=polar(theta_rot_s5, radius_s5,'r');
set(h4, 'LineWidth', 2)

h5=polar(theta_rot_s6(index_s6), radius_s6(index_s6),'r');
set(h5, 'LineWidth', 2)
polar(theta_rot_s6(index_ms6), radius_s6(index_ms6),'+k');

h8=polar(theta_rot_s7(index_s7), radius_s7(index_s7),'r');
set(h8, 'LineWidth', 2)

h6=polar(theta_rot_i3(index_i3), radius_i3(index_i3));
set(h6, 'LineWidth', 3)

hold off

title({'Injection and Sampling Tracks in Eddy Coordinates';'Oceanus 415-2, 23Jul05-1Aug05'})
%; 'Up is NORTH at 0200, 3Aug05 (214.2)'})

text(-27, 20, 'S2')
text(7, 10, 'S6')
text(7, -7, 'S3')
text(-14, -16, 'S5')
text(6, -12, 'S4')
text(-9, 17, 'I3')
text(9, -9, 'S7')
legend('Sampling', 'Injection', 'Medusa')

% subplot(331)
% polar(theta_i3(index_i3), radius_i3(index_i3), '.')
% hold on
% polar(theta_rot_i3(index_i3), radius_i3(index_i3),'.r')
% hold off
%
% subplot(332)
% polar(theta_i4, radius_i4, '.')
% hold on
% polar(theta_rot_i4, radius_i4,'.r')
% hold off
%
% % subplot(334)
% % polar(theta_s1, radius_s1, '.')
% % hold on
% % polar(theta_rot_s1, radius_s1,'.r')
% % hold off
%
% subplot(335)
% polar(theta_s2, radius_s2, '.')
% hold on
% polar(theta_rot_s2, radius_s2,'.r')
% hold off
%
% subplot(336)
% polar(theta_s3, radius_s3, '.')
% hold on
% polar(theta_rot_s3, radius_s3,'.r')
% hold off
%
% subplot(337)
% polar(theta_s4, radius_s4, '.')
% hold on
% polar(theta_rot_s4, radius_s4,'.r')
% hold off
%
% subplot(338)
% polar(theta_s5, radius_s5, '.')
% hold on
% polar(theta_rot_s5, radius_s5,'.r')
% hold off
%
% subplot(339)
% polar(theta_s6(index_s6), radius_s6(index_s6), '.')
% hold on
% polar(theta_rot_s6(index_s6), radius_s6(index_s6),'.r')
% hold off