function [x_df, z_df, x, z, z_mean_high] = dunefoot_plot(tr,year)
%% dunefoot_2nd_deriv
%clear;
%close all;

%% Get JARKUS

url = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/transect.nc';
nargin=1; %change 1

if nargin     == 1
    yrs       = 1:length(nc_varget(url, 'time'));
elseif nargin == 2
    yrs       = find(str2num(datestr(nc_varget(url, 'time')+datenum('01-01-1970','dd-mm-yyyy'),'yyyy'))==year);
    id.time   = ismember(1970 + floor(nc_varget(url, 'time')/365),year);
    id_idtime = find(id.time);
end

% pre-allocate output
x_df = NaN(length(yrs),1); 
z_df = NaN(length(yrs),1);


kustvak         = jarkus_areaCode2AreaName(floor(tr/1e6));
id.alongshore   = nc_varget(url, 'id'  ) == tr;

x               = nc_varget(url, 'cross_shore');
id.cross_shore  = ones(size(x));
id_idalongshore = find(id.alongshore);

z = NaN(length(yrs),1,length(id.cross_shore));

if nargin == 2 %(with years)
    for yy = 1:length(yrs)
        z(yy,1,:) = nc_varget(url, 'altitude', [id_idtime(yy)-1, id_idalongshore-1, 0], [1,1,-1]);
    end
elseif nargin == 1
    z(:,1,:) = nc_varget(url, 'altitude', [0, id_idalongshore-1, 0], [-1,1,-1]);
end

% nodata for that profile, all years
if sum(~isnan(squeeze(z))) == 0
    %1
    return
end

% Water level data
url_WL         = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/MHW_MLW/MHW_MLW.nc';
z_mean_high    = nc_varget(url_WL, 'mean_high_water',id_idalongshore-1,1);

if isempty(z_mean_high)
    disp('No water level data')
    return
end

id_high        = NaN(length(yrs),1);
%% Boundary conditions
profile = z;
profile = squeeze(profile);
height_of_peaks=2.4;

for yy=year:year

    testing_profile=profile(yy,:);
    if nansum(testing_profile(:))==0; %if there are no measured data for the profile
        %3
        continue
    end

    %Remove NaN values
    idnan   =find(isnan(testing_profile), 1);
    x1      =x;
    if ~isempty(idnan)
        x1(isnan(testing_profile))=[];
        testing_profile(isnan( testing_profile))=[];
    end
    
    if isempty(testing_profile) %maybe we can do it with percentage
        %disp('Profile with NaN values')
        continue
    end
    
    [z_prof2, x2, locs] = version1(yy, testing_profile, x1, z_mean_high, height_of_peaks );
    if isempty(z_prof2)
        continue
    end
    
    %Example
    figure()
    subplot(4, 1, [1 2]);

    plot(x, profile(yy,:),'--k')
    hold on
    plot(x1(locs), testing_profile(locs),'ok')
    plot(x2, z_prof2,'k')
    xlim  ([-300 200]);  
    %% First derivative
    
    [z_deriv1, x_deriv1]    = derivative(z_prof2,x2);
    thresholdZderiv1        = 1e-3;                                 % how arbitrary can it be?
    z_deriv1(abs(z_deriv1)  < thresholdZderiv1) = 0;                % "Calibration parameter" 1, if slope smaller than threshold it's zero
    
    subplot(4, 1, 3);
    grid on;
    plot(x_deriv1,1*z_deriv1,'.k') 
    
    % Filter slope
    thresholdFlatDuration = 5; % if there are flat stretches (zderiv1 = 0) longer than 5 units, get rid of the profil more seaward than that
    degree=1;
    [z_deriv1, x_deriv1] = filter_slope(thresholdFlatDuration, z_deriv1, x_deriv1, degree);
    
    %Repeat steps
    [z_deriv1, x_deriv1] = repeat_steps(z_deriv1, x_deriv1, z_prof2, x2, locs, yy, x1, x, profile,testing_profile,z_mean_high, loc_land, height_of_peaks,thresholdFlatDuration);
   
    hold on
    plot(x_deriv1,1*z_deriv1,'*k') 
    xlabel('x [m]','FontSize',14); ylabel('z[m]','FontSize',14)
    xlabel('crossshore [m]')
    ylabel('x [m]','FontSize',14); ylabel('z[m]','FontSize',14)
    ylabel('1st derivative')
    legend('1st derivative', '1st derivative after threshold','Location', 'NorthEast','FontSize',14)
    xlim  ([-300 200]);

    %% Second derivative
    
    [z_deriv2, x_deriv2]    = derivative(z_deriv1,x_deriv1);
    
    subplot(4, 1, 4);
    grid on;
    xlim  ([-300 200]);
    thresholdZderiv2    = 1e-2; 
    x_derivfig          = x_deriv2(z_deriv2<thresholdZderiv2);
    z_derivfig          = z_deriv2(z_deriv2<thresholdZderiv2);
    plot(x_derivfig,1*z_derivfig,'LineStyle','none','Marker','hexagram','Color',[0 0 0])
    

    z_deriv2(z_deriv2<thresholdZderiv2) = 0;                    % "Calibration parameter" 2, very important. 1e-2 proved very good results
    %z_deriv2(abs(z_deriv2)<thresholdZderiv2) = 0;              % should we consider negative derivatives as "breaking points"? before flat areas should we erase the negative values-- or after?
    z_deriv2(isnan(z_deriv2))           = 0;
   % [z_deriv2, x_deriv2] = zero_2nd_derivatives(z_deriv2, x_deriv2);
    
    hold on;
    %plot(x_deriv2,1*z_deriv2,'.r')
    degree=2;
    [z_deriv2, x_deriv2] = filter_slope(thresholdFlatDuration, z_deriv2, x_deriv2, degree);
    z_deriv2_0           = z_deriv2(z_deriv2==0);
    x_deriv2_0           = x_deriv2(z_deriv2==0);
    z_deriv2_no0         = z_deriv2(z_deriv2~=0);
    x_deriv2_no0         = x_deriv2(z_deriv2~=0);
    
    plot(x_deriv2_no0,1*z_deriv2_no0,'LineStyle','none','Marker','square', 'MarkerFaceColor',[0 0 0],'Color',[0 0 0])
    xlabel('x [m]','FontSize',14); ylabel('z[m]','FontSize',14)
    xlabel('crossshore [m]')
    ylabel('x [m]','FontSize',14); ylabel('z[m]','FontSize',14)
    ylabel('2nd derivative')
    legend('2nd derivative', '2nd derivative after threshold','Location', 'NorthEast','FontSize',14)

    xlim  ([-300 200]);
    % if the second derivative is zero but there were values before, then
    % we will decrease the threshold +++++++++++++++++++++++++++++++++
    
    

    dummy                        = z_deriv2;
    dummy(dummy > 0)             = 1;
    dummy(dummy < 0)             = 1;
    id_zero_neighbor             = find(diff(dummy) == -1);
    
    % skip it if it didn't find a second derivative. Handled better with
    % further research
    if isempty(id_zero_neighbor)
        %disp('no neigh')
        %yy
        continue
    end
    
    id_dfp                       = id_zero_neighbor(end);
    x_df(yy)                     = x_deriv2(id_dfp);
    
    z_df(yy)                     = interp1(x(~isnan(z(yy,:,:))),squeeze(z(yy,:,~isnan(z(yy,:,:)))),x_df(yy));

    subplot(4, 1, [1 2]);
    grid on;
    hold on
    plot( x_df(yy), z_df(yy),'ks','MarkerFaceColor','k','MarkerSize',6 )
    xlim  ([-300 200]);
    ylim  ([-3 15]);
    xlabel('x [m]','FontSize',14); ylabel('z[m]','FontSize',14)
    xlabel('crossshore [m]')
    ylabel('z [m]')
    legend('Dune profile', 'Dune peaks','Testing profile', 'Dune foot - New method','Location', 'NorthEast','FontSize',14)
    
end