function [fAH,fFH,fBH,fMH,ksiAH,tintH,TintH,Tcrit] = heatingcycle(T_step,TexH,texH,Tex2,tex2,fA1,fF1,fB1,fM1,TFeq_inp_unique,fFeq_inp_unique,AC1H,AC3H,TFS)
% Heating cycle computes the amount of austenite that is formed when the
% temperature of the thermal cycle exceeds AC1H

% This script computes the avrami coefficients based on austenite volume
% fraction measurements during a heating experiment with a heating rate of
% 20 degC/s instead of using a TTT diagram as in the cooling cycle. This
% was done as no other data was available

istep = 1; % the ith neighbouring point used to compute 'n' and 'k' as we do not have a CHT diagram

% shift data to dit current CCT diagram
AC1H = AC1H;
AC3H = AC3H;

TexHc = [Tex2(Tex2>=AC1H);Tex2(find(Tex2<=AC1H,1,'first'))];
texHc = [tex2(Tex2>=AC1H);tex2(find(Tex2<=AC1H,1,'first'))];

TpeakH = max(TexH);
TH = linspace(AC1H,TpeakH,ceil((TpeakH-AC1H)/T_step)+1)';
TintH = [TexH(TexH<AC1H);TH];
tintH = interp1(TexH,texH,TintH,'linear');
TintHc = linspace(TpeakH,AC1H,ceil((TpeakH-AC1H)/T_step)+1)';
tintHc = interp1(TexHc,texHc,TintHc,'linear');
%
%remove the first data point as it overlaps with the previous dataset
TintHc = TintHc(1:end);
tintHc = tintHc(1:end);

% heating austenite
load('fA_heating20.mat');
Tmeas = fA_heating20(:,1);
tmeas = fA_heating20(:,2);
fAmeas = fA_heating20(:,3);

% Convert the data to a stepwise function, as this prevents heavy
% fluctuations in the Avrami parameters

for i = 1:5:length(fAmeas)-5
    fAstep(i:i+5,1) = mean(fAmeas(i:i+5));
end
nmeas = zeros(length(Tmeas),1);
kmeas = zeros(length(Tmeas),1);
nstep = zeros(length(Tmeas),1);
kstep = zeros(length(Tmeas),1);

for i  = 2: length(fAmeas)-istep
    nmeas(i) =  log(log(1-fAmeas(i))/log(1-fAmeas(i+istep)))/log(tmeas(i)/tmeas(i+istep));
    kmeas(i) = -log(1-fAmeas(i))/(tmeas(i)^nmeas(i));
end
for i  = 2:5:length(fAstep)-5
    nstep(i:i+5) =  log(log(1-fAstep(i))/log(1-fAstep(i+5)))/log(tmeas(i)/tmeas(i+5));
    kstep(i:i+5) = -log(1-fAstep(i))/(tmeas(i)^nstep(i));
end

nstep(end+1-5:end) = nstep(end-5);
kstep(end+1-5:end) = kstep(end-5);
%% interpolate the k and n over the temperature
nH = 1*interp1(Tmeas,nstep,TintH,'linear');
kH = 1*interp1(Tmeas,kstep,TintH,'linear');
fAeqH = 1-interp1(TFeq_inp_unique,fFeq_inp_unique,TintH,'linear','extrap');

nHc = 1*interp1(Tmeas,nstep,TintHc,'linear');
kHc = 1*interp1(Tmeas,kstep,TintHc,'linear');
fAeqHc = 1-interp1(TFeq_inp_unique,fFeq_inp_unique,TintHc,'linear','extrap');
% combine both temperature vectors, this is done later on, as otherwise
% there are two time points for a given temperature which causes errors in
% the simulation

nH = [nH;nHc];
kH = [kH;kHc];
fAeqH = [fAeqH;fAeqHc];
TintH = [TintH;TintHc];
tintH = [tintH;tintHc];

for i=1:length(nH)
    if isnan(nH(i)) == 1 && TintH(i)<AC1H
        nH(i) = 0;
        kH(i) = 0;
    end
    if isnan(nH(i)) == 1 && TintH(i)>AC3H
        nH(i) = nH(i-1);
        kH(i) = kH(i-1);
    end
end

%% Compute volume fractions
% initial conditions
fAH0 = fA1(end);
fFH0 = fF1(end);
fBH0 = fB1(end);
fMH0  = fM1(end);
fAH = zeros(1,length(TintH));
fFH = zeros(1,length(TintH));
fBH = zeros(1,length(TintH));
fMH = zeros(1,length(TintH));
ftotAH = zeros(1,length(TintH));
ksiAH = zeros(1,length(TintH));
Tcrit = zeros(1,length(TintH));
i_crit = zeros(1,length(TintH));

fAH(1) = fAH0;
fFH(1) = fFH0;
fBH(1) = fBH0;
fMH(1) = fMH0;
ftotAH(1) = fAH0+fFH0+fBH0+fMH0;
if fAH(1) > fAeqH(1)
    ksiAH(1) = (-log((fAH(1) - ftotAH(1)*fAeqH(1))/(ftotAH(1)*(1-fAeqH(1))))/kH(1+1))^(1/nH(1+1));
else
    ksiAH(1) = (log(1-(fAH(1)/fAeqH(1+1))/ftotAH(1))/-kH(1+1))^(1/nH(1+1));
end

%% compute the amount of austenite that is transformed back
for i = 2 : length(TintH)-1
    % Compute Austenite
    if TintH(i) < AC1H % below AC1 the phases do not transform to austenite
        fAH(i) = fAH(i-1);
        fFH(i) = fFH(i-1);
        fBH(i) = fBH(i-1);
        fMH(i) = fMH(i-1);
        ftotAH(i) = fAH(i) + fFH(i) + fBH(i) + fMH(i);
        
        if fAH(i) > fAeqH(i) % based on this statement we pick a ksi equation that results in a real bounded time value
            ksiAH(i) = (-log((fAH(i) - ftotAH(i)*fAeqH(i))/(ftotAH(i)*(1-fAeqH(i))))/kH(i+1))^(1/nH(i+1));
        else
            ksiAH(i) = (log(1-(fAH(i)/fAeqH(i+1))/ftotAH(i))/-kH(i+1))^(1/nH(i+1));
        end
    elseif TintH(i) > AC3H % above AC3 the amount of austenite should be 100%
        fAH(i) = 1;
        fFH(i) = ((1-fAH(i))/(1-fAH(i-1)))*fFH(i-1);
        fBH(i) = ((1-fAH(i))/(1-fAH(i-1)))*fBH(i-1);
        fMH(i) = ((1-fAH(i))/(1-fAH(i-1)))*fMH(i-1);
        ftotAH(i) = fAH(i) + fFH(i) + fBH(i) + fMH(i); % all available phase that is present in this transformation
        ksiAH(i) = (-log((fAH(i) - ftotAH(i)*fAeqH(i))/(ftotAH(i)*(1-fAeqH(i))))/kH(i+1))^(1/nH(i+1));
        
        if TintH(i) == AC1H
            if fAH(i) > fAeqH(i)  % based on this statement we pick a ksi equation that results in a real bounded time value
                ksiAH(i) = (-log((fAH(i) - ftotAH(i)*fAeqH(i))/(ftotAH(i)*(1-fAeqH(i))))/kH(i+1))^(1/nH(i+1));
            else
                ksiAH(i) = (log(1-(fAH(i)/fAeqH(i+1))/ftotAH(i))/-kH(i+1))^(1/nH(i+1));
            end
        end
        
    else % between AC1 and AC3 we use te JMAK equation
        fAH(i) = ftotAH(i-1)*(1-exp(-kH(i)*(ksiAH(i-1) + tintH(i)-tintH(i-1))^nH(i)))*fAeqH(i);
        if fAH(i) > 1*fAeqH(i)
            fAH(i) = fAeqH(i);
            Tcrit(i) = TintH(i);
            i_crit(i) = i
        end
        fFH(i) = ((1-fAH(i))/(1-fAH(i-1)))*fFH(i-1);
        fBH(i) = ((1-fAH(i))/(1-fAH(i-1)))*fBH(i-1);
        fMH(i) = ((1-fAH(i))/(1-fAH(i-1)))*fMH(i-1);
        ftotAH(i) = fAH(i) + fFH(i) + fBH(i) + fMH(i); % all available phase that is present in this transformation
        
        if fAH(i) > fAeqH(i) % based on this statement we pick a ksi equation that results in a real bounded time value
            ksiAH(i) = (-log((fAH(i) - ftotAH(i)*fAeqH(i))/(ftotAH(i)*(1-fAeqH(i))))/kH(i+1))^(1/nH(i+1));
        else
            ksiAH(i) = (log(1-(fAH(i)/fAeqH(i+1))/ftotAH(i))/-kH(i+1))^(1/nH(i+1));
        end
    end
end
fAH(end) = fAH(end-1);
fFH(end) = fFH(end-1);
fBH(end) = fBH(end-1);
fMH(end) = fMH(end-1);

Tcrit = nonzeros(Tcrit);
i_crit = nonzeros(i_crit);

% only consider the first index at which fA >= fAeq and select the data
% accordingly
if isempty(Tcrit) == 0
    Tcrit = Tcrit(1);
    i_crit = i_crit(1);
    tintH = tintH(1:i_crit);
    TintH = TintH(1:i_crit);
    fAH = fAH(1:i_crit);
    fFH = fFH(1:i_crit);
    fBH = fBH(1:i_crit);
    fMH = fMH(1:i_crit);
end

% When there is no Tcrit (T at which fA >= fAeq) Tcrit equals the last
% temperature of the interval

if isempty(Tcrit) == 1
    Tcrit = TintH(end);
end
end