close all
clc
clear all

load Ht.mat
H = Htfirst(:,2);
t = Htfirst(:,1);

% H is the time series including the values of H
% t is the corresponding time vector
% Ref: 

Hmax = std(H)/5;
Hmin = 0.01*Hmax;
amin = log10(Hmin);
amax = log10(Hmax);
da = (amax-amin)/25;
Veca = amin:da:amax;
VecH = 10.^Veca;
NH = length(VecH);

FPTsum = zeros(NH,NH);
Occurrence = zeros(NH,NH);

%% Characterization of the time serie: ascending or descending zone

H = [0; H];                                                                 % First index added for derivative determination
t = [t(1); t];
N = length(t);

type = zeros(1,N);
nextindex = zeros(1,N);

for i = 2:N-1
    if (H(i-1)<H(i) && H(i)<H(i+1)) || (H(i-1)>H(i) && H(i)<H(i+1))         % Ascending zone or min
        type(i) = 1;
    elseif (H(i-1)>H(i) && H(i)>H(i+1)) || (H(i-1)<H(i) && H(i)>H(i+1))     % Descend or max
        type(i) = 2;
    end
end

% Determine in ascending zone the index of the next maximum of H and in descending zone the
% index of the next higher value of H

start1 = 0;                                                                 % start1 = 1 if in type 1 zone
for i = 1:N-1
    if type(i) == 1 && start1 == 0                                          % entering an ascending zone
        start1 = 1;
        indexstart1 = i;                                                    % save of the first index of the ascending zone
        if i == N-1
            nextindex(i) = NaN;                                             % detection of th end of the time serie: NaN
        end
    elseif (type(i) == 2 && start1 == 1) || (i==N-1 && start1==1)           % entering a descending zone
        start1 = 0;
        nextindex(indexstart1:i-1) = i;                                     % this index is the end of the ascending zone
        if i == N-1
            nextindex(i) = NaN;
        end
    end
    
    if type(i) == 2                                                         % in a descending zone, looking for a next higher value
        reached = 0;                                                        % reach = 1 when a new higher value of H is met
        j = i+1;
        while reached == 0 && j<=N
            if H(j)>H(i)
                nextindex(i) = j;
                reached = 1;
            end
            if j==N-1 || j==N
                nextindex(i) = NaN;
            end
            j = j+1;
        end
    end
end
%% Determination of the main envelope: Presence registers the indexes of the values in H_ or 0 if not present

i0 = 2;
Hc = H(i0);                                                                 % Hc a fragment a H to be build
t1 = t(i0);                                                                 % t1 a fragment of t to be build
ilocal = i0;

Presence = zeros(1,N);
Presence(i0)=1;

while not(isnan(nextindex(ilocal)))                                         % until the end of the signal
    [Hclocal, t1local, inext, indexes] = local(ilocal,H,t,type,nextindex);
    Presence(indexes) = length(t1)+1:length(t1)+length(indexes);            % Presence registers at each index the index of that point in the
    % envelope if it is in it; index = 0 if it is not in the envelope
    Hc = [Hc; Hclocal];
    t1 = [t1; t1local];
    ilocal = inext;
end
Hcmain = Hc;
t1main = t1;

%%  Determination of the partial envelope and completion of the FPT chart

for i0 = 2:N-1
    H0 = H(i0);
    if H0 <= Hmax                                                           % Hmax is excluded as that bin will include all the values of H(t) higher than Hmax
        
        a0 = log10(H0);
        indexH0 = max(ceil((a0-amin)/da),1);
        
        if Presence(i0)~=0                                                   % If the concerned H0 is part of the main envelope, we take the envelope from that point
            Hc = Hcmain(Presence(i0):end);
            t1 = t1main(Presence(i0):end);
        else                                                                % If it is not, we reconstruct the envelope until we reach the envelope again
            Hc = H(i0);                                                                 
            t1 = t(i0);                                                                 
            ilocal = i0;
            Main = 0;                                                       % We are not on the main envelope yet
            
            while not(isnan(nextindex(ilocal)))                             % until the end of the signal or until the main envelope is met
                [Hclocal, t1local, inext, ~] = local(ilocal,H,t,type,nextindex);
                Hc = [Hc; Hclocal];
                t1 = [t1; t1local];
                ilocal = inext;
                if Presence(inext)~=0                                       % when the main envelope is reached we take the main envelope from that point until the end
                    Hc = [Hc; Hcmain(Presence(ilocal)+1:end)];
                    t1 = [t1; t1main(Presence(ilocal)+1:end)];
                    ilocal = N-1;
                end
            end
        end
        t1 = t1-t(i0);                                                      % time vector shifted to start from zero
                
        DH = Hc-H0;
        Da = log10(DH);
        I = find(Hc <= Hmax);
        indexDH = max(ceil((Da-amin)/da),1);
        
        FPTsum(indexH0,indexDH(I)) = FPTsum(indexH0,indexDH(I)) + transpose(t1(I));    % first passage times are summed in FPT and counted in Occurrence to get the average at the end
        Occurrence(indexH0,indexDH(I)) = Occurrence(indexH0,indexDH(I)) + 1;
    end
end

FPTMean = FPTsum./Occurrence;

% minU = min(min(FPTMean));
% maxU = max(max(FPTMean));
% delta = maxU-minU;
% VecU = minU+0.25*delta:delta/20:maxU-0.25*delta;
% figure;hold on;[C h] = contour(Veca,Veca,transpose(FPTMean),VecU);xlabel('H_0');ylabel('\Delta H');title('t_1')



