function [ipeak,tpeak,Tpeak,iheat,theat,Theat,dTheat,i_ccycle_start,i_ccycle_end,i_hcycle_start,i_hcycle_end,ncc] = peaktemperatures(t,T,T_bound_p,T_bound_h,AC3H,AC1H,N_l)
%This function gives the time and temperatures of the peaks of a thermal history
% also it looks at the times at which the sample gets heated up
j1 = 1;
% j2 = 1;
ipeak = zeros(N_l,1);
Tpeak = zeros(N_l,1);
tpeak = zeros(N_l,1);
iheat = zeros(N_l,1);
Theat = zeros(N_l,1);
theat = zeros(N_l,1);
dTheat = zeros(N_l,1);

for i = 2:length(t)-1
    if T(i)>T_bound_p && (T(i)-T(i-1))>0 && (T(i)-T(i+1))>0
        ipeak(j1) = i;
        Tpeak(j1) = T(i);
        tpeak(j1) = t(i);
        j1 = j1+1;
    end
end

% remove the empty rows
ipeak = nonzeros(ipeak);
Tpeak = nonzeros(Tpeak);
tpeak = nonzeros(tpeak);

% Find the heating points (minima) inbetween the peaks. The first heating
% point is not of interest, as the temperature of the first peak always
% surpasses the melting point

for i = 1:length(ipeak)-1
    iheat(i+1) = find(T==min(T(ipeak(i):ipeak(i+1))));
    Theat(i+1) = T(iheat(i+1));
    theat(i+1) = t(iheat(i+1));
    dTheat(iheat(i+1)) = ((T(iheat(i+1)+1)-T(iheat(i+1)))/(t(iheat(i+1)+1)-t(iheat(i+1))));
end

% remove the empty rows, introduce row for the first peak ( the exact value
% of this peak does not matter, but there needs to be a row to simulate the
% presence of a first heating peak. The last 'heating peak' is the end of
% the simulation
iheat = [1;nonzeros(iheat);length(t)];
Theat = [1;nonzeros(Theat);T(end)];
theat = [1;nonzeros(theat);t(end)];

% find the interesting domains to study, so how many cooling and heating
% cycles are there
i_start = find(Tpeak>AC3H,1,'last'); % So the analysed domain of the thermal history starts at the last peak in which T>AC3
i_end = find(Tpeak>AC1H,1,'last'); % The analysed domain of the thermal history ends after the last cooling cycle with T>AC1H
%i_end = find(Tpeak<AC1H,1,'first'); % The analysed domain of the thermal history ends just before a reheating cycle in which T<AC1

ncycles = i_end-i_start+1
disp(['ncycles = ', ncycles])
if ncycles>1
    i_ccycle_start = zeros(ncycles,1);
    i_ccycle_end   = zeros(ncycles,1);
    i_hcycle_start = zeros(ncycles-1,1);
    i_hcycle_end   = zeros(ncycles-1,1);
    
    for i = 1: ncycles
        i_ccycle_start(i) = ipeak(i_start+i-1);
        i_ccycle_end(i) = iheat(i_start+i);
        if i < ncycles % we always have one heatcycle less than cooling cycles
            i_hcycle_start(i) = iheat(i_start+i);
            i_hcycle_end(i) = ipeak(i_start+i);
        end
    end
else
    i_ccycle_start = ipeak(i_start);
    i_ccycle_end   = iheat(i_start+1);
    i_hcycle_start = [];
    i_hcycle_end   = [];
end

ncc = length(i_ccycle_start);   % number of cooling cycles of interest
end