% Here is a Matlab script to plot the two Wattmeters in stress-grain size
% space from Tokle and Hirth (2021). 

clc; close all; clear all
%% Deformation Conditions
T = 900+273; % Temperature in K
%fug = 5030;  % Water fugacity in MPa
P = 1.5e9; % pressure in Pa

%% recrystallized grain size
d = 10;

%% Basic Parameters
sigma = linspace(1,5000,1e3); % stress from 1 to 5000 MPa

R = 8.314; % gas constant
 
%% Flow law parameters (Tokle et al. 2019)
n4 = 4;
Q4 = 125000;
r4 = 1;
A4 = 1.75e-12;

n3 = 3;
Q3 = 115000;
r3 = 1.2;
A3 = 8e-13;

%% Wattmeter parameters
gamma = 1; 
beta = 1;
c = 3.1416; 

lamda=0.015;  % RMS grain size
%lamda=0.04;  % arithmetic mean grain size
%% Water fugacity (Shinevar et al. 2015)
AA1 = 5521;
AA2 = 31.28e3;
AA3 = -2.009e-5;

f = AA1*exp(-(AA2+P.*AA3)/(R*T));

%% Grain growth parameters
% Tokle Hirth  2020
p = 3;

Qg = 134e3;
rg = 1.38;
Ag = 0.261;

Kg = Ag*(f.^rg);


%% wattmeter n = 4

edotdis4 = A4.*sigma.^n4.*f.^r4.*exp(-Q4./R./T);
WM4 = (((Kg.*exp(-Qg./R./T).*p.^(-1).*c.*gamma)./((sigma).*lamda.*beta.*edotdis4)).^(1./(1+p)));

%% wattmeter n = 3

edotdis3 = A3.*sigma.^n3.*f.^r3.*exp(-Q3./R./T);
WM3 = (((Kg.*exp(-Qg./R./T).*p.^(-1).*c.*gamma)./((sigma).*lamda.*beta.*edotdis3)).^(1./(1+p)));

%% Solving for stress (in MPa) with a known grain size (in um)
% n = 4
k4 = (n4+1)./(p+1);
r_4 = (rg-r4)./(p+1);
Q_4 = (Qg-Q4)./(p+1);
C4 = ((Ag.*c.*gamma)./(A4.*beta.*lamda.*p)).^(1./(p+1));

W4_stress = (d./(C4.*(f.^r_4).*exp(-Q_4./(R.*T)))).^(-(p+1)./(n4+1))

% n = 3
k3 = (n3+1)./(p+1);
r_3 = (rg-r3)./(p+1);
Q_3 = (Qg-Q3)./(p+1);
C3 = ((Ag.*c.*gamma)./(A3.*beta.*lamda.*p)).^(1./(p+1));

W3_stress = (d./(C3.*(f.^r_3).*exp(-Q_3./(R.*T)))).^(-(p+1)./(n3+1))

% Stipp and Tullis (2003)
ST03_stress = (d/3631)^(-1/1.26)

% Cross et al (2017) [Sliding RMS]
C17_stress = (d/10^4.22)^(-1/1.59)

%% Plotting
figure
hold all

% Cross (2017) sliding piezometer
C17 = 10^4.22*(sigma.^-1.59);
% plot(sigma,C17,'y','Linewidth',2);

% Stipp and Tullis (2003) piezometer
ST = 3631*(sigma.^-1.26);
plot(sigma,ST,'k','Linewidth',2);

% Wattmeter n = 4
plot(sigma,WM4,'r','LineWidth',2);

% Wattmeter n = 3
plot(sigma,WM3,'b','LineWidth',2);

%% Figure details
box on
pbaspect([1 1 1])
set(gca,'Xscale','log');
set(gca,'Yscale','log');
axis([10 5000 0.5 100]);
xticks([10 30 100 300 1000 3000]);
yticks([1 3 10 30 100 ]);
set(gca,'fontsize',18);
set(gca,'TickLength',[0.02, 0.01])
xlabel('Equivalent Stress (MPa)')
ylabel('Grain size (um)')
grid on