function [fA,fF,fB,fM,ksiF,ksiB] = coolingcycle(Tint,tint_shift,fA0,fF0,fB0,fM0,frA,kF,nF,kB,nB,fFeq,gamma,TFS,TFF,TBS,TBF,TMS,TMF)
% Cooling cycle
% This function computes the volume fractions of the phases that transforms
% from austenite during a cooling cycle. This is done with the JMAK model and
% Koistinen-Marburger model.

%% Compute volume fractions
% initial condition
ftot0 = fA0+fF0+fB0;
fA = zeros(1,length(Tint));
fF = zeros(1,length(Tint));
fB = zeros(1,length(Tint));
fM = zeros(1,length(Tint));
ftotF = zeros(1,length(Tint));
ftotB = zeros(1,length(Tint));
ksiF = zeros(1,length(Tint));
ksiB = zeros(1,length(Tint));

fA(1) = fA0;
fF(1) = fF0;
fB(1) = fB0;
fM(1) = fM0;
ftotF(1) = fA0+fF0;
ftotB(1) = fA0+fB0;
fAeq = 1-fFeq;
fBeq = fFeq;
if fA0 == 1
    ksiF(1) = 0;
else
    if fA(1) > fAeq(1) % the first formula becomes imaginary otherwise
        ksiF(1) = (-log((fA(1) - ftotF(1)*fAeq(1))/(ftotF(1)*(1-fAeq(1))))/kF(1+1))^(1/nF(1+1));
    else
        ksiF(1) = (log(1-(fF(1)/fFeq(1+1))/ftotF(1))/-kF(1+1))^(1/nF(1+1));
    end
    if isreal(ksiF(1)) == 0
        disp(['imaginary ksiF neglected at index:',num2str(1)])
        ksiF(1) = real(ksiF(1));
    end
end
ksiB(1) = 0;

%% apply JMAK and MK to compute ferrite, bainite and martensite content
for i = 2 : length(Tint)
    % Compute Ferrite content
    if Tint(i) <= TFF  % if the temperature is below the transformation finish temperature and the temperature is not increasing, the volume fraction stays constant
        fF(i) = fF(i-1);
        fA(i) = fA(i-1);
    else
        fF(i) = ftotF(i-1)*(1-exp(-kF(i)*(ksiF(i-1) + tint_shift(i)-tint_shift(i-1))^nF(i)))*fFeq(i);
        fB(i) = fB(i-1);
        fM(i) = fM(i-1);
        fA(i) = fA(i-1)-(fF(i)-fF(i-1))-(fB(i)-fB(i-1))-(fM(i)-fM(i-1));
        ftotF(i) = fA(i) + fF(i); % all available phase that is present in this transformation
        ksiF(i) = (-log((fA(i) - ftotF(i)*fAeq(i))/(ftotF(i)*(1-fAeq(i))))/kF(i+1))^(1/nF(i+1));
        if isreal(ksiF(i)) == 0
            disp(['imaginary ksiF neglected at index:',num2str(i)])
            ksiF(i) = real(ksiF(i));
        end
        
    end
    % Compute Bainite content
    if Tint(i) == TBS
        if kB(i-1) == 0 && nB(i-1) == 0 % otherwise we do not get results as 1/0 is infinity
            disp('kB(i-1) == 0 && nB(i-1) == 0')
            kB(i-1) = kB(i);
            nB(i-1) = nB(i);
        end
        ftotB(i-1) = fA(i-1) + fB(i-1);
        ksiB(i-1) = (-log((fA(i-1) - ftotB(i-1)*fAeq(i-1))/(ftotB(i-1)*(1-fAeq(i-1))))/kB(i))^(1/nB(i));
    end
    
    if Tint(i) > TBS || Tint(i) <= TBF % if the temperature is above the starting temperature the volume fraction is equal to the volume fraction at the previous timestep
        fB(i) = fB(i-1);
        ftotB(i) = fA(i) + fB(i);
    else
        fB(i) = ftotB(i-1)*(1-exp(-kB(i)*(ksiB(i-1) + tint_shift(i)-tint_shift(i-1))^nB(i)))*fBeq(i);
        fM(i) = fM(i-1);
        fA(i) = fA(i-1)-(fF(i)-fF(i-1))-(fB(i)-fB(i-1))-(fM(i)-fM(i-1));
        ftotB(i) = fA(i) + fB(i);
        ksiB(i) = (-log((fA(i) - ftotB(i)*fAeq(i))/(ftotB(i)*(1-fAeq(i))))/kB(i+1))^(1/nB(i+1));
        if isreal(ksiB(i)) == 0
            disp(['imaginary ksiB neglected at index:',num2str(i)])
            ksiB(i) = real(ksiB(i));
        end
    end
    
    if Tint(i) == TMS
        i0M = i;
    end
    
    % Compute martensite content
    if Tint(i) >= TMS && (Tint(i)-Tint(i-1)) < 0 % if the temperature is above the starting temperature and the temperature is decreasing the volume fraction is equal to the volume fraction at the previous timestep
        fM(i) = fM(i-1);
    else
        fM(i) = fM(i0M-1) + (fA(i0M))*(1-exp(-gamma*(TMS-Tint(i))));
        fA(i) = fA(i-1)-(fM(i)-fM(i-1))-(fF(i)-fF(i-1))-(fB(i)-fB(i-1));
    end
end
end