function [ttotal,Ttotal,fAtotal,fFtotal,fBtotal,fMtotal,Htot,HM,HB,HF] = CompMS(t,T,n_l,N_l)
% Calculate the Microstructure phase fractions based on the thermal history
% at a point for S690 steel

% Description of input and output
% Output
% ttotal = time vector that contains the selection of the thermal history
% that was used to compute the microstructure phases
% Ttotal = temperature vector that contains the selection of the thermal history
% that was used to compute the microstructure phases
% fAtotal = vector with austenite volume fraction at each timestep of the
% selection of the thermal history
% fFtotal = vector with ferrite volume fraction at each timestep of the
% selection of the thermal history
% fBtotal = vector with bainite volume fraction at each timestep of the
% selection of the thermal history
% fMtotal = vector with martensite volume fraction at each timestep of the
% selection of the thermal history
% Htot = hardness of the material in evaluated point
% HM = Martensite hardness in evaluated point
% HB = Bainite hardness in evaluated point
% HF = Ferrite hardness in evaluated point

% Input
% t = time vector of the thermal history of evaluated point
% T = Temperature vector of the thermal history of evaluated point
% n_l = layer number from which the evaluated point was taken
% N_l = total layers present in the part from which the point was taken

% Calculate the Microstructure phase fractions based on the thermal history
% at a point for S690 steel
%% constants/ parameters
% this parameter determines the maximum temperature difference between each
% timestep
T_step = 5;         % Temperature interval step

% TTT temperatures
TFS = 800.49;       % Ferrite start transformation temperature
TFF = 536.49;       % Ferrite finish transformation temperature
TBS = 533.55;       % Bainite start transformation temperature
TBF = 419.55;       % Bainite finish transformation temperature
TMS = 400;          % Martensite start transformation temperature
TMF = 280;          % Martensite finish transformation temperature

% Start and end temperatures during the heating cycle (heating rate dependent)
AC1H = 731;         % Starting temperature of the austenite formation
AC3H = 863;         % Finish temperature of the austenite formation

% TTT start and end volume fraction
fFs = 0.01;         % Ferrite start volume fraction
fFf = 0.99;         % Ferrite finish volume fraction
fBs = 0.01;         % Bainite start volume fraction
fBf = 0.99;         % Bainite finish volume fraction
fMs = 0.02;         % Martensite start volume fraction (guessed value, checked with Aravind)
fMf = 0.98;         % Martensite finish volume fraction (guessed value, checked with Aravind)

%% load TTT and feq data
load('TTT_Ferrite.mat')
load('TTT_Bainite.mat')
load('Ffeq_data.mat')

% Ferrite
TFs_inp = TTT_ferrite(:,1);
tFs_inp = TTT_ferrite(:,2);
tFf_inp = TTT_ferrite(:,3);
TFf_inp = TFs_inp;

% Bainite
TBs_inp = TTT_Bainite(:,1);
tBs_inp = TTT_Bainite(:,2);
tBf_inp = TTT_Bainite(:,3);
TBf_inp = TBs_inp; 

% Ferrite equilibrium data
TFeq_inp  = Ffeq_data(:,1);
fFeq_inp  = Ffeq_data(:,2);

% remove iterative data points of equilirium data
[TFeq_inp_unique, ia, ic] = unique(TFeq_inp);
fFeq_inp_unique = fFeq_inp(ia);
%%
T_bound_p = 300;    % temperature (degC) should be higher than this temperature to be a peak temperature
T_bound_h = 300;    % temperature (degC) should be lower than this temperature to be a heating temperature
Tcrit = [];         % critical temperature at which transformation fA->fF starts

% compute index, times and temperatures of the peaks and minimum
% temperatures
[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);

% tell how much cooling cycles are considered
disp(['#cc for layer',num2str(n_l),' = ', num2str(ncc)])

fontsize= 12;
figure
plot(t,T,'LineWidth', 1.5)
%yline(TFS,'k--','F_s','FontSize', fontsize);
yline(AC1H,'k--','AC1','FontSize', fontsize);
yline(AC3H,'k--','AC3','FontSize', fontsize);
% for i = 1:length(tpeak)
% xline(tpeak(i),'k--');
% xline(theat(i),'k--');
% end
for i = 1:length(i_ccycle_start)
xline(t(i_ccycle_start(i)),'k--');
xline(t(i_ccycle_end(i)),'k--');
end
title(['Thermal history of layer ',num2str(n_l)], 'FontSize', fontsize)
xlabel('t (s)', 'FontSize', fontsize)
ylabel('T (\circ C)', 'FontSize', fontsize)
% xlim([0,510])
ax = gca;
ax.XAxis.FontSize = fontsize;
ax.YAxis.FontSize = fontsize;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%% 1st cooling cycle %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% cooling cycle 1
% Select the temperature of the first cooling cycle
Tex = T(i_ccycle_start(1):i_ccycle_end(1));  
tex = t(i_ccycle_start(1):i_ccycle_end(1));

%% compute the cooling rate of the cooling cycle and plot it against the temperature
CR=abs(diff(Tex)./diff(tex));
% figure
% plot(CR,Tex(1:end-1))
% set ( gca, 'xdir', 'reverse')
% title(['cooling rate over the temperature cycle 1 L',num2str(n_l)])
% xlabel('cooling rate {\circ}C/s')
% ylabel('Temperature {\circ}C')
% yline(TFS,'k--','F_s');
% yline(TBS,'k--','B_s');
% yline(TBF,'k--','M_s');
% ylim([100 900])

%% split cooling cycle into isothermal temperature intervals (and plot it)
[TF,TB,TM,Tint,tint,tint_shift,Tiso,tiso] = isoTint(Tex,tex,T_step,TFS,TFF,TBS,TBF,TMS,TMF,Tcrit);
%% load TTT data of Ferrite, and Bainite and compute n and k and obtain fFeq for each temperature
[tFs_int,tFf_int,fFeq,nF,kF,nB,kB,gamma] = comp_nkfeqgamma(Tint,TF,TB,tFs_inp,TFs_inp,tFf_inp,TFf_inp,tBs_inp,TBs_inp,tBf_inp,TBf_inp,TFeq_inp_unique,fFeq_inp_unique,TMS,TMF,fFs,fFf,fBs,fBf,fMs,fMf);
%% Compute volume fractions of cooling cycle 1
% initial conditions
fA10 = 1;
fF10 = 0;
fB10 = 0;
fM10  = 0;
frA = 0;

%JMAK and MK
[fA1,fF1,fB1,fM1,ksiF1,ksiB1] = coolingcycle(Tint,tint_shift,fA10,fF10,fB10,fM10,frA,kF,nF,kB,nB,fFeq,gamma,TFS,TFF,TBS,TBF,TMS,TMF);

% print table on workspace
varNames = {'f_A','f_F','f_B','f_M'};
disp(['Microstructure after 1st cooling cycle L',num2str(n_l)])
Table = table(fA1(end),fF1(end),fB1(end),fM1(end),'VariableNames',varNames)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% 1st heating cycle and 2nd cooling cycle %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ncc>1 % (ncc = # cooling cycles) we need to compute the heating and 2nd cooling cycle if ncc>=2
%% Heating cycle
TexH = T(i_hcycle_start(1):i_hcycle_end(1));
texH = t(i_hcycle_start(1):i_hcycle_end(1));

% Select the temperature of the second cooling cycle
Tex2 = T(i_ccycle_start(2):i_ccycle_end(2));  
tex2 = t(i_ccycle_start(2):i_ccycle_end(2));

% compute phase fractions during heating cycle
[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);
disp(['Tcrit for Layer ',num2str(n_l),' = ',num2str(Tcrit),' degrees Celsius'])

disp(['Microstructure after 1st heating cycle L',num2str(n_l)])
varNames = {'f_A','f_F','f_B','f_M'};
Table = table(fAH(end),fFH(end),fBH(end),fMH(end),'VariableNames',varNames)

%% compute the cooling rate of the cooling cycle and plot it against the
% temperature
% CR2=abs(diff(Tex2)./diff(tex2));
% figure
% plot(CR2,Tex2(1:end-1))
% set ( gca, 'xdir', 'reverse')
% title(['cooling rate over the temperature 2nd cooling cycle L',num2str(n_l)])
% xlabel('cooling rate {\circ}C/s')
% ylabel('Temperature {\circ}C')
% yline(TFS,'k--','F_s');
% yline(TBS,'k--','B_s');
% yline(TBF,'k--','M_s');
% ylim([100 900])

%% split cycle into isothermal temperature intervals
[TF2,TB2,TM2,Tint2,tint2,tint_shift2,Tiso2,tiso2] = isoTint(Tex2,tex2,T_step,TFS,TFF,TBS,TBF,TMS,TMF,Tcrit);

%% load TTT data of Ferrite, compute n and k and obtain fFeq for each temperature
[tFs_int2,tFf_int2,fFeq2,nF2,kF2,nB2,kB2,gamma] = comp_nkfeqgamma(Tint2,TF2,TB2,tFs_inp,TFs_inp,tFf_inp,TFf_inp,tBs_inp,TBs_inp,tBf_inp,TBf_inp,TFeq_inp_unique,fFeq_inp_unique,TMS,TMF,fFs,fFf,fBs,fBf,fMs,fMf);

%% Compute volume fractions of second cooling cycle
% initial conditions
fA20 = fAH(end);
fF20 = fFH(end);
fB20 = fBH(end);
fM20 = fMH(end);

disp(['Microstructure after 2nd cooling cycle L',num2str(n_l)])
[fA2,fF2,fB2,fM2,ksiF2,ksiB2] = coolingcycle(Tint2,tint_shift2,fA20,fF20,fB20,fM20,frA,kF2,nF2,kB2,nB2,fFeq2,gamma,TFS,TFF,TBS,TBF,TMS,TMF);

varNames = {'f_A','f_F','f_B','f_M'};
Table = table(fA2(end),fF2(end),fB2(end),fM2(end),'VariableNames',varNames)

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% 2nd heating cycle and 3rd cooling cycle %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ncc>2 % (ncc = # cooling cycles) we need to compute the 2nd heating and 3rd cooling cycle if ncc>= 3 
%% 2nd Heating cycle
TexH2 = T(i_hcycle_start(2):i_hcycle_end(2));
texH2 = t(i_hcycle_start(2):i_hcycle_end(2));

% Select the temperature of the third cooling cycle
Tex3 = T(i_ccycle_start(3):i_ccycle_end(3));  
tex3 = t(i_ccycle_start(3):i_ccycle_end(3));

% compute phase fractions during heating cycle
[fAH2,fFH2,fBH2,fMH2,ksiAH2,tintH2,TintH2,Tcrit] = heatingcycle(T_step,TexH2,texH2,Tex3,tex3,fA1,fF1,fB1,fM1,TFeq_inp_unique,fFeq_inp_unique,AC1H,AC3H,TFS);
disp(['Tcrit for L',num2str(n_l),' = ',num2str(Tcrit),' degrees Celsius'])

disp(['Microstructure after 2nd heating cycle L',num2str(n_l)])
varNames = {'f_A','f_F','f_B','f_M'};
Table = table(fAH2(end),fFH2(end),fBH2(end),fMH2(end),'VariableNames',varNames)


%% compute the cooling rate of the cooling cycle and plot it against the
% temperature
% CR3=abs(diff(Tex3)./diff(tex3));
% figure
% plot(CR3,Tex3(1:end-1))
% set ( gca, 'xdir', 'reverse')
% title(['cooling rate over the temperature 3rd cooling cycle L',num2str(n_l)])
% xlabel('cooling rate {\circ}C/s')
% ylabel('Temperature {\circ}C')
% yline(TFS,'k--','F_s');
% yline(TBS,'k--','B_s');
% yline(TBF,'k--','M_s');
% ylim([100 900])

%% split cycle into isothermal temperature intervals
[TF3,TB3,TM3,Tint3,tint3,tint_shift3,Tiso3,tiso3] = isoTint(Tex3,tex3,T_step,TFS,TFF,TBS,TBF,TMS,TMF,Tcrit);

%% load TTT data of Ferrite, compute n and k and obtain fFeq for each temperature
[tFs_int3,tFf_int3,fFeq3,nF3,kF3,nB3,kB3,gamma] = comp_nkfeqgamma(Tint3,TF3,TB3,tFs_inp,TFs_inp,tFf_inp,TFf_inp,tBs_inp,TBs_inp,tBf_inp,TBf_inp,TFeq_inp_unique,fFeq_inp_unique,TMS,TMF,fFs,fFf,fBs,fBf,fMs,fMf);

%% Compute volume fractions of second cooling cycle
% initial conditions
fA30 = fAH2(end);
fF30 = fFH2(end);
fB30 = fBH2(end);
fM30 = fMH2(end);

disp(['Microstructure after 3rd cooling cycle L',num2str(n_l)])
[fA3,fF3,fB3,fM3,ksiF3,ksiB3] = coolingcycle(Tint3,tint_shift3,fA30,fF30,fB30,fM30,frA,kF3,nF3,kB3,nB3,fFeq3,gamma,TFS,TFF,TBS,TBF,TMS,TMF);

varNames = {'f_A','f_F','f_B','f_M'};
Table = table(fA3(end),fF3(end),fB3(end),fM3(end),'VariableNames',varNames)

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% Combine all cycles of interest %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Combine all cycles of interest
if ncc == 1 % 1 cooling cycle
    Ttotal = Tint;
    ttotal = tint;
    fAtotal = fA1;
    fFtotal = fF1;
    fBtotal = fB1;
    fMtotal = fM1;
    
elseif ncc>1 && ncc <= 2 % 2 cooling cycles
    Ttotal = [Tint,TintH',Tint2];
    ttotal = [tint,tintH',tint2];
    fAtotal = [fA1, fAH, fA2];
    fFtotal = [fF1, fFH, fF2];
    fBtotal = [fB1, fBH, fB2];
    fMtotal = [fM1, fMH, fM2];
    
elseif ncc>2 && ncc <= 3 % 3 cooling cycles
    Ttotal = [Tint,TintH',Tint2,TintH2',Tint3];
    ttotal = [tint,tintH',tint2,tintH2',tint3];
    fAtotal = [fA1, fAH, fA2, fAH2, fA3];
    fFtotal = [fF1, fFH, fF2, fFH2, fF3];
    fBtotal = [fB1, fBH, fB2, fBH2, fB3];
    fMtotal = [fM1, fMH, fM2, fMH2, fM3];  
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%% Compute Hardness %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% compute the hardness of each phase
% MS composition in wt%
pC      = 0.08;
pSi     = 0.44;
pMn     = 1.7;
pNi     = 1.35;
pCr     = 0.23;
pMo     = 0.3;

% Find cooling rate at 700 degrees in degrees per hour
CRH = interp1(Tex(1:end-1),CR,700)*3600;

% Empirical hardness equation to compute hardness of each separate phase
HM = 127 + 949*pC + 27*pSi + 11*pMn + 8*pNi + 16*pCr + 21*log10(CRH);
HB = -323 + 185*pC + 330*pSi + 153*pMn + 65*pNi + 144*pCr + 191*pMo +...
    (89 + 53*pC - 55*pSi - 22*pMn - 10*pNi - 20*pCr - 33*pMo)*log10(CRH);
HF = 42 + 223*pC + 53*pSi + 30*pMn;

% Compute total hardness
Htot = fMtotal(end)*HM + fBtotal(end)*HB + fFtotal(end)*HF; 

%% If there is 1 cooling cycles
%% plot all the different cooling cycles in one figure
if ncc == 1
 figure
subplot(2,1,1)
plot(ttotal,Ttotal,'b')
title(['Thermal history L',num2str(n_l)])
xlabel('time (s)')
ylabel('Temperature ({\circ}C)')
hold on
plot(tiso,Tiso,'r')
plot(tint,Tint,'r*','MarkerSize',1.5)
yline(TFS,'k--','F_s');
yline(TBS,'k--','B_s');
yline(TBF,'k--','M_s');
legend('temperature','isothermal steps')
subplot(2,1,2);
plot(ttotal,fAtotal,ttotal,fFtotal,ttotal,fBtotal,ttotal,fMtotal,'LineWidth',1.5)
title(['Microstructure phase composition L',num2str(n_l)])
xlabel('time (s)')
ylabel('Volume fraction')
legend('Austenite','Ferrite','Bainite','Martensite')
ylim([0,1])   

%% If there are 2 cooling cycles
elseif ncc>1 && ncc<=2
%% plot all the different cooling cycles in one figure
figure
subplot(2,1,1)
plot(ttotal,Ttotal,'b')
title(['Thermal history L',num2str(n_l)])
xlabel('time (s)')
ylabel('Temperature ({\circ}C)')
hold on
plot(tiso,Tiso,'r')
plot(tint,Tint,'r*','MarkerSize',1.5)
plot(tiso2,Tiso2,'r')
plot(tint2,Tint2,'r*','MarkerSize',1.5)
yline(TFS,'k--','F_s');
yline(TBS,'k--','B_s');
yline(TBF,'k--','M_s');
legend('temperature','isothermal steps')
subplot(2,1,2);
plot(ttotal,fAtotal,ttotal,fFtotal,ttotal,fBtotal,ttotal,fMtotal,'LineWidth',1.5)
title(['Microstructure phase composition L',num2str(n_l)])
xlabel('time (s)')
ylabel('Volume fraction')
legend('Austenite','Ferrite','Bainite','Martensite')
ylim([0,1])

%% If there are 3 cooling cycles
elseif ncc>2 && ncc<=3
%% plot all the different cooling cycles in one figure
fontsize = 12;
figure
subplot(2,1,1)
plot(ttotal,Ttotal,'b')
title(['Thermal history L',num2str(n_l)], 'FontSize', fontsize)
xlabel('time (s)', 'FontSize', fontsize)
ylabel('Temperature ({\circ}C)', 'FontSize', fontsize)
hold on
plot(tiso,Tiso,'r')
plot(tint,Tint,'r*','MarkerSize',1.5)
plot(tiso2,Tiso2,'r')
plot(tint2,Tint2,'r*','MarkerSize',1.5)
plot(tiso3,Tiso3,'r')
plot(tint3,Tint3,'r*','MarkerSize',1.5)
yline(TFS,'k--','F_s');
yline(TBS,'k--','B_s');
yline(TBF,'k--','M_s');
legend('temperature','isothermal steps')
ax = gca;
ax.XAxis.FontSize = fontsize;
ax.YAxis.FontSize = fontsize;
subplot(2,1,2);
% plot(ttotal,(fAtotal+fFtotal+fBtotal+fMtotal),ttotal,fAtotal,ttotal,fFtotal,ttotal,fBtotal,ttotal,fMtotal,'LineWidth',1.5)
plot(ttotal,fAtotal,ttotal,fFtotal,ttotal,fBtotal,ttotal,fMtotal,'LineWidth',1.5)
title(['Microstructure phase composition L',num2str(n_l)], 'FontSize', fontsize)
xlabel('time (s)', 'FontSize', fontsize)
ylabel('Volume fraction', 'FontSize', fontsize)
% legend('ftot','Austenite','Ferrite','Bainite','Martensite')
legend('Austenite','Ferrite','Bainite','Martensite')
ylim([0,1])
ax = gca;
ax.XAxis.FontSize = fontsize;
ax.YAxis.FontSize = fontsize;


figure
subplot(2,1,1)
plot(ttotal,Ttotal,'b')
title(['Thermal history L',num2str(n_l)])
xlabel('time (s)')
ylabel('Temperature ({\circ}C)')
hold on
plot(tiso,Tiso,'r')
plot(tint,Tint,'r*','MarkerSize',1.5)
plot(tiso2,Tiso2,'r')
plot(tint2,Tint2,'r*','MarkerSize',1.5)
plot(tiso3,Tiso3,'r')
plot(tint3,Tint3,'r*','MarkerSize',1.5)
yline(TFS,'k--','F_s');
yline(TBS,'k--','B_s');
yline(TBF,'k--','M_s');
legend('temperature','isothermal steps')
subplot(2,1,2);
plot(ttotal,fAtotal,ttotal,fFtotal,ttotal,fBtotal,ttotal,fMtotal,'LineWidth',1.5)
title(['Microstructure phase composition L',num2str(n_l)])
xlabel('time (s)')
ylabel('Volume fraction')
legend('Austenite','Ferrite','Bainite','Martensite')
ylim([0,1])
end  
end