If you want to calculate delta_rfO for all combinations of the specified values of T, Q, and Beta, you can permute (or transpose as the case may be) your T, Q, and Beta vectors appropriately; then the delta_rfO calculation will result in a 3-dimensional array.
Here I've decided to have T correspond to the first dimension of delta_rfO, Q to the second dimension, and Beta to the third dimension, because delta_rfH depends on T and Q only and doesn't depend on Beta. But you can set it up however you like. (Also, I changed some of the lengths of the T, Q, and Beta vectors so that they're all different because it's more clear to see which dimension of the result corresponds to which vector when they're not all the same length.)
NT = 100;
NQ = 75;
NB = 50;
T = linspace(573.15,1173.15,NT); % Temperatures used in calculation
Beta = linspace(0,1,NB); % Variable relating cation content of mineral
Q = linspace(0,3,NQ); % Ratio value for exponential
T = T.';
Beta = permute(Beta,[1 3 2]);
T is a column vector, Q is a row vector, and Beta is a vector that goes in the third dimension (a "page vector"?):
whos T Beta Q
delta_winO = 0;
delta_riO = 6;
delta_winH = 0;
delta_riH = -60;
f_w = 0.2;
C_H = 0.0045;
C_O = 0.5;
QoH = f_w+(1-f_w)*C_H;
QoO = f_w+(1-f_w)*C_O;
Delta_rwO = (2.91-0.76.*Beta) .* (10^6./T.^2) - 3.41 - 0.41.*Beta; % Fractionation plag-water
Delta_rwH = 9.3 * (10^6./T(:).^2) - 61.9;
delta_rfO = Delta_rwO + delta_winO + (delta_riO - Delta_rwO - delta_winO).*exp(-Q./QoO);
delta_rfH = Delta_rwH + delta_winH + (delta_riH - Delta_rwH - delta_winH).*exp(-Q./QoH);
whos delta_rf*
You can see that delta_rfO is now a 3D array of size NT-by-NQ-by-NB; thus it contains a result for all combinations of the specified T, Q, and Beta values. Similarly delta_rfH is a NT-by-NQ matrix, containing a result for each combination of T and Q (it doesn't depend on Beta).
You can plot slices from these against each other, using indexing. For example:
figure
hold on
plot(delta_rfO(1,:,1), delta_rfH(1,:), 'b','DisplayName',sprintf('T=%g, all Q, Beta=%g',T(1), Beta(1)));
plot(delta_rfO(end,:,1), delta_rfH(end,:),'r','DisplayName',sprintf('T=%g, all Q, Beta=%g',T(end),Beta(1)));
plot(delta_rfO(end,:,end),delta_rfH(end,:),'g','DisplayName',sprintf('T=%g, all Q, Beta=%g',T(end),Beta(end)));
plot(delta_rfO(1,:,end), delta_rfH(1,:), 'k','DisplayName',sprintf('T=%g, all Q, Beta=%g',T(1), Beta(end)));
legend('Location','NorthWest')
xlabel('delta\_rfO')
ylabel('delta\_rfH')
Another example:
figure
hold on
plot(delta_rfO(:,9,1), delta_rfH(:,9), 'b','DisplayName',sprintf('all T, Q=%g, Beta=%g',Q(9), Beta(1)));
plot(delta_rfO(:,end,1), delta_rfH(:,end),'r','DisplayName',sprintf('all T, Q=%g, Beta=%g',Q(end),Beta(1)));
plot(delta_rfO(:,end,end),delta_rfH(:,end),'g','DisplayName',sprintf('all T, Q=%g, Beta=%g',Q(end),Beta(end)));
plot(delta_rfO(:,9,end), delta_rfH(:,9), 'k','DisplayName',sprintf('all T, Q=%g, Beta=%g',Q(9), Beta(end)));
legend('Location','NorthWest')
xlabel('delta\_rfO')
ylabel('delta\_rfH')
If you want to be able to specify a range of values and plot from that - say based on a min and max T and a specific given Beta value - you can do something like:
Tmin = 1080;
Tmax = 1120;
Beta_given = 0.4;
Tidx = find(T >= Tmin & T <= Tmax)
[~,Bidx] = min(abs(Beta-Beta_given))
figure
hold on
for ii = 1:numel(Tidx)
for jj = 1:numel(Bidx)
plot(delta_rfO(Tidx(ii),:,Bidx(jj)),delta_rfH(Tidx(ii),:), ...
'DisplayName',sprintf('T=%g, Beta=%g',T(Tidx(ii)),Beta(Bidx(jj))));
end
end
legend('Location','Best')
xlabel('delta\_rfO')
ylabel('delta\_rfH')
