% fl_pos_fit_blair  float position fit for Sargasso Eddy
% =========================================================================
% fl_pos_fit_blair  Version 1.0 02-Sep-2005
%
% Usage: 
%   fl_pos_fit_blair
%
% Description:
%   This function sets the initial conditions for fitting the drogue drifter
%   data for the Sargasso Eddy A4.  This eddy is propagting to the southwest
%   at about 4-5 km/day and has a rotation period of approximately 6-7 days.
%   
% Input:
%
% Output:
%   Plot of original data and results of fit
%
% Called by:
%   n/a
%
% Calls:
%   n/a
%
% Author:
%   Cindy Sellers WHOI
%   Blair Greenan Bedford Institute of Oceanography
%   James Girton APL, University of Washington
%   Matlab 7
%   greenanb@mar.dfo-mpo.gc.ca
% =========================================================================
%

% History
% Version 1.0 02-Sep-2005

drogue2use = '52484';
drogue2use = '52486';
drogue2use = '52485';
drogues = {'52484' '52486' '52485'};

now = datenum(2005,09,08,00,00,00);

figure(1);
%clf;

for ii = 1:length(drogues);
  drogue2use = drogues{ii};

% Set an intial latitude and longitude in the vicinity of the eddy for reference
lonref=-69;
latref=30;

switch lower(drogue2use)
  case '52484'
    % Initial latitude and longitude of the eddy center

    lon0 = -69.220;
    lat0 = 29.733;

    % guess radius of drogue drifter path in meters (m)
    %r=8.4e3; 
    r=1.5e4; 
    
    % guess initial angle (rad ccw from E)
    %th0 = 0;
    th0 = -0.8;
        
    % translation velocities (m/s)
    u = -0.05; % positive to the east
    v = 0.01; % positive to the north
    % rotation rate (rad/s)
    omega = -1.1e-5;
  case '52486'
    lon0 = -69.180;
    lat0 = 29.783;
    r=1.5e4; 
    th0 = -0.5;
    u = -0.05; % positive to the east
    v = 0.00; % positive to the north
    omega = -1.1e-5; % 7.16 days to complete one rotation
  case '52485'
    % initial location of center
    lat0 = 29.76;
    lon0 = -69.40;

    % guess radius of drogue drifter path in meters (m)
    %r=8.4e3; 
    r=2e4; 

    % guess initial angle (rad)
    th0 = -0.3;
    
    % guess angular velocity of eddy (rad/s) - 1.4e-5 = 5.25 day period of rotation
    %omega = -1.4e-5;
    omega = -1.1e-5;
    
    % guess initial velocity of eddy center (m/s)
    u = -0.07;
    v = 0.02; 

end % switch lower(drogue2use)

phi = lat0*pi/180;
x00 = 111000/cos(phi)*(lon0 - lonref); % convert from degrees to distance in m
y00 = 111000*(lat0-latref);


% create a vector to pass to the fitting routine
Z0=[x00 y00 r th0 omega u v];

% set options for the fminsearch function and call this function
options = optimset('MaxFunEvals',1e20,'MaxIter',1e20);
[X,fval,exitflag,output]=fminsearch(@my_fun_blair,Z0,options,drogue2use);

% Display results of the fit on the command line
xout = X(1);
yout = X(2);
rout = X(3);
th0_out = X(4);
omega_out = X(5);
uout = X(6);
vout = X(7);

lon_out = lonref + xout/111000*cos(phi);
lat_out = yout/111000 + latref;


%%%%%%%%%%%%%% PLOT RESULTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

phi = lat0*pi/180; % rad
x=111000/cos(phi)*(lon - lonref); % m
y=111000*(lat-latref); % m
% make a single value
z=x+i*y;  

% reference time
t0 = t(1)-.001; % days
% convert to seconds
tsec_ref = (t-t0)*3600*24;

figure(ii+1)
clf
zfun = ((xout+i*yout) + rout * exp(i*(th0_out+omega_out*tsec_ref)) + (uout+i*vout)*tsec_ref ); 
plot([x real(zfun)]',[y imag(zfun)]','y.-'); hold on;
plot(real(zfun),imag(zfun),'.-',x,y,'.-'); grid on
set(gca,'dataaspectratio',[1 cos(lat0*pi/180) 1]);
hold off



nlon = (xout+uout*(now-t0-datenum(2005,1,1,0,0,0))*86400)/111000*cos(phi) + lonref;
nlat = (yout+vout*(now-t0-datenum(2005,1,1,0,0,0))*86400)/111000 + latref;
disp(['Drogue ' drogue2use]);
disp(['Center now (' datestr(now) '): ' num2str(nlat) ' N, ' num2str(-nlon) ' W']);
disp(['Rotation period: ' num2str(2*pi/(-omega_out*86400)) ' days']);
disp(['Eddy speed: U=' num2str(uout,2) ' m/s, V=' num2str(vout,2) ' m/s']);

figure(1)
ph = plot(nlon,nlat,'p');
grid on
cord = get(gca,'colororder');
col = cord(ii,:);
set(ph,'color',col);
hold on;


end % for drogue2use = drogues,
