%%Simulation Code for Temperature Programmed Desorption in MATLAB by Omar A.
%%Abdelrahman, UMass Amherst, Released 2021
clear % Clear any Previous Values in Workspace
clc % Clear Command Window
Tspan = linspace (300,700,1000); %Temperature span of TPD experiment in degrees Kelvin
Theta0 = 1; % Initial fractional coverage of surface adsorbate, must be between 0 and 1
n = 1; % Desorption Order, 0,1 and 2 valid options
dSact = 44; % Entropy of Activation J/mol.K
dHact = 165; % Enthalpy of Activation kJ/mol
R = 8.314; % Universal Gas Constant J/mol.K
B = 10; % Temperature Ramp Rate K/min
kb = 1.38064852 * 10^-23; % Boltzmann Constant m^2 kg s^-2 K^-1
h = 6.62607004 * 10^-34; % Planck's constant m^2 kg / s
[T,Theta] = ode15s(@(T,Theta) -(kb*T/h)*(60/B)*exp(dSact/R)*exp(-dHact*1000/R/T)*Theta^n, Tspan, Theta0); % Polanyi-Wigner Equation,pre-exponential factor and activation energy assumed to not be a function of coverage
X = [T,Theta]; % Temperature and Coverage Matrix
Rate = (kb*X(:,1)./h).*X(:,2).^n*(60/B)*exp(dSact/R).*exp(-dHact*1000/R./X(:,1)); % Rate of desorption in mol Product/mol Site/s
[M,I]=max(Rate); % Identifying maximum rate and index value
Tpeakx = [T(I),T(I)]; % Desorption Peak Temperature Location, for plotting purposes
Tpeaky = [0,M*1.5]; % Height of Desorption Peak, for plotting purposes

%%Plotting TPD results 
plot(T,Rate,'k','LineWidth',2) % Plot the rate of desorption 
title('Temperature Programmed Desorption') % Plot Title
xlabel ('Temperature (K)','fontweight','bold') % X-axis label
ylabel ('Rate of Desorption (mol Product/mol site/s)','fontweight','bold') % Y-axis label
ylim([-0.0001 inf]) % Y Axis limits
hold on 
plot (Tpeakx, Tpeaky,'--k') % Plotting the Peak Temperature Location of the Rate of Desorption
txt = ['\leftarrow T_p_e_a_k = ' num2str(round(T(I))) ' K'];
text(T(I),M*1.25,txt)
hold on
yyaxis right % Adding Secondary Y-axis for Coverage
set(gca, 'YColor','b') % Secondary Y-axis Color
ylabel ('Coverage (\Theta)','color','b','fontweight','bold') % Secondary Y-axis label
ylim([-0.0001 1.0001]) % Secondary Y-axis limits
plot(X(:,1),X(:,2),'--b','LineWidth',2) % Plotting Surface Fractional Coverage