% PLOTFITS  Plot results of eddy fits to drifter tracks (from eddyfit)

do_plotpr = 0;
do_plotuv = 0;
do_plotpos = 0;
do_plotfits = 1;
do_addsubj = 1;

ef1 = load('edfit_52484.mat');
ef2 = load('edfit_52486.mat');
ef3 = load('edfit_52485.mat');

  mld1 = datenum(2005,8,28,0,0,0);
%  mld1 = datenum(2005,8,8,0,0,0);
%  mld2 = datenum(2005,9,8,23,0,0);
  mld2 = now+3/24;

if do_plotpr
subplot(211)
plot(ef1.mld_mid,-2*pi./ef1.omega_fit/86400,'.-');
grid on
hold on
plot(ef2.mld_mid,-2*pi./ef2.omega_fit/86400,'r.-');
plot(ef3.mld_mid,-2*pi./ef3.omega_fit/86400,'k.-');
legend('52484','52486','52485',2);
ylabel('Rotation period (days)');
axis([mld1 mld2 4 8.5]);
datetick('x','keepticks','keeplimits')

subplot(212)
plot(ef1.mld_mid,ef1.rad_fit*1e-3,'.-');
grid on;
hold on
plot(ef2.mld_mid,ef2.rad_fit*1e-3,'r.-');
plot(ef3.mld_mid,ef3.rad_fit*1e-3,'k.-');
ylabel('Drifter Radius (km)');
xlabel('Date in 2005)');
  axis([mld1 mld2 5 35]);
datetick('x','keepticks','keeplimits')
end % if do_plotpr

if do_plotuv
subplot(211)

[dd1,th1] = sw_dist(ef1.lat_fit,ef1.lon_fit,'km');
u1 = dd1.*cos(th1*pi/180)./diff(ef1.mld_mid)/86400*1e3;
v1 = dd1.*sin(th1*pi/180)./diff(ef1.mld_mid)/86400*1e3;
t1 = (ef1.mld_mid(1:end-1)+ef1.mld_mid(2:end))/2;
[dd2,th2] = sw_dist(ef2.lat_fit,ef2.lon_fit,'km');
u2 = dd2.*cos(th2*pi/180)./diff(ef2.mld_mid)/86400*1e3;
v2 = dd2.*sin(th2*pi/180)./diff(ef2.mld_mid)/86400*1e3;
t2 = (ef2.mld_mid(1:end-1)+ef2.mld_mid(2:end))/2;
[dd3,th3] = sw_dist(ef3.lat_fit,ef3.lon_fit,'km');
u3 = dd3.*cos(th3*pi/180)./diff(ef3.mld_mid)/86400*1e3;
v3 = dd3.*sin(th3*pi/180)./diff(ef3.mld_mid)/86400*1e3;
t3 = (ef3.mld_mid(1:end-1)+ef3.mld_mid(2:end))/2;

plot(ef1.mld_mid,ef1.u_fit,'b.-')
grid on;
hold on
plot(ef2.mld_mid,ef2.u_fit,'r.-');
plot(ef3.mld_mid,ef3.u_fit,'k.-');
legend('52484','52486','52485',2);
ylabel('E-W Eddy Velocity (m/s)');

plot(t1,u1,'b^--');
plot(t2,u2,'r^--');
plot(t3,u3,'k^--');

title('Dots: fit velocity; Triangles: fit position differences');
  axis([mld1 mld2 -0.15 0.05]);
datetick('x','keepticks','keeplimits')


subplot(212)
plot(ef1.mld_mid,ef1.v_fit,'b.-');
grid on;
hold on
plot(ef2.mld_mid,ef2.v_fit,'r.-');
plot(ef3.mld_mid,ef3.v_fit,'k.-');
ylabel('N-S Eddy Velocity (m/s)');
xlabel('Date in 2005)');

plot(t1,v1,'b^--');
plot(t2,v2,'r^--');
plot(t3,v3,'k^--');

axis([mld1 mld2 -0.15 0.05]);
datetick('x','keepticks','keeplimits')

end % if do_plotuv

if do_plotpos

  subplot(221)

  if do_addsubj
    % Add subjective eddy center from nightly ADCP examination
  [dstr1,dstr2,clat,clon] = textread('../adcp/positions/trajectory.dat',...
      '%s     %s    %7.4f %8.4f','headerlines',1);
  clear cmld
  for ind = 1:length(dstr1),
    cmld(ind) = datenum([dstr1{ind} ' ' dstr2{ind}(1:2)...
	  ':' dstr2{ind}(3:4) ':00'],0);
  end
  plot(cmld,clat,'ko','markerfacecolor','y'); hold on
  end % if do_addsubj

  plot(ef1.mld_mid,ef1.lat_fit,'b.-'); grid on
  hold on
  plot(ef2.mld_mid,ef2.lat_fit,'r.-');
  plot(ef3.mld_mid,ef3.lat_fit,'k.-');
  axis([mld1 mld2 29.6 30]);

  datetick('x','keepticks','keeplimits')
  if do_addsubj
    legend('Manual Picks','52484','52486','52485',1);
  else
    legend('52484','52486','52485',1);
  end
  ylabel('lat');

  
  subplot(223)
  if do_addsubj
    plot(cmld,clon,'ko','markerfacecolor','y'); hold on
  end
  plot(ef1.mld_mid,ef1.lon_fit,'b.-'); grid on
  hold on
  plot(ef2.mld_mid,ef2.lon_fit,'r.-');
  plot(ef3.mld_mid,ef3.lon_fit,'k.-');
  xlabel('UTC')
  axis([mld1 mld2 -69.6 -69]);
  datetick('x','keepticks','keeplimits')
  ylabel('lon');

  
  subplot(122)
  if do_addsubj
    plot(cmld,clon,'ko','markerfacecolor','y'); hold on
  end
  plot(-ef1.lon_fit,ef1.lat_fit,'b.'); grid on
  hold on; set(gca,'xdir','reverse');
  plot(-ef2.lon_fit,ef2.lat_fit,'r.');
  plot(-ef3.lon_fit,ef3.lat_fit,'k.');
  mlat = mean(ef1.lat_fit);
  axis([69 69.6 29.6 30]);
  set(gca,'dataaspectratio',[1 cos(mlat*pi/180) 1]);
  xlabel('longitude (\circ W)'); ylabel('latitude (\circ N)');
  title('Eddy Center');

  if do_addsubj
  plot(-clon,clat,'ko','markerfacecolor','y');
   end 
 end

if do_plotfits
  % Go through and plot 7-day model trajectory for each fit

  drogues = {'52484' '52486' '52485'};
  %drogues = {'52486' '52484' '52485'};
  filenames = {'l_argos_pos.txt' 'drft_shp.txt' 'drift_52485.txt'};
  timewin = 7; % number of days in window
  cols = {{'b' 'y'} {'r' 'k'} {'c' 'm'}};
  
  if 0
  plot(-ef1.lon_fit,ef1.lat_fit,'bp'); grid on
  hold on; set(gca,'xdir','reverse');
  plot(-ef2.lon_fit,ef2.lat_fit,'rp');
  plot(-ef3.lon_fit,ef3.lat_fit,'cp');
  end

  mlat = mean(ef1.lat_fit);
  %axis([69 69.6 29.6 30]);
  %set(gca,'dataaspectratio',[1 cos(mlat*pi/180) 1]);
  %title('Eddy Center');

  for ii = 1:length(drogues);
    drogue2use = drogues{ii}
    eval(['ef = ef' num2str(ii) ';']);

    dat = load(['Drogue_' drogue2use '_smoothed.txt']);
    t = dat(:,1);
    lon=dat(:,2);
    lat= dat(:,3);

    subplot(3,1,ii);
    plot(-lon,lat,'k.-');
    %plot(-lon,lat,'.-','color',col{1});
    xlabel('longitude (\circ W)'); ylabel('latitude (\circ N)');
    hold on; set(gca,'xdir','reverse');
    

    col = cols{ii};
    ic = 1; nc = length(col);
    for jj = 1:length(ef.mld_mid),

      
      % Convert fit parameters into modeled drifter position (km)
      tsec_ref = [-timewin/2:1/24:timewin/2]*86400;
      phi = ef.lat_fit(jj)*pi/180; % rad

      zfun = ef.rad_fit(jj) * exp(i*(ef.theta_fit(jj)+ef.omega_fit(jj)*tsec_ref)) + ... 
	  (ef.u_fit(jj)+i*ef.v_fit(jj))*tsec_ref ; 

      % Convert position back into lat,lon
      dlon = real(zfun)/1.11e5*cos(phi) + ef.lon_fit(jj);
      dlat = imag(zfun)/1.11e5 + ef.lat_fit(jj);
      plot(-dlon,dlat,'-','linewidth',0.5,'color',col{ic});

      
      ic = ic+1;
      if ic>nc,
	ic = ic - nc;
      end
      %pause
    end %     for jj = 1:length(ef.mld_mid),
  set(gca,'dataaspectratio',[1 cos(mlat*pi/180) 1]);
    plot(-lon,lat,'k.-');
    plot(-ef.lon_fit,ef.lat_fit,'kd','markerfacecolor',col{1},...
	'markersize',7); grid on
    %plot(-ef.lon_fit,ef.lat_fit,'kp'); grid on
  
  if do_addsubj
    % Add subjective eddy center from nightly ADCP examination
    [dstr1,dstr2,clat,clon] = textread('../adcp/positions/trajectory.dat',...
	'%s     %s    %7.4f %8.4f','headerlines',1);
    plot(-clon,clat,'ko','markerfacecolor','y');
   end 
  end %  for ii = 1:length(drogues);
  
end % if do_plotfits
