%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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);

%
% type = 'asl';
% elevation = [10,50,100,150,500,600,700,800,900,1000,1500,2000,2500];
type = 'agl';
% elevation = [600,700,800,900];
elevation = [10,50,100,150,500,600,700,800,900,1000,1500,2000,2500];

nelev = length(elevation);
% load RUWRF Files
for iz=1:nelev
    model=strcat(num2str(elevation(iz)),'m',type);
    switch model
        case '10masl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_10m_asl(fnames);
        case '50masl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_50m_asl(fnames);
        case '100masl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_100m_asl(fnames);
        case '150masl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_150m_asl(fnames);
        case '500masl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_500m_asl(fnames);
        case '1000masl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_1000m_asl(fnames);
        case '1500masl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_1500m_asl(fnames);
        case '2000masl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_2000m_asl(fnames);
        case '2500masl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_2500m_asl(fnames);
        case '10magl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_10m_agl(fnames);
        case '50magl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_50m_agl(fnames);
        case '100magl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_100m_agl(fnames);
        case '150magl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_150m_agl(fnames);
        case '500magl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_500m_agl(fnames);
        case '600magl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_600m_agl(fnames);
        case '700magl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_700m_agl(fnames);
        case '800magl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_800m_agl(fnames);
        case '900magl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_900m_agl(fnames);
        case '1000magl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_1000m_agl(fnames);
        case '1500magl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_1500m_agl(fnames);
        case '2000magl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_2000m_agl(fnames);
        case '2500magl'
            [Lon,Lat,U,V,TT] = loadRUWRFFile_2500m_agl(fnames);
    end
    
    
    %%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 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'),...
            '_',model,conf.Traj.FileSuffix));
        save(filename,'TRAJ');
    catch
        fprintf('Saving trajectory failed because:\n')
        res = lasterror;
        fprintf('%s\n',res.message)
    end
    
    %% LCS stuff
    try
        
        % This section calculates LCS
        % default variables
        startFrame = 1;
        finalFrame = size(TRAJ.TimeStamp,2);
        o = [1,0]; % direction [1,0]
        
        % convert from degree to kilometer
        % define an arbitrary reference origin
        lon_orig=TRAJ.InitialLonLat(1,1);
        lat_orig=TRAJ.InitialLonLat(1,2);
        
        [east, north] = lonlat2km(lon_orig,lat_orig,TRAJ.Lon,TRAJ.Lat);
        
        [...
            AD,...
            D,...
            MAD,...
            XMaxD,...
            XMinD,...
            XMAD,...
            XMean,...
            YMaxD,...
            YMinD,...
            YMAD,...
            YMean,...
            XYMaxD,...
            XYMinD,...
            XYMAD,...
            XYMean] = compute_met(...
            east,...
            north,...
            TRAJ.TimeStamp,...
            'startFrame',startFrame,...
            'finalFrame',finalFrame,...
            'direction',o);
        
        % set the last frame
        EAST = east(:,end);
        NORTH = north(:,end);
        
        % Relative Dispersion
        % generates a rectangular grid
        [ln,lt] = meshgrid_vector_data(TRAJ.InitialLonLat(:,1),TRAJ.InitialLonLat(:,2));
        [d1,d2] = size(ln);
        
        EAST = reshape(EAST,d1,d2);
        NORTH = reshape(NORTH,d1,d2);
        % convert from degree to kilometer
        RD = compute_rd(EAST,NORTH);
        
        InitialLonLat=TRAJ.InitialLonLat;
        TimeStamp=TRAJ.TimeStamp(end);
        ad(:,iz)=AD(:,end);
        mad(:,iz)=MAD(:,end);
        xmad(:,iz)=XMAD(:,end);
        ymad(:,iz)=YMAD(:,end);
        rd(:,iz)=RD(:);
        
        
    catch
        fprintf(' LCS calculation failed because:\n')
        res = lasterror;
        fprintf('%s\n',res.message)
    end
end

% Saving LCS
try
    if ~exist( conf.LCS.BaseDir, 'dir' )
        fprintf('%s is creating the directory:\n%s\n',mfilename,conf.LCS.BaseDir);
        mkdir(conf.LCS.BaseDir);
    end
    
    filename = fullfile(...
        conf.LCS.BaseDir,...
        strcat(conf.LCS.FilePrefix,...
        datestr(startdate,'yyyy_mm_dd_HHMM'),...
        '_',type,conf.LCS.FileSuffix));
    save(filename,'ad','mad','xmad','ymad','rd','InitialLonLat','TimeStamp');
catch
    fprintf('Saving LCS failed because:\n')
    res = lasterror;
    fprintf('%s\n',res.message)
end

% Plotting
try
    
    figure
    for iz=1:nelev
        lcs = reshape(rd(:,iz),d1,d2);
        surf(ln,lt,0.*lcs+elevation(iz),lcs)
        hold on
    end
    zlim([1 elevation(nelev)])
    grid on
    shading interp
    zlabel('depth m')
    ylabel('latitude (^o)');
    xlabel('longitude (^o)');
%     set(gca,'zdir','reverse')
    view(-20,10)
    set(gca,'ztick',0:500:elevation(nelev))
    
catch
    fprintf('plot LCS failed because:\n')
    res = lasterror;
    fprintf('%s\n',res.message)
end