%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function performs the particle tracking from Antartica HFR current data.
%
% References:
% http://www.mathworks.com/matlabcentral/answers/93554-how-can-i-rotate-a-set-of-points-in-a-plane-by-a-certain-angle-about-an-arbitrary-point
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 	$Id: PLDP_HFRPTrajectory.m 001 2015-06-09 10:55:10Z ef $
%
% Copyright (C) 2014-2015 Erick Fredj
% Licence: GPL (Gnu Public License)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
clear all
clc

thisDir = fileparts( mfilename('fullpath') );

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Basic info on the Total files.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
conf = MABSB_conf;

% time period
startdate=datenum(2013,04,27,18,0,0);
finaldate = datenum(2013,04,27,23,0,0);
dtime = startdate:1/24:finaldate;

% Generate all filenames at once using repmat appropriately
fnames=datenum_to_directory_filename(conf.WRFOUT.BaseDir,dtime,...
                                     conf.WRFOUT.FilePrefix,conf.WRFOUT.FileSuffix,conf.WRFOUT.MonthFlag);

% load RUWRF Files
[Lon,Lat,U,V,TT] = loadRUWRFFile_500(fnames);

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do some particle tracking on the grid.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp( 'Doing some particle tracking on a rectangular Lon,Lat grid.' );

% Domain: MAB: d01: SEA_BREEZE (MABSB)
% generate the drifter release points
lonlims = [-79.52 -68.48]; % (degree)
latlims = [35.37 43.52]; % (degree)
heightlims=[11 3000]; % (meter)

% % Domain: MAB: d02: SEA_BREEZE (MABSB)
% % generate the drifter release points
% lonlims = [-75.98 -72.95]; % (degree)
% latlims = [38.53 40.97]; % (degree)
% heightlims=[11 3000]; % (meter)

dlonlat=[10.0,10.0];
% dlonlat=[1.0, 1.0];
% dlonlat=[0.5,0.5];
% dlonlat=[0.10, 0.10];
du='km';
%
[lon_MABSB_tracer,lat_MABSB_tracer] = LonLat_grid([min(lonlims),min(latlims)],[max(lonlims),max(latlims)],dlonlat,du );
lon_MABSB_tracer=lon_MABSB_tracer(:);
lat_MABSB_tracer=lat_MABSB_tracer(:);
%
ln = lon_MABSB_tracer;
lt = lat_MABSB_tracer;
LL = [ ln(:), lt(:) ];


% Get some good option values for tracking - otherwise output is junk.
abs_tol = 1.0e-3; % Not sure about this
rel_tol = 1.0e-3; % Not sure about this
maxstep = 1/24; % 1 hour  
options = odeset('RelTol',rel_tol,'AbsTol',abs_tol,'MaxStep',maxstep);

tspan = startdate:1/24/6:finaldate; % every hour

% This calculates trajectories and puts in a TRAJ structure with some
% extra metadata (but not much).
TRAJ = ptrack2TRAJstruct('particle_track_ode_grid',Lon,Lat,...
                    U,V,TT,tspan,LL,options);
                
TRAJ.TrajectoryDomain = conf.WRFOUT.DomainName;
TRAJ.OtherMetadata.ptrack2TRAJstruct.options = options;

% Useful for plotting
LonLims = [ min(Lon(1,:)), max(Lon(1,:)) ];
LatLims = [ min(Lat(:,1)), max(Lat(:,1)) ];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set up basemap
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clf
coastFileName='coastMABSB';

plotBasemap( LonLims, ...
			LatLims , ...
			coastFileName,'mercator','patch','g');
hold on

th = title( 'Select location of distance bar', 'fontsize', 16 );

[dh,dth,dl] = m_distance_bar( 10,[-79,42] );
set(dh,'linewidth',2)
set(dth,'fontsize',14)

ppp = get(dth,'position');
ppp(3) = 1000;
set(dth,'position',ppp)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot particle tracks
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ph2 = m_plot( TRAJ.Lon', TRAJ.Lat', 'linewidth', 4, 'color', 'k' );

th = title( { 'Total Currents (color), Particle tracks (black) ', ...
	'Click to see Particle Tracks' }, 'fontsize', 16 );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Saving trajectory
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
try
    if ~exist( conf.Traj.BaseDir, 'dir' )
        fprintf('%s is creating the directory:\n%s\n',mfilename,conf.Traj.BaseDir);
        mkdir(conf.Traj.BaseDir);
    end
    
    filename = fullfile(...
        conf.Traj.BaseDir,...
        strcat(conf.Traj.FilePrefix,...
        datestr(startdate,'yyyy_mm_dd_HHMM'),...
        '_500m_',conf.Traj.FileSuffix));
    save(filename,'TRAJ');
catch
    fprintf('Saving trajectory failed because:\n')
    res = lasterror;
    fprintf('%s\n',res.message)
end

