% Simulates glass volume change in cooling/reheating cycles using the free
% volume theory of relaxation
% JJ Hu, 09/25/2015

clear all;
CTE_SCL = 1E-4;     % CTE of supercooled liquid in /K
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

% Cooling
V_0 = 1;            % Initial volume (normalized)
V_f0 = 0;           % Free volume of initial state (normalized)
T_start = 1000;     % Starting temperature in K
T_end = 600;        % 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 = 100;         % Cooling rate in K/s
dt = abs(dT / dT_dt);   % Time step in s

ii = 1;
V_f(ii) = V_f0;
V(ii) = V_0;

for T = T_start : dT : T_end
    tau_re = tau0 * exp(Delta_E_a/k_B/T);
    dV = V_f(ii) * (exp(-dt/tau_re) - 1);
    % dV = - dt * V_f(ii) / tau_re;
    V(ii + 1) = V(ii) + dV;
    V_f(ii + 1) = V_f(ii) + dV - dT * CTE_SCL;
    ii = ii + 1;
end
hold on
plot([T_start : dT : T_end + dT], V, 'b')

% Reheating
V_0 = V(ii);    % Initial volume (normalized)
V_f0 = V_f(ii); % Free volume of initial state (normalized)
clear V V_f;
T_start = T_end;    % 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 = 1;          % Reheating rate in K/s
dt = abs(dT / dT_dt);	% Time step in s

ii = 1;
V_f(ii) = V_f0;
V(ii) = V_0;

for T = T_start : dT : T_end
    tau_re = tau0 * exp(Delta_E_a/k_B/T);
    dV = V_f(ii) * (exp(-dt/tau_re) - 1);
    % dV = - dt * V_f(ii) / tau_re;
    V(ii + 1) = V(ii) + dV;
    V_f(ii + 1) = V_f(ii) + dV - dT * CTE_SCL;
    ii = ii + 1;
end

hold on
plot([T_start : dT : T_end + dT], V, 'b')
plot([T_start : dT : T_end + dT], 1 - CTE_SCL*abs(dT)*(length(V)-1) : CTE_SCL*abs(dT) : 1, 'k--')
set(gca, 'box', 'on')
xlabel('Temperature (K)', 'Fontsize', 12)
ylabel('Normalized volume', 'Fontsize', 12)
axis([min(T_start, T_end) max(T_start, T_end) 1-(1-min(V))*1.25 1])