Use specific variable from index to plot answer of equation
Show older comments
I have multiple equations using variables that are ranges of values, and some that aren't. Essentially I would like help indexing the variables properly so I can calculate the equations using the full range of variables and then also pull out specific numbers from the range and re-run the equations using those, without having to rewrite the equations.
Ultimately I want to run the calcultion for all possible values of delta_rfO and delta_rfH, but plot only the values that use T=573.15, B=0, and the full range of Q using indexing, for example. Is there a way to set this up with an index so I can choose the values (or ranges of values withing a range) in the plot command?
I tried in this line, specifing the value in Q from row 1 column 9
"plot(delta_rfO(Q(1,9)),delta_rfH(Q(1,9)))" obviously this doesn't work but it is an example of what I would like to do.
clc
clear
T=linspace(573.15,1173.15,100)% Temperatures used in calculation
n1=length(T)
Beta=linspace(0,1,100)%Variable relating cation content of mineral
n2=length(Beta)
delta_winO=0
delta_riO=6
delta_winH=0
delta_riH=-60
Q=linspace(0,3,100)%Ratio value for exponential
n3=length(Q)
f_w=0.2
C_H=0.0045%Concentration value
C_O=0.5%Concentration value
QoH=(f_w+(1-f_w)*C_H)%oxygen concentration
QoO=(f_w+(1-f_w)*C_O)%oxygen concentration
Delta_rwO=(2.91-0.76.*Beta) .* (10^6./T(:,1).^2) - 3.41 - 0.41.*Beta; %Fractionation plag-water
Delta_rwH=9.3 * (10^6./T(:,1).^2) - 61.9;
delta_rfO = Delta_rwO + delta_winO + (delta_riO - Delta_rwO - delta_winO).*exp(-Q./QoO) %rock final
delta_rfH = Delta_rwH + delta_winH + (delta_riH - Delta_rwH - delta_winH).*exp(-Q./QoH) %rock final
plot(delta_rfO,delta_rfH,'b')
plot(delta_rfO(Q(1,9)),delta_rfH(Q(1,9)))
hold on
Delta_rwO=(2.91-0.76.*Beta) .* (10^6./T(:,end).^2) - 3.41 - 0.41.*Beta; %Fractionation plag-water
Delta_rwH=9.3 * (10^6./T(:,end).^2) - 61.9;
delta_rfO = Delta_rwO + delta_winO + (delta_riO - Delta_rwO - delta_winO).*exp(-Q./QoO) %rock final
delta_rfH = Delta_rwH + delta_winH + (delta_riH - Delta_rwH - delta_winH).*exp(-Q./QoH) %rock final
plot(delta_rfO,delta_rfH,'b')
Accepted Answer
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')
10 Comments
Erin Summerlin-Donofrio
on 16 Oct 2024
Voss
on 16 Oct 2024
You're welcome! Any questions, let me know. Otherwise, please Accept this answer. Thanks!
Erin Summerlin-Donofrio
on 17 Oct 2024
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Calculation all the same as before %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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]);
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);
Here's a surface as requested, using the delta_rfH and delta_rfO values for all T and for the 1st Beta value (i.e., index 1 in the 3rd dimension of delta_rfO). If you want to use a different Beta value, you can find its index Bidx as shown in the last part of my answer, and use that Bidx here instead of 1.
% plot a surface
figure
Bidx = 1;
surf(delta_rfO(:,:,Bidx),delta_rfH,repmat(Q,NT,1),'EdgeColor','none')
xlabel('delta\_rfO')
ylabel('delta\_rfH')
zlabel('Q')
Thank you so much! That wil definitely be helpful for visualizing things on my end. For presentations though it might be clearer if I show it in 2D. Is there a similar method for plotting a shape/fill in 2D between the area of two curves and a straight line? For example the areas in blue and red here? I've plotted contours of equal Q as the dotted lines. Is there a way to fill in the shape made by the two blue curves and the Q(end) line? I am trying to show the separate fields as well as the overlap of phase space.
% % %Water Rock Ratio Open System Calculation (after Taylor, 1977 and Albarede,1995)
% %
%Setting up arrays
NT1 = 100;
NQ1 = 75;
NB1 = 50;
T1 = linspace(573.15,1173.15,NT1); % Temperatures used in calculation
Q1 = linspace(0,3,NQ1); % Variable relating cation content of mineral
Beta1 = linspace(0,1,NB1);% Ratio value for exponential
T1 = T1.';
%T is a column vector, Q is a row vector, and Beta is a vector that goes in
%the third dimension
Beta1 = permute(Beta1,[1 3 2]);
whos T1 Beta1 Q1
delta_winO = -4.5;
delta_riO = 6;
delta_winH = -26;
delta_riH = -85;
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.*Beta1) .* (10^6./T1.^2) - 3.41 - 0.41.*Beta1; % Fractionation plag-water
Delta_rwH = 9.3 * (10^6./T1(:).^2) - 61.9;
delta_rfO = Delta_rwO + delta_winO + (delta_riO - Delta_rwO - delta_winO).*exp(-Q1./QoO);
delta_rfH = Delta_rwH + delta_winH + (delta_riH - Delta_rwH - delta_winH).*exp(-Q1./QoH);
whos delta_rf*
hold on
% 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).
%Plotting albite-anorthite phase space for high T and low T, all Q
%Albite low T
plot(delta_rfO(1,:,1), delta_rfH(1,:), 'b','DisplayName',sprintf('T1=%g, all Q, Beta1=%g',T1(1), Beta1(1)));
%Albite high T
plot(delta_rfO(end,:,1), delta_rfH(end,:),'b','DisplayName',sprintf('T1=%g, all Q, Beta1=%g',T1(end),Beta1(1)));
%Contouring Water-Rock Ratios
%Q=0.125,0.25,0.5,1,infinity
plot([delta_rfO(1,4,1),delta_rfO(end,4,1)],[delta_rfH(1,4),delta_rfH(end,4)],'b', LineStyle=':' ,LineWidth=1)
plot([delta_rfO(1,7,1),delta_rfO(end,7,1)],[delta_rfH(1,7),delta_rfH(end,7)],'b', LineStyle=':' ,LineWidth=1)
plot([delta_rfO(1,13,1),delta_rfO(end,13,1)],[delta_rfH(1,13),delta_rfH(end,13)],'b', LineStyle=':' ,LineWidth=1)
plot([delta_rfO(1,26,1),delta_rfO(end,26,1)],[delta_rfH(1,26),delta_rfH(end,26)],'b', LineStyle=':' ,LineWidth=1)
plot([delta_rfO(1,end,1),delta_rfO(end,end,1)],[delta_rfH(1,end),delta_rfH(end,end)],'b')
%Anorthite low T
plot(delta_rfO(1,:,end), delta_rfH(1,:), 'r','DisplayName',sprintf('T1=%g, all Q, Beta1=%g',T1(1), Beta1(1)));
%Anorthite high T
plot(delta_rfO(end,:,end), delta_rfH(end,:),'r','DisplayName',sprintf('T1=%g, all Q, Beta1=%g',T1(end),Beta1(1)));
%Contouring Water-Rock Ratios
%Q=0.125,0.25,0.5,1,infinity
plot([delta_rfO(1,4,end),delta_rfO(end,4,end)],[delta_rfH(1,4),delta_rfH(end,4)],'r', LineStyle=':' ,LineWidth=1)
plot([delta_rfO(1,7,end),delta_rfO(end,7,end)],[delta_rfH(1,7),delta_rfH(end,7)],'r', LineStyle=':' ,LineWidth=1)
plot([delta_rfO(1,13,end),delta_rfO(end,13,end)],[delta_rfH(1,13),delta_rfH(end,13)],'r', LineStyle=':' ,LineWidth=1)
plot([delta_rfO(1,26,end),delta_rfO(end,26,end)],[delta_rfH(1,26),delta_rfH(end,26)],'r', LineStyle=':' ,LineWidth=1)
plot([delta_rfO(1,end,end),delta_rfO(end,end,end)],[delta_rfH(1,end),delta_rfH(end,end)],'r')
hold on
% %
% % %Water Rock Ratio Open System Calculation (after Taylor, 1977 and Albarede,1995)
% %
%Setting up arrays
NT = 100;
NQ = 75;
NB = 50;
T = linspace(573.15,1173.15,NT); % Temperatures used in calculation
Q = linspace(0,3,NQ); % Variable relating cation content of mineral
Beta = linspace(0,1,NB);% Ratio value for exponential
T = T.';
%T is a column vector, Q is a row vector, and Beta is a vector that goes in
%the third dimension
Beta = permute(Beta,[1 3 2]);
whos T Beta Q
delta_winO = 0;
delta_riO = 6;
delta_winH = 0;
delta_riH = -50;
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*
hold on
% 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).
%Plotting albite-anorthite phase space for high T and low T, all Q
%Albite low T
plot(delta_rfO(1,:,1), delta_rfH(1,:), 'b','DisplayName',sprintf('T=%g, all Q, Beta=%g',T(1), Beta(1)));
%Albite high T
plot(delta_rfO(end,:,1), delta_rfH(end,:),'b','DisplayName',sprintf('T=%g, all Q, Beta=%g',T(end),Beta(1)));
%Contouring Water-Rock Ratios
%Q=0.125,0.25,0.5,1,infinity
plot([delta_rfO(1,4,1),delta_rfO(end,4,1)],[delta_rfH(1,4),delta_rfH(end,4)],'b', LineStyle=':' ,LineWidth=1)
plot([delta_rfO(1,7,1),delta_rfO(end,7,1)],[delta_rfH(1,7),delta_rfH(end,7)],'b', LineStyle=':' ,LineWidth=1)
plot([delta_rfO(1,13,1),delta_rfO(end,13,1)],[delta_rfH(1,13),delta_rfH(end,13)],'b', LineStyle=':' ,LineWidth=1)
plot([delta_rfO(1,26,1),delta_rfO(end,26,1)],[delta_rfH(1,26),delta_rfH(end,26)],'b', LineStyle=':' ,LineWidth=1)
plot([delta_rfO(1,end,1),delta_rfO(end,end,1)],[delta_rfH(1,end),delta_rfH(end,end)],'b')
%Anorthite low T
plot(delta_rfO(1,:,end), delta_rfH(1,:), 'r','DisplayName',sprintf('T=%g, all Q, Beta=%g',T(1), Beta(1)));
%Anorthite high T
plot(delta_rfO(end,:,end), delta_rfH(end,:),'r','DisplayName',sprintf('T=%g, all Q, Beta=%g',T(end),Beta(1)));
%Contouring Water-Rock Ratios
%Q=0.125,0.25,0.5,1,infinity
plot([delta_rfO(1,4,end),delta_rfO(end,4,end)],[delta_rfH(1,4),delta_rfH(end,4)],'r', LineStyle=':' ,LineWidth=1)
plot([delta_rfO(1,7,end),delta_rfO(end,7,end)],[delta_rfH(1,7),delta_rfH(end,7)],'r', LineStyle=':' ,LineWidth=1)
plot([delta_rfO(1,13,end),delta_rfO(end,13,end)],[delta_rfH(1,13),delta_rfH(end,13)],'r', LineStyle=':' ,LineWidth=1)
plot([delta_rfO(1,26,end),delta_rfO(end,26,end)],[delta_rfH(1,26),delta_rfH(end,26)],'r', LineStyle=':' ,LineWidth=1)
plot([delta_rfO(1,end,end),delta_rfO(end,end,end)],[delta_rfH(1,end),delta_rfH(end,end)],'r')
xlabel('delta_rfO','fontweight','bold','fontsize',20)
ylabel('delta_rfH','fontweight','bold','fontsize',20)
box on
%Setting up arrays
NT = 100;
NQ = 75;
NB = 50;
T = linspace(573.15,1173.15,NT); % Temperatures used in calculation
Q = linspace(0,3,NQ); % Variable relating cation content of mineral
Beta = linspace(0,1,NB);% Ratio value for exponential
T = T.';
Beta = permute(Beta,[1 3 2]);
delta_winO = -4.5;
delta_riO = 6;
delta_winH = -26;
delta_riH = -85;
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);
figure
hold on
% plot stuff:
plot_stuff(delta_rfO,delta_rfH)
% change some variables:
delta_winO = 0;
delta_riO = 6;
delta_winH = 0;
delta_riH = -50;
% recalculate only what needs to be recalculated:
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);
% plot the new results too:
plot_stuff(delta_rfO,delta_rfH)
xlabel('\Delta_{rf}^O','fontweight','bold','fontsize',20)
ylabel('\Delta_{rf}^H','fontweight','bold','fontsize',20)
box on
function plot_stuff(delta_rfO,delta_rfH)
fill( ...
[delta_rfO(1,:,end) delta_rfO([1 end],end,end).' delta_rfO(end,end:-1:1,end)], ...
[delta_rfH(1,:) delta_rfH([1 end],end).' delta_rfH(end,end:-1:1)], ...
'r','EdgeColor','r','FaceAlpha',0.25)
fill( ...
[delta_rfO(1,:,1) delta_rfO([1 end],end,1).' delta_rfO(end,end:-1:1,1)], ...
[delta_rfH(1,:) delta_rfH([1 end],end).' delta_rfH(end,end:-1:1)], ...
'b','EdgeColor','b','FaceAlpha',0.25)
idx = [4 7 13 26];
%Contouring Water-Rock Ratios
%Q=0.125,0.25,0.5,1,infinity
for ii = 1:numel(idx)
plot(delta_rfO([1 end],idx(ii),1),delta_rfH([1 end],idx(ii)),'b', LineStyle=':' ,LineWidth=1)
end
plot(delta_rfO([1 end],end,1),delta_rfH([1 end],end),'b')
%Contouring Water-Rock Ratios
%Q=0.125,0.25,0.5,1,infinity
for ii = 1:numel(idx)
plot(delta_rfO([1 end],idx(ii),end),delta_rfH([1 end],idx(ii)),'r', LineStyle=':' ,LineWidth=1)
end
plot(delta_rfO([1 end],end,end),delta_rfH([1 end],end),'r')
end
So I made some changes to your code besides just adding the fills. These include:
1. Simplifying expressions like
[delta_rfO(1,4,1),delta_rfO(end,4,1)]
to
delta_rfO([1 end],4,1)
That is, you can concatenate indices inside an indexing expression, rather than concatenating the results of multiple indexing expressions, and get the same result. Another (simpler) example:
x = [1 2 3];
[x(1) x(end)]
x([1 end])
2. Replacing sequences of similar lines of code with for loops, as in
plot([delta_rfO(1,4,1),delta_rfO(end,4,1)],[delta_rfH(1,4),delta_rfH(end,4)],'b', LineStyle=':' ,LineWidth=1)
plot([delta_rfO(1,7,1),delta_rfO(end,7,1)],[delta_rfH(1,7),delta_rfH(end,7)],'b', LineStyle=':' ,LineWidth=1)
plot([delta_rfO(1,13,1),delta_rfO(end,13,1)],[delta_rfH(1,13),delta_rfH(end,13)],'b', LineStyle=':' ,LineWidth=1)
plot([delta_rfO(1,26,1),delta_rfO(end,26,1)],[delta_rfH(1,26),delta_rfH(end,26)],'b', LineStyle=':' ,LineWidth=1)
replaced by
idx = [4 7 13 26];
for ii = 1:numel(idx)
plot([delta_rfO(1,idx(ii),1),delta_rfO(end,idx(ii),1)],[delta_rfH(1,idx(ii)),delta_rfH(end,idx(ii))],'b', LineStyle=':' ,LineWidth=1)
end
but with the additional simplification listed in point #1.
3. Avoiding duplicate code by creating a plotting function that does the plotting, which can be called twice (or as many times as you need), rather than copy-pasting the lines of code that do the plotting.
4. Avoiding duplicate code by avoiding defining redundant variables T1, Q1, etc., which are the same as T, Q, etc. When delta_winO, delta_riO, etc., need to change for the second set of calculations, you can get away with only recalculating delta_rfO and delta_rfH, since those are the only things that depend on the new variable values. Thus, you only have to re-run 2 lines of calculation code instead of the whole script. (And you could make a function that does the calculation as well, but since it was only two lines of code, I left it alone.)
5. Removing the plotting of the lines comprising the edges of the patches created by fill, and instead specifying the EdgeColor in the fill call. If needed, you can modify the LineStyle and LineWidth properties in that same fill call to modify the properties of the edge.
Erin Summerlin-Donofrio
on 17 Oct 2024
Erin Summerlin-Donofrio
on 18 Oct 2024
It makes sense to me to modify the plotting function to include creation of the texts. An example is below. Note that I had to include the variable Q as an input to the function. And I made the figure bigger, so the texts don't overlap too much. Adjust as needed.
%Setting up arrays
NT = 100;
NQ = 75;
NB = 50;
T = linspace(573.15,1173.15,NT); % Temperatures used in calculation
Q = linspace(0,3,NQ); % Variable relating cation content of mineral
Beta = linspace(0,1,NB);% Ratio value for exponential
T = T.';
Beta = permute(Beta,[1 3 2]);
delta_winO = -4.5;
delta_riO = 6;
delta_winH = -26;
delta_riH = -85;
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);
figure('Position',[10 10 1000 1000])
hold on
% plot stuff:
plot_stuff(delta_rfO,delta_rfH,Q)
% change some variables:
delta_winO = 0;
delta_riO = 6;
delta_winH = 0;
delta_riH = -50;
% recalculate only what needs to be recalculated:
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);
% plot the new results too:
plot_stuff(delta_rfO,delta_rfH,Q)
xlabel('\Delta_{rf}^O','fontweight','bold','fontsize',20)
ylabel('\Delta_{rf}^H','fontweight','bold','fontsize',20)
box on
function plot_stuff(delta_rfO,delta_rfH,Q)
fill( ...
[delta_rfO(1,:,end) delta_rfO([1 end],end,end).' delta_rfO(end,end:-1:1,end)], ...
[delta_rfH(1,:) delta_rfH([1 end],end).' delta_rfH(end,end:-1:1)], ...
'r','EdgeColor','r','FaceAlpha',0.25)
fill( ...
[delta_rfO(1,:,1) delta_rfO([1 end],end,1).' delta_rfO(end,end:-1:1,1)], ...
[delta_rfH(1,:) delta_rfH([1 end],end).' delta_rfH(end,end:-1:1)], ...
'b','EdgeColor','b','FaceAlpha',0.25)
idx = [4 7 13 26];
%Contouring Water-Rock Ratios
%Q=0.125,0.25,0.5,1,infinity
for ii = 1:numel(idx)
plot(delta_rfO([1 end],idx(ii),1),delta_rfH([1 end],idx(ii)),'b', LineStyle=':' ,LineWidth=1)
text(delta_rfO(end,idx(ii),1),delta_rfH(end,idx(ii)),sprintf('Q=%0.3f',Q(idx(ii))),'Color','b','HorizontalAlignment','center')
end
text(delta_rfO(end,end,1),delta_rfH(end,end),'Q=Inf','Color','b','HorizontalAlignment','center')
%Contouring Water-Rock Ratios
%Q=0.125,0.25,0.5,1,infinity
for ii = 1:numel(idx)
plot(delta_rfO([1 end],idx(ii),end),delta_rfH([1 end],idx(ii)),'r', LineStyle=':' ,LineWidth=1)
text(delta_rfO(1,idx(ii),end),delta_rfH(1,idx(ii)),sprintf('Q=%0.3f',Q(idx(ii))),'Color','r','HorizontalAlignment','center')
end
text(delta_rfO(1,end,end),delta_rfH(1,end),'Q=Inf','Color','r','HorizontalAlignment','right')
end
More Answers (1)
Walter Roberson
on 16 Oct 2024
plot(delta_rfO(Q(1,9)),delta_rfH(Q(1,9)))
You are trying to use Q(1,9) as a subscript to delta_rf0 but Q(1,9) does not happen to be a positive integer.
Probably what you want is
plot(delta_rfO(1,9),delta_rfH(1,9), 'o')
The 'o' will cause a 'o' marker to be drawn at the single point designated. Without the 'o' option, plot() will only draw lines when there are two adjacent finite points, and since you are only drawing a single point there are no two adjacent finite points so no line would be drawn.
Categories
Find more on Creating and Concatenating Matrices in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
