Use specific variable from index to plot answer of equation

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

T = 1×100

573.1500 579.2106 585.2712 591.3318 597.3924 603.4530 609.5136 615.5742 621.6348 627.6955 633.7561 639.8167 645.8773 651.9379 657.9985 664.0591 670.1197 676.1803 682.2409 688.3015 694.3621 700.4227 706.4833 712.5439 718.6045 724.6652 730.7258 736.7864 742.8470 748.9076

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

n1=length(T)

n1 = 100

Beta=linspace(0,1,100)%Variable relating cation content of mineral

Beta = 1×100

0 0.0101 0.0202 0.0303 0.0404 0.0505 0.0606 0.0707 0.0808 0.0909 0.1010 0.1111 0.1212 0.1313 0.1414 0.1515 0.1616 0.1717 0.1818 0.1919 0.2020 0.2121 0.2222 0.2323 0.2424 0.2525 0.2626 0.2727 0.2828 0.2929

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

n2=length(Beta)

n2 = 100

delta_winO=0

delta_winO = 0

delta_riO=6

delta_riO = 6

delta_winH=0

delta_winH = 0

delta_riH=-60

delta_riH = -60

Q=linspace(0,3,100)%Ratio value for exponential

Q = 1×100

0 0.0303 0.0606 0.0909 0.1212 0.1515 0.1818 0.2121 0.2424 0.2727 0.3030 0.3333 0.3636 0.3939 0.4242 0.4545 0.4848 0.5152 0.5455 0.5758 0.6061 0.6364 0.6667 0.6970 0.7273 0.7576 0.7879 0.8182 0.8485 0.8788

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

n3=length(Q)

n3 = 100

f_w=0.2

f_w = 0.2000

C_H=0.0045%Concentration value

C_H = 0.0045

C_O=0.5%Concentration value

C_O = 0.5000

QoH=(f_w+(1-f_w)*C_H)%oxygen concentration

QoH = 0.2036

QoO=(f_w+(1-f_w)*C_O)%oxygen concentration

QoO = 0.6000

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_rfO = 1×100

6.0000 5.9715 5.9417 5.9108 5.8790 5.8462 5.8127 5.7784 5.7435 5.7081 5.6722 5.6359 5.5993 5.5623 5.5252 5.4878 5.4503 5.4127 5.3750 5.3372 5.2995 5.2617 5.2240 5.1863 5.1488 5.1113 5.0739 5.0366 4.9995 4.9625

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

delta_rfH = Delta_rwH + delta_winH + (delta_riH - Delta_rwH - delta_winH).*exp(-Q./QoH) %rock final

delta_rfH = 1×100

-60.0000 -56.3477 -53.2005 -50.4885 -48.1516 -46.1378 -44.4025 -42.9072 -41.6187 -40.5083 -39.5515 -38.7271 -38.0166 -37.4044 -36.8768 -36.4222 -36.0305 -35.6930 -35.4021 -35.1514 -34.9354 -34.7493 -34.5889 -34.4507 -34.3317 -34.2290 -34.1406 -34.0644 -33.9987 -33.9422

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

plot(delta_rfO,delta_rfH,'b')

plot(delta_rfO(Q(1,9)),delta_rfH(Q(1,9)))

Array indices must be positive integers or logical values.

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')

Answer 1

[Voss]

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

  Name        Size              Bytes  Class     Attributes

  Beta        1x1x50              400  double              
  Q           1x75                600  double              
  T         100x1                 800  double              

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*

  Name             Size                 Bytes  Class     Attributes

  delta_rfH      100x75                 60000  double              
  delta_rfO      100x75x50            3000000  double              

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')

(You could also potentially plot surfaces.)

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)

Tidx = 7×1

85 86 87 88 89 90 91

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[~,Bidx] = min(abs(Beta-Beta_given))

Bidx = 21

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')

[Erin Summerlin-Donofrio]

Yes! This is exactly what I was trying to accomplish! Thank you!

[Voss]

You're welcome! Any questions, let me know. Otherwise, please Accept this answer. Thanks!

[Erin Summerlin-Donofrio]

I would actually like to plot a surface too, with delta_rfO as x, delta_rfH as y, and Q as z. How would I specify the indexing to do this? I assume I would need to pick the same T and Beta for all calculations?

[Voss]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 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')

[Erin Summerlin-Donofrio]

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

  Name         Size              Bytes  Class     Attributes

  Beta1        1x1x50              400  double              
  Q1           1x75                600  double              
  T1         100x1                 800  double              

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*

  Name             Size                 Bytes  Class     Attributes

  delta_rfH      100x75                 60000  double              
  delta_rfO      100x75x50            3000000  double              

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

  Name        Size              Bytes  Class     Attributes

  Beta        1x1x50              400  double              
  Q           1x75                600  double              
  T         100x1                 800  double              

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*

  Name             Size                 Bytes  Class     Attributes

  delta_rfH      100x75                 60000  double              
  delta_rfO      100x75x50            3000000  double              

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

[Voss]

%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

[Voss]

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)]

ans = 1×2

1 3

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

x([1 end])

ans = 1×2

1 3

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

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]

AMAZING!!! Thank you!!

The indexing information is so helpful! I really appreciate it. I don't know why I didn't think to just restate the variables and the plotting function.

[Erin Summerlin-Donofrio]

One last question I think! How would I put a label/annotation/text at the end of a chosen plotted line, for example, at either end of one of the Q=0.25 lines? Is that all modified within the function? I've not worked with those before. I attempted to put a text line after the first plot (line 26 or 37) command but it only labeled one of the lines. How would I code the annotation in the for loop?

[Voss]

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.

https://www.mathworks.com/help/matlab/ref/matlab.graphics.primitive.text-properties.html

%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

Answer 2

[Walter Roberson]

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.