function [x_df, z_df, x, z, z_mean_high] = Dune_foot_2nd_derivative(tr, year, height_of_peaks, thresholdFlatDuration, thresholdZderiv1, thresholdZderiv2)
%% Dune_foot_2nd_deriv_method

%Calculates the dune foot position based on the geometry of the profiles. The method uses the JARKUS dataset. 

% input:
% tr                      =     id of the transect based on jarkus data
% year                    =     year of interest (1965-2016)
% height_of_peaks         =     the heigh of the peaks that will be considered
%                               as dune peaks(in m)
% thresholdFlatDuration   =     if there are flat stretches longer than this
%                               threshold, the profil more seaward than that is deleted
% thresholdZderiv1        =     "Calibration parameter" 1, for the 1st derivative 
% thresholdZderiv2        =     "Calibration parameter" 2, for the 2nd derivative 
% 
% output:
% x_df                    =     cross-shore position of the dune foot
% z_df                    =     height of the dune foot position

%   --------------------------------------------------------------------

% $Id: $
% $Date: $
% $Author: $
% $Revision: $
% $HeadURL: $
% $Keywords: $
%   --------------------------------------------------------------------

%clear;
%close all;
%% Get JARKUS

%url = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/transect.nc';
url = 'C:\Users\fuentesm\Marilu\Deltares\Projects\Coastviewer\jarkus\scripts\transect_r20240405.nc';
if nargin<=2
    % define default values
    height_of_peaks         =   2.4;
    thresholdFlatDuration   =   5;
    thresholdZderiv1        =   1e-3; 
    thresholdZderiv2        =   1e-2;
end

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

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

if sum(~isnan(squeeze(z))) == 0
    %disp('No data for this profile, for all years')
    return
end

% Water level data
%url_WL         =    'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/MHW_MLW/MHW_MLW.nc';
url_WL         =    'C:\Users\fuentesm\Marilu\Deltares\Projects\Coastviewer\MHW_MLW\processed\MHW_MLW_r20240414.nc';
z_mean_high    =    nc_varget(url_WL, 'mean_high_water',id_idalongshore-1,1);

if isempty(z_mean_high)
    error('No data for water level')
    return
end
%id_high        = NaN(length(yrs),1);

%% Spatial constrains
profile = z;
profile = squeeze(profile);

for yy=1:length(yrs)
    if nargin == 2
        testing_profile = profile;
    else
        testing_profile = profile(yy,:);
    end
    
    if nansum(testing_profile(:)) == 0;   
        %disp('No measured data for the profile')
        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--how?? Is it the same as before?
        %disp('Profile with NaN values')
        continue
    end
    
    [z_prof2, x2, locs, loc_land] = version1(yy, testing_profile, x1, z_mean_high, height_of_peaks );
  
    if isempty(z_prof2)
        continue
    end
      
    %% First derivative
    
    [z_deriv1, x_deriv1]    = derivative(z_prof2, x2);                              
    z_deriv1(abs(z_deriv1)  < thresholdZderiv1) = 0;                % "Calibration parameter" 1: if the slope is smaller than the threshold, it takes a zero value
   
    
    % Filter slope
    % if there are flat stretches (zderiv1 = 0) longer than the thresholdFlatDuration, get rid of the profil more seaward than that
    degree=1;                   % of the derivative
    [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, thresholdZderiv1);
     
    %% Second derivative
    
    [z_deriv2, x_deriv2]    = derivative(z_deriv1,x_deriv1);                           
    z_deriv2(z_deriv2<thresholdZderiv2) = 0;       
    z_deriv2(isnan(z_deriv2))           = 0;

    degree=2;
    [z_deriv2, x_deriv2] = filter_slope(thresholdFlatDuration, z_deriv2, x_deriv2, degree);
    
    dummy                        =  z_deriv2;
    dummy(dummy > 0)             =  1;
    dummy(dummy < 0)             =  1;
    id_zero_neighbor             =	find(diff(dummy) == -1);
    

    if isempty(id_zero_neighbor)
        %disp('no neigh')
        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));  
    
end