function [varargout] = hovmoller(Xdata,Ydata,Zdata,Tdata,Data,varargin)
%==========================================================================
%
% DESCRIPTION:
% A Hovmoller diagram is a commonly used way of plotting meteorological
% data to highlight the role of waves. The axes of a Hovmoller diagram are
% typically longitude or latitude (abscissa or x-axis) and time (ordinate
% or y-axis) with the value of some field represented through color or
% shading.
%
% INPUT:
% Tdata: time
% Xdata: Longitude
% Ydata: Latitude
% Zdata: Elevation
% Data: Data(X,Y,Z,T)
% Options
%
% OUTPUT:
% transect
%
% Usage:
% [transect] =
% Author(s)
%   Fredj Erick (The Jerusalem College of Technology)
%
% Revision history:
%   2020-03-30 First version.
%
%==========================================================================

global ftbverbose

% global ftbverbose
subname = 'hovmoller';
if ftbverbose
    fprintf('\nbegin : %s\n', subname)
end

% default parameter
lonlims = [-180 180];
latlims  = [-90 90];
statebdryfile   = '';
boemfile        = '';
transectfile     = fullfile(pwd,'db_transect.csv');
coastfile = '/Users/seroka/Research/code/geo/eastcoastmask_new.mat';
longitude_coastal_roi_vertices = []; % longitude-coordinates of polygon vertices
latitude_coastal_roi_vertices = []; % latitude-coordinates of polygon vertices
method = 'L1';
k=1;
while k<length(varargin)
    switch lower(varargin{k})
        case 'coastfile'
            coastfile=varargin{k+1};
            varargin([k k+1])=[];
        case 'lonlims'
            lonlims = varargin{k+1};
            varargin([k k+1])=[];
        case 'latlims'
            latlims = varargin{k+1};
            varargin([k k+1])=[];
        case 'statebdryfile'
            statebdryfile = varargin{k+1};
            varargin([k k+1])=[];
        case 'boemfile'
            boemfile = varargin{k+1};
            varargin([k k+1])=[];
        case 'transectfile'
            transectfile = varargin{k+1};
            varargin([k k+1])=[];
        case 'coastfile'
            coastfile = varargin{k+1};
            varargin([k k+1])=[];
        case 'longitude_coastal_roi_vertices'
            longitude_coastal_roi_vertices = varargin{k+1};
            varargin([k k+1])=[];
        case 'latitude_coastal_roi_vertices'
            latitude_coastal_roi_vertices = varargin{k+1};
            varargin([k k+1])=[];
        case 'method'
            method = varargin{k+1};
            varargin([k k+1])=[];
        otherwise
            k=k+2;
    end
end

% Load transect dataset
if ~exist(transectfile,'file')
    error('transect cvs file into %s is missing', subname)
else
    warning('default transect cvs is %s', transectfile);
end

% Load coastline file
coast = load(coastfile);

%Load federal state boundary
if exist(statebdryfile,'file')
    load(statebdryfile);
end

%Load BOEM call area
if exist(boemfile,'file')
    load(boemfile)
    boem_lon = X2013_boem_osw_call_area(:,2);
    boem_lat = X2013_boem_osw_call_area(:,4);
end

% extract elevation data
[etopo2_elev,etopo2_long,etopo2_lat]=m_etopo2([lonlims,latlims]);

% Data dimension
[nx,ny,nz,nt] = size(Data);
[X,Y] = meshgrid(Xdata,Ydata);

% Load transect dataset
dbtransect = importTransectFile(transectfile, [2, Inf]);
number_of_transect = size(dbtransect,1);

for i = 1:number_of_transect
    
    transect_reference_latitude = dbtransect.LATITUDE(i);
    transect_reference_longitude = dbtransect.LONGITUDE(i);
    transect_location = dbtransect.LOCATION(i);
    transect_distance = dbtransect.DISTANCE(i);
    transect_bearing = dbtransect.BEARING(i);
    transect_spacing = dbtransect.SPACING(i);
    
    nwaypoints = round(transect_distance/transect_spacing);
    
    transect_wp = [transect_reference_latitude,transect_reference_longitude]; % starting point
    for n = 2: nwaypoints
        transect_wp(n,:)=course2ll([transect_wp(n-1,1) transect_wp(n-1,2) transect_bearing transect_spacing]);
    end
    
    % extract the elevation, and the Data extract transect.
    transect_wp_elev = NaN(nwaypoints,1);
    T3d = NaN(nwaypoints,nz,nt);
    for n = 1: nwaypoints
        transect_wp_lat = transect_wp(n,1);
        transect_wp_long = transect_wp(n,2);
        % Elevation
        [Point,~] = find_nearest_point(transect_wp_long,transect_wp_lat,etopo2_long(:),etopo2_lat(:));
        % Subscripts from linear index
        [ilat,ilon] = ind2sub(size(etopo2_elev),Point);
        transect_wp_elev(n,1) = etopo2_elev(ilat,ilon);
        % Data
        [Point,~] = find_nearest_point(transect_wp_long,transect_wp_lat,X(:),Y(:));
        [ilat,ilon] = ind2sub(size(X),Point);
        T3d(n,:,:) = Data(ilon,ilat,:,:);
    end
    
    % the intersection points of two curves transect, and coastline.
    %P = InterX([transect_wp(:,2)';transect_wp(:,1)'],[coast.long';coast.lat']);
        P = InterX([transect_wp(:,2)';transect_wp(:,1)'],[coast.lon';coast.lat']);

    if ~isempty(longitude_coastal_roi_vertices)||~isempty(latitude_coastal_roi_vertices)
        [in,~]=inpolygon(P(1,:),P(2,:),longitude_coastal_roi_vertices,latitude_coastal_roi_vertices);
        coastline_reference_point_longitude = P(1,in);
        coastline_reference_point_latitude = P(2,in);
    else
        coastline_reference_point_longitude = P(1,end);
        coastline_reference_point_latitude = P(2,end);
    end
    
    m_proj('mercator','lat',latlims,'lon',lonlims);
    transect_wpkm = NaN(nwaypoints,1);
    transect_wp_longitude=NaN(nwaypoints,1);
    transect_wp_latitude=NaN(nwaypoints,1);
    
    for p=1: nwaypoints
        [range,ln,lt]=m_lldist([coastline_reference_point_longitude transect_wp(p,2)],[coastline_reference_point_latitude transect_wp(p,1)],40);
        transect_wpkm(p) = range(end);
        transect_wp_longitude(p) = ln(end);
        transect_wp_latitude(p) = lt(end);
    end
    s=sign(diff(transect_wpkm));
    final_transect = s.*transect_wpkm(1:length(s));
    npt = length(s);
    
    transect(i).offshore_distance = final_transect;
    transect(i).longitude=transect_wp_longitude(1:npt);
    transect(i).latitude =  transect_wp_latitude(1:npt);
    transect(i).bathy = transect_wp_elev(1:npt);
    transect(i).coast_reference_longitude = coastline_reference_point_longitude;
    transect(i).coast_reference_latitude= coastline_reference_point_latitude;
    transect(i).name = transect_location;
    transect(i).bearing = transect_bearing;
    transect(i).data = T3d(1:npt,:,:);
    transect(i).elevation = Zdata;
    transect(i).time = Tdata;
    
    [hovmoller_data] = worker_hovmoller_func(transect(i).data,transect(i).elevation,method);
    transect(i).hovmoller_data = hovmoller_data;
    transect(i).hovmoller_method = method;
    
    %%%%%%%%%%%     Plot    %%%%%%%%%%%%%%%%
    figure(1);
    set(gcf,'Paperposition',[0 0 9 9]);
    
    % Plot Basemap
    subplot(4,1,1:3)
    m_proj('mercator','lons',lonlims,'lat',latlims);
    [cs,h]=m_etopo2('contourf',[-5000:500:0],'edgecolor','none');
    %m_gshhs_h('patch',[.5 .8 0],'edgecolor','none');
    m_gshhs_l('patch',[.5 .8 0],'edgecolor','none');
    m_grid('linewi',2,'layer','top');
    m_contfbar(.92,[.2 .5],cs,h,'endpiece','no','axfrac',.02);
    colormap(m_colmap('blue'));
    m_northarrow(-78,41,1,'type',4,'linewi',.5);
    m_plot(transect(i).longitude,transect(i).latitude,'marker','o','color','y','markersize',8,'markerfacecolor','y')
    m_plot(transect(i).coast_reference_longitude,transect(i).coast_reference_latitude,'+')
    m_text(transect(i).coast_reference_longitude,transect(i).coast_reference_latitude,transect(i).name)
    xlabel('Longitude ^o')
    ylabel('Latitude ^o')
    
    %Plot bathymetry along cross section
    subplot(4,1,4)
    plot(transect(i).offshore_distance,transect(i).bathy,'.-b');
    ylabel('Terrain/Water Depth (m)');
    xlabel(sprintf('Distance Offshore in kilometer from lat: %3.2f ^o, lon: %3.2f ^o ',transect(i).coast_reference_latitude,...
        transect(i).coast_reference_longitude));
    
    
    figure(2); clf
    [TT,OD] = meshgrid(transect(i).time,transect(i).offshore_distance);
    contourf(OD,TT,transect(i).hovmoller_data-10)
    xlabel(sprintf('Distance Offshore in kilometer from lat: %3.2f ^o, lon: %3.2f ^o ',transect(i).coast_reference_latitude,...
        transect(i).coast_reference_longitude));
    ylabel('Time');
    datetick('y');
    caxis([-6 6]);
    colormap redblue
    h = colorbar;
    %set(get(h,'title'),'string','hovmoller');
    set(get(h,'title'),'string','RD-10km (km/hr)');
    set(gca, 'Ydir','reverse')
    grid on
end

varargout{1} = transect;
if ftbverbose
    fprintf(['end   : ' subname '\n'])
end

end

%==========================================================================
function [hovmoller_data] = worker_hovmoller_func(data,z,method,varargin)
% This function will be used to actually evaluate the hovmoller diagram.
%
% INPUT:
% data: 
% z: elevation
% atmospheric Layer have 3 options 
% L1: layer range in meter [50 150]
% L2: layer range in meter [0 500]
% L3: layer range in meter [0 2000]
% 
% OUTPUT:
% hovmoller_data
% 
% Author(s)
%   Fredj Erick (The Jerusalem College of Technology)
%
% Revision history:
%   2020-04-01 First version.
%
%==========================================================================

% size
[nx,nz,nt] = size(data);

switch (upper(method))
    case 'L1'
        subElevation = [50 150]; % L1 - Elevation range in meter
    case 'L2'
        subElevation = [0 500];
    case 'L3'
        subElevation = [0 2000];
    otherwise
        error('elevation range %s is unknown', method);
end
zidx=z>=subElevation(1) & z<=subElevation(2);

% memory allocation
mean_hovmoller_data = NaN(nx,nt);
for i=1:nx
    for t=1:nt
        mean_hovmoller_data(i,t) = mean(data(i,zidx,t));
    end
end

hovmoller_data = mean_hovmoller_data;

end