% routine to get MHW and MLW
%

%% settings
% jarkus database - dataset 1965 - 2010
OPT.ncpath2      = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/jarkus/profiles/';
OPT.ncpath2      = 'F:\JARKUSdata\';
OPT.Filename2    = 'transect.nc';

%% jarkus database - dataset 1965 - 2010
%% Read data
id       = nc_varget([OPT.ncpath2 OPT.Filename2],'id');
MHWm     = nc_varget([OPT.ncpath2 OPT.Filename2],'mean_high_water');
MLWm     = nc_varget([OPT.ncpath2 OPT.Filename2],'mean_low_water');
lat      = nc_varget([OPT.ncpath2 OPT.Filename2],'lat');
lon      = nc_varget([OPT.ncpath2 OPT.Filename2],'lon');
cross    = nc_varget([OPT.ncpath2 OPT.Filename2],'cross_shore');
along    = nc_varget([OPT.ncpath2 OPT.Filename2],'alongshore');
areaname = nc_varget([OPT.ncpath2 OPT.Filename2],'areaname');
areacode = nc_varget([OPT.ncpath2 OPT.Filename2],'areacode');
time     = nc_varget([OPT.ncpath2 OPT.Filename2],'time');

rsp_x    = nc_varget([OPT.ncpath2 OPT.Filename2],'rsp_x');
rsp_y    = nc_varget([OPT.ncpath2 OPT.Filename2],'rsp_y');
x        = nc_varget([OPT.ncpath2 OPT.Filename2],'x');
y        = nc_varget([OPT.ncpath2 OPT.Filename2],'y');
angle    = nc_varget([OPT.ncpath2 OPT.Filename2],'angle');

%% Evaluation of MHW

% counters
[kv, m, ~] = unique(areacode);
n          = vertcat(1, m(1:end-1)+1);

% definition
MHW_lat  = NaN(length(areacode),length(time)); MHW_lon  = NaN(length(areacode),length(time));
MHW_rsp  = NaN(length(areacode),length(time));

%% MHW
for i = 1:length(kv) % instead of test 2:2 Ameland
    num_tr  = m(i)-n(i)+1;
i
    for j = 1:num_tr
        altitude = nc_varget([OPT.ncpath2 OPT.Filename2],'altitude',[0 n(i)+j-2 0],[-1 1 -1]);
        
        for k = 1:length(time)
            if ~isempty(altitude(k,~isnan(altitude(k,:))))
                
            dummyint = findCrossings((cross(~isnan(altitude(k,:))))',altitude(k,~isnan(altitude(k,:))),...
                [cross(1) cross(end)],[2.5 2.5]); % dummy intersection, to be sure to evaluate MW after the dune part
            if ~isempty(dummyint)
                maxdummyint = max(dummyint);
                dummycross = intersect(cross(cross>=maxdummyint) , cross(~isnan(altitude(k,:))));
                dumindex   = ismember(cross,dummycross);
                
                if sum(dumindex) >= 2 % findCrossings needs at least 2 values
                    dumMHW   = min(findCrossings((dummycross)',altitude(k,dumindex),...
                    [cross(1) cross(end)],[MHWm(n(i)+j-1) MHWm(n(i)+j-1)]));

                end
                
            else
            dumMHW   = min(findCrossings((cross(~isnan(altitude(k,:))))',altitude(k,~isnan(altitude(k,:))),...
                [cross(1) cross(end)],[MHWm(n(i)+j-1) MHWm(n(i)+j-1)])); % it was max instead of min
            end
%             
%             if dumMHW < 0; dumMHW = []; end % and without this cycle 
            if ~isempty(dumMHW)
                MHWlatdum = interp1(cross',lat(n(i)+j-1,:),dumMHW);
                MHWlondum = interp1(cross',lon(n(i)+j-1,:),dumMHW);
                MHW_lat(n(i)+j-1,k) = MHWlatdum; MHW_lon(n(i)+j-1,k) = MHWlondum;
            else
                MHWlatdum = NaN; MHWlondum = NaN;
                MHW_lat(n(i)+j-1,k) = MHWlatdum; MHW_lon(n(i)+j-1,k) = MHWlondum;
            end
            
            end
        end
    end

%     MHW per each kustvak
    MHW_lat_kv  = MHW_lat(n(i):m(i),:);   MHW_lon_kv  = MHW_lon(n(i):m(i),:);

end

% MHW in RD coordinates
[MHW_x, MHW_y] = convertCoordinates(MHW_lon,MHW_lat,'CS1.code',4326,'CS2.code',28992);

% MHW in cross-shore coordinates wrt RSP
for j = 1:length(id) 
    j
    for k = 1:length(time)
        if (angle(j) >= 180 && angle(j) <= 360) % || (angle >= 90 && angle <= 180)
            if (MHW_x(j,k) - rsp_x(j)) >= 0
                MHW_rsp(j,k) = - sqrt((MHW_x(j,k) - rsp_x(j))^2 + (MHW_y(j,k) - rsp_y(j))^2);
            else
                MHW_rsp(j,k) = sqrt((MHW_x(j,k) - rsp_x(j))^2 + (MHW_y(j,k) - rsp_y(j))^2);
            end
        else
            if (MHW_x(j,k) - rsp_x(j)) >= 0
                MHW_rsp(j,k) = sqrt((MHW_x(j,k) - rsp_x(j))^2 + (MHW_y(j,k) - rsp_y(j))^2);
            else
                MHW_rsp(j,k) = - sqrt((MHW_x(j,k) - rsp_x(j))^2 + (MHW_y(j,k) - rsp_y(j))^2);
            end
        end
    end
end

save('MHW_rsp','MHW_rsp')

%% Evaluation of MLW

% counters
[kv, m, ~] = unique(areacode);
n          = vertcat(1, m(1:end-1)+1);
areakml    = cell(1,length(kv));

% definition
MLW_lat  = NaN(length(areacode),length(time)); MLW_lon  = NaN(length(areacode),length(time));
MLW_rsp  = NaN(length(areacode),length(time));

%% MLW
for i = 1:length(kv) % instead of test 2:2 Ameland
    num_tr  = m(i)-n(i)+1;
i
    for j = 1:num_tr
        altitude = nc_varget([OPT.ncpath2 OPT.Filename2],'altitude',[0 n(i)+j-2 0],[-1 1 -1]);
        
        for k = 1:length(time)
            if ~isempty(altitude(k,~isnan(altitude(k,:))))
            dummyint = findCrossings((cross(~isnan(altitude(k,:))))',altitude(k,~isnan(altitude(k,:))),...
                [cross(1) cross(end)],[2.5 2.5]); % dummy intersection, to be sure to evaluate MW after the dune part
            if ~isempty(dummyint)
                maxdummyint = max(dummyint);
                dummycross = intersect(cross(cross>=maxdummyint) , cross(~isnan(altitude(k,:))));
                dumindex   = ismember(cross,dummycross);
                
                if sum(dumindex) >= 2 % findCrossings needs at least 2 values           
                    dumMLW   = min(findCrossings((dummycross)',altitude(k,dumindex),...
                    [cross(1) cross(end)],[MLWm(n(i)+j-1) MLWm(n(i)+j-1)]));
                
                end
                
            else
            dumMLW   = min(findCrossings((cross(~isnan(altitude(k,:))))',altitude(k,~isnan(altitude(k,:))),...
                [cross(1) cross(end)],[MLWm(n(i)+j-1) MLWm(n(i)+j-1)])); % it was max instead of min
            end
%             
%             if dumMLW < 0; dumMLW = []; end % and without this cycle 
            if ~isempty(dumMLW)
                MLWlatdum = interp1(cross',lat(n(i)+j-1,:),dumMLW);
                MLWlondum = interp1(cross',lon(n(i)+j-1,:),dumMLW);
                MLW_lat(n(i)+j-1,k) = MLWlatdum; MLW_lon(n(i)+j-1,k) = MLWlondum;
            else
                MLWlatdum = NaN; MLWlondum = NaN;
                MLW_lat(n(i)+j-1,k) = MLWlatdum; MLW_lon(n(i)+j-1,k) = MLWlondum;
            end
            
            end
        end
    end
    
    
%     MLW per each kustvak
    MLW_lat_kv  = MLW_lat(n(i):m(i),:);   MLW_lon_kv  = MLW_lon(n(i):m(i),:);

end

% MLW in RD coordinates
[MLW_x, MLW_y] = convertCoordinates(MLW_lon,MLW_lat,'CS1.code',4326,'CS2.code',28992);

% MLW in cross-shore coordinates wrt RSP
for j = 1:length(id) 
    j
    for k = 1:length(time)
        if (angle(j) >= 180 && angle(j) <= 360) % || (angle >= 90 && angle <= 180)
            if (MLW_x(j,k) - rsp_x(j)) >= 0
                MLW_rsp(j,k) = - sqrt((MLW_x(j,k) - rsp_x(j))^2 + (MLW_y(j,k) - rsp_y(j))^2);
            else
                MLW_rsp(j,k) = sqrt((MLW_x(j,k) - rsp_x(j))^2 + (MLW_y(j,k) - rsp_y(j))^2);
            end
        else
            if (MLW_x(j,k) - rsp_x(j)) >= 0
                MLW_rsp(j,k) = sqrt((MLW_x(j,k) - rsp_x(j))^2 + (MLW_y(j,k) - rsp_y(j))^2);
            else
                MLW_rsp(j,k) = - sqrt((MLW_x(j,k) - rsp_x(j))^2 + (MLW_y(j,k) - rsp_y(j))^2);
            end
        end
    end
end

save('MLW_rsp','MLW_rsp')

