% MATLAB Code file in support of "Controls on Iron-Redox State in Martian
% Magmas Quantified by Mössbauer Spectroscopy, Colorimetric Wet Chemistry,
% and XANES Spectroscopy," by S. P. Aithala, R. A. Lange, and M. M. Hirschmann

% Submitted to JGR-Planets 

% Code for regressing linear models for Eqns. 9 and 10 from
% experimental data from the present study, Righter et al. (2013), and Matzen
% et al. (2022)

% Must be run with regressiondata3.4.csv in the same file path

testData = readtable('regressiondata3.4.csv');

% Fe3+/FeT data from experiments, same data recast as log(Fe3+/Fe2+),
% uncertainties, weights (1/sigma^2), and experimental conditions
Fer_Tot = table2array(testData(1:41,2));
sigmaFer_Tot = table2array(testData(1:41,7));
logFeFe = table2array(testData(1:41,5));
weights = table2array(testData(1:41,9));
logfO2 = table2array(testData(1:41,10));
invT = table2array(testData(1:41,13));

%compositional data
XSiO2 = table2array(testData(1:41,33));
XTiO2 = table2array(testData(1:41,34));
XAl2O3 = table2array(testData(1:41,35));
XMgO = table2array(testData(1:41,36));
XCaO = table2array(testData(1:41,37));
XNa2O = table2array(testData(1:41,38));
XK2O = table2array(testData(1:41,39));
XP2O5 = table2array(testData(1:41,40));
XSiXAl = table2array(testData(1:41,41));
XSiXMg = table2array(testData(1:41,42));
nFeO = table2array(testData(1:41,26));
nTotal = table2array(testData(1:41,32));
XFeO = nFeO./nTotal;

%%
% stepwise linear regeression for predicting log(XFe3+/XFe2+) from logfO2, inverse
% temperature, and compositional term (Eq. 9)

tot_mat = [XSiO2 XTiO2 XAl2O3 XMgO XCaO XNa2O XK2O XP2O5 XSiXAl...
    XSiXMg XFeO];
test = [logfO2 invT tot_mat];
E9 = stepwiselm(test, logFeFe, 'Weights', weights)
E9anova = anova(E9, 'summary')

% Calculate Fe3+/FeT ratios predicted by the Eq. 9 model and compare it to
% observed Fe3+/FeT and quantify agreement with reduced chi-squared (rcs)
% statistic
E9FF = 10.^(E9.Fitted);
E9Fer_tot = E9FF./(1+E9FF);
E9wsr = sigmaFer_Tot.^-2.*(E9Fer_tot-Fer_Tot).^2;
E9_rcs = sum(E9wsr)/E9anova.DF(3)

%%
% Iterative regression of all possible linear regression models with 11
% compositional components
tot_mat = [XSiO2 XTiO2 XAl2O3 XMgO XCaO XNa2O XK2O XP2O5 XSiXAl...
   XSiXMg XFeO];
P = dec2bin(0:(2^11)-1)-'0';
i=1;
for(i=1:2048)
    mat = tot_mat.*P(i,:);
    test2 = [logfO2 invT mat];
    E10 = fitlm(test2, logFeFe, 'Weights', weights);
    wsr(i) = sum((sigmaFer_Tot.^-2).*((10.^E10.Fitted)./(1+10.^E10.Fitted) - Fer_Tot).^2);
    redchisq(i) = wsr(i)./(39-3-sum(P(i,:)));
    k_val(i) = E10.Coefficients.Estimate(2);
    h_val(i) = E10.Coefficients.Estimate(3);
   
end
%determine the iteration of compositional components that has a reduced chi
%square between calculated and observed Fe3+/FeT closest to 1.
[x1 y1] = min(abs(1-redchisq))


%%
% Take the iteration of compositional components that optimizes the reduced
% chi-square between calculated and observed Fe3+/FeT to define Eq. 10
matfinal = tot_mat.*P(y1,:);
test3 = [logfO2 invT matfinal];

E10 = fitlm(test3, logFeFe, 'Weights', weights)
E10_anova = anova(E10, 'summary')

%%
% Generate figure 10 - some aspects of the plot are further modified using
% illustrative software to aid with visual clarity
figure 
%plot all possible combinations of k and h values
scatter(k_val, h_val, 3, 'b','filled', 'MarkerFaceAlpha', 0.2, 'MarkerEdgeAlpha', 0.2)
hold on

%plot k and h values for Eq. 9 and 10 with uncertainties
e = errorbar([E9.Coefficients.Estimate(2) E10.Coefficients.Estimate(2)]...
    ,[E9.Coefficients.Estimate(3) E10.Coefficients.Estimate(3)],...
    [E9.Coefficients.SE(3) E10.Coefficients.SE(3)],...
    [E9.Coefficients.SE(3) E10.Coefficients.SE(3)],...
    [E9.Coefficients.SE(2) E10.Coefficients.SE(2)],...
    [E9.Coefficients.SE(2) E10.Coefficients.SE(2)]...
    ,'s','MarkerFaceColor', 'r', 'MarkerEdgeColor', 'r')
e.Color = [0 0 0];
e.LineWidth = 2;
% plot place holders for k and h coefficients for Kress and Carmichael
% (1991), Jayasuriya et al. (2004), Borisov et al. (2018), and Righter et
% al. (2013)
scatter([.196 .1967 .207 .22], [11492/2.303 12420/2.303 4633.3...
3800/2.303],50,'MarkerFaceColor', [0 0 0], 'MarkerEdgeColor',[0 0 0])
xlim([0.18 0.23]);
k_valT = k_val.';
h_valT = h_val.';
% draw purple polygon
k = boundary(k_valT, h_valT,0);
hold on
plot(k_val(k), h_val(k),'LineWidth',1)
ax = gca;
box on
ax.BoxStyle = 'full';
ax.LineWidth = 1;
ax.TickDir = 'in';
ax.FontSize = 10;
ax.XMinorTick = 'on';
ax.YMinorTick = 'on';
xlabel('k (*logfO2)')
ylabel('h (*1/T)')

%% KC91
% Linear regression models for Table 4 based off Kress and Carmichael
% (1991)
%
KC91_mat=[logfO2 invT XAl2O3 XCaO XFeO XNa2O XK2O];
KC91mdl=fitlm(KC91_mat, logFeFe, 'Weights', weights.^2)
KC91anova = anova(KC91mdl, 'summary');

KC91FF = 10.^(KC91mdl.Fitted);
KC91Fer_tot = KC91FF./(1+KC91FF);
KC91wsr = sigmaFer_Tot.^-2.*(KC91Fer_tot-Fer_Tot).^2;
KC91_rcs = sum(KC91wsr)/KC91anova.DF(3)

%% B18
% Linear regression models for Table 4 based off Borisov et al. (2018)

b18_mat = [logfO2 invT XSiO2 XTiO2 XMgO XCaO XNa2O XK2O XP2O5 XSiXAl XSiXMg];
b18mdl=fitlm(b18_mat, logFeFe, 'Weights', weights.^2)
b18anova=anova(b18mdl, 'summary');

b18FF = 10.^(b18mdl.Fitted);
b18Fer_tot = b18FF./(1+b18FF);
b18wsr = sigmaFer_Tot.^-2.*(b18Fer_tot-Fer_Tot).^2;
b18_rcs = sum(b18wsr)/b18anova.DF(3)
 
%% J04
% Linear regression models for Table 4 based off Jaysuriya et al. (2004)

J04_mat = [logfO2 invT XMgO XCaO XNa2O XK2O XAl2O3 XP2O5 XFeO];
J04mdl=fitlm(J04_mat, logFeFe, 'Weights', weights.^2);
J04anova = anova(J04mdl, 'summary');

J04FF = 10.^(J04mdl.Fitted);
J04Fer_tot = J04FF./(1+J04FF);
J04wsr = sigmaFer_Tot.^-2.*(J04Fer_tot-Fer_Tot).^2; 
J04_rcs = sum(J04wsr)/J04anova.DF(3)

%% R13
% Linear regression models for Table 4 based off Righter et al. (2013)

R13_mat = [logfO2 invT XFeO XAl2O3 XCaO XNa2O XK2O XP2O5];
R13mdl=fitlm(R13_mat, logFeFe, 'Weights', weights.^2);
R13anova = anova(R13mdl, 'summary');

R13FF = 10.^(R13mdl.Fitted);
R13Fer_tot = R13FF./(1+R13FF);
R13wsr = sigmaFer_Tot.^-2.*(R13Fer_tot-Fer_Tot).^2;
R13_rcs = sum(R13wsr)/R13anova.DF(3)