% Simulates DSC curve (heat capacity) of glass undergoing glass transition
% JJ Hu, 11/27/2015

clear all;
C_SCL = 1E-3;       % Heat capacity of supercooled liquid normalized to
% enthalpy of the starting state H_0/1K
Delta_E_a = 3;      % Activation energy in eV
k_B = 8.62E-5;      % Boltzmann constant in eV/K
tau0 = 1E-18;       % Pre-exponential factor for relaxation time in s

% Heating
H_0 = 1;            % Initial enthalpy (normalized)
H_e0 = 0.2;         % Excess enthalpy of initial state with respect to the
% supercooled liquid state (normalized)
T_start = 600;      % Starting temperature in K
T_end = 1000;       % End point temperature in K
dT = 0.1 * sign(T_end - T_start);
% Temperature step in modeling in K; takes positive value for heating and
% negative value for cooling
dT_dt = 10;         % Heating rate in K/s
dt = abs(dT / dT_dt);   % Time step in s

% Solving enthalpy change during temperature change using the first-order
% relaxation model in a step-wise manner
ii = 1;
H_e(ii) = H_e0;
H(ii) = H_0;
for T = T_start : dT : T_end
    tau_re = tau0 * exp(Delta_E_a/k_B/T);
    dH = H_e(ii) * (exp(-dt/tau_re) - 1);
    H(ii + 1) = H(ii) + dH;
    H_e(ii + 1) = H_e(ii) + dH - dT * C_SCL;
    ii = ii + 1;
end

% Plot enthalpy of glass and supercooled liquid
% hold on
% plot(T_start : dT : T_end + dT, H, 'b')
% plot(T_start : dT : T_end + dT, H_0 - H_e0 : C_SCL*abs(dT) : H_0 - H_e0 + C_SCL*abs(dT)*(length(H)-1), 'k--')
% set(gca, 'box', 'on')
% xlabel('Temperature (K)', 'Fontsize', 12)
% ylabel('Normalized enthalpy', 'Fontsize', 12)
% axis([min(T_start, T_end) max(T_start, T_end) max(H)-(max(H)-min(H))*1.25 min(H)+(max(H)-min(H))*1.25])

% Calculate heat absorbed or released
% Assume a constant heating rate
dQ = diff(H);       % Heat flux, reference to a reference sample which has
% a constant heat capacity equaling to that of the glass in its glass state
% figure
hold on
plot(T_start : dT : T_end, dQ, 'r')
set(gca, 'box', 'on')
xlabel('Temperature (K)', 'Fontsize', 16)
ylabel('Heat flow', 'Fontsize', 16)