Use specific variable from index to plot answer of equation

3 views (last 30 days)
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')

Accepted Answer

Voss
Voss on 16 Oct 2024
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')
  10 Comments
Erin Summerlin-Donofrio
Erin Summerlin-Donofrio on 18 Oct 2024
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
Voss 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

Sign in to comment.

More Answers (1)

Walter Roberson
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 Surfaces, Volumes, and Polygons in Help Center and File Exchange

Products


Release

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!