Why my code find a second intersection of curve with line and not the first?
2 views (last 30 days)
Show older comments
I am trying to find a intersection of very complex curve's equation with line parallel with X-axis, but Matlab gives me a result which is suitable for the second intersection. I need first one.
Plot code for better understand is also included. Magenta line should reach and stop at first intersection (pH=6.4 not the second one at pH=13,38). What am I doing wrong?
Oh yes, and I need a exact value of X-coordinate of the first intersection and also a graphical solution as it can be seen.
Here is a MWE:
clc;
close all;
%There are input constants
cP=0.05;
cN=0.05;
cMg=0.005;
logb0=2.88;
logb1=1.17;
logb2=4.85;
logb3=2.58;
logKs1=-12.6;
pKa1=-2.143;
pKa2=-7.205;
pKa3=-12.34;
pKb1NH3=-4.753355133;
%there is a range of x-axis (0-14)
pH=0:0.01:14;
%There are calculations of cH2PO4 cHPO4 cPO4, they needs to be calculated for final curve equation (they are variables of final curve equation and depend on pH value)
Ka1=power(10,pKa1);
Ka2=power(10,pKa2);
Ka3=power(10,pKa3);
cH2PO4=cP.*Ka1.*(power(10,-pH)).^(2)./((power(10,-pH).^(3))+Ka1.*(power(10,-pH).^(2))+Ka1.*Ka2.*(power(10,-pH).^(1))+Ka1.*Ka2.*Ka3);
cHPO4=cP.*Ka1.*Ka2.*power(10,-pH)./((power(10,-pH).^(3))+Ka1.*(power(10,-pH).^(2))+Ka1.*Ka2.*(power(10,-pH).^(1))+Ka1.*Ka2.*Ka3);
cPO4=cP.*Ka1.*Ka2.*Ka3./((power(10,-pH).^(3))+Ka1.*(power(10,-pH).^(2))+Ka1.*Ka2.*(power(10,-pH).^(1))+Ka1.*Ka2.*Ka3);
%There is calculations of cNH4, it have to be calculated for final curve equation (it is variable of final curve equation and depends on pH value also)
pOH=14-pH;
Kb1=power(10,pKb1NH3);
cNH4=cN.*Kb1./(power(10,-pOH)+Kb1);
%there are some partial calculations for the final curve equation
b0=power(10,logb0);
b1=power(10,logb1);
b2=power(10,logb2);
b3=power(10,logb3);
Ks1=power(10,logKs1);
cOH=power(10,pH-14);
%this is the final curve equation
logcNP=Ks1./(cPO4.*cNH4).*(1+b1.*cH2PO4+b0.*cHPO4+b2.*cPO4+b3*cOH);
%I am trying to calculate a intersection of the curve above with y value of 0.005 (cMg at the start of the code).
crossectionY=cMg
[crossectionY, indexAtCrossectionY] = min(abs(logcNP-crossectionY));
crossectionX = pH(indexAtCrossectionY(1))
%This is plot of the final curve equation (black) and calculated intersection with magenta line
semilogy(pH,logcNP,'-', 'color', 'k', 'MarkerFaceColor','w', 'MarkerEdgeColor', 'm');
box off;
hold on;
plot([0 crossectionX],[cMg cMg],'-', 'color', 'm', 'LineWidth',0.5);
box off;
hold on;
axis([0 14 0.000000000000000001 1000000000000000000]);
xlabel('pH (-)');
ylabel('c {(mol/dm^{3})}');
axis square;
grid on;
axis on;
box off;
hold off;
Than you for answers. :-)
0 Comments
Accepted Answer
Star Strider
on 14 May 2023
See if the lines oif copde I added do what you want —
%There are input constants
cP=0.05;
cN=0.05;
cMg=0.005;
logb0=2.88;
logb1=1.17;
logb2=4.85;
logb3=2.58;
logKs1=-12.6;
pKa1=-2.143;
pKa2=-7.205;
pKa3=-12.34;
pKb1NH3=-4.753355133;
%there is a range of x-axis (0-14)
pH=0:0.01:14;
%There are calculations of cH2PO4 cHPO4 cPO4, they needs to be calculated for final curve equation (they are variables of final curve equation and depend on pH value)
Ka1=power(10,pKa1);
Ka2=power(10,pKa2);
Ka3=power(10,pKa3);
cH2PO4=cP.*Ka1.*(power(10,-pH)).^(2)./((power(10,-pH).^(3))+Ka1.*(power(10,-pH).^(2))+Ka1.*Ka2.*(power(10,-pH).^(1))+Ka1.*Ka2.*Ka3);
cHPO4=cP.*Ka1.*Ka2.*power(10,-pH)./((power(10,-pH).^(3))+Ka1.*(power(10,-pH).^(2))+Ka1.*Ka2.*(power(10,-pH).^(1))+Ka1.*Ka2.*Ka3);
cPO4=cP.*Ka1.*Ka2.*Ka3./((power(10,-pH).^(3))+Ka1.*(power(10,-pH).^(2))+Ka1.*Ka2.*(power(10,-pH).^(1))+Ka1.*Ka2.*Ka3);
%There is calculations of cNH4, it have to be calculated for final curve equation (it is variable of final curve equation and depends on pH value also)
pOH=14-pH;
Kb1=power(10,pKb1NH3);
cNH4=cN.*Kb1./(power(10,-pOH)+Kb1);
%there are some partial calculations for the final curve equation
b0=power(10,logb0);
b1=power(10,logb1);
b2=power(10,logb2);
b3=power(10,logb3);
Ks1=power(10,logKs1);
cOH=power(10,pH-14);
%this is the final curve equation
logcNP=Ks1./(cPO4.*cNH4).*(1+b1.*cH2PO4+b0.*cHPO4+b2.*cPO4+b3*cOH);
%I am trying to calculate a crosssection of the curve above with y value of 0.005 (cMg at the start of the code).
crossectionY=cMg
[crossectionY, indexAtCrossectionY] = min(abs(logcNP-crossectionY));
crossectionX = pH(indexAtCrossectionY(1))
xcidx = find(diff(sign(logcNP - 0.005))); % ADDED: Approximate Crossing In dices For 'logcNP' At 0.005
for k = 1:numel(xcidx)
idxrng = max(1,xcidx(k)-1) : min(numel(pH),xcidx(k)+1); % ADDED: Index Range For Interpolation
pHx(k) = interp1(logcNP(idxrng),pH(idxrng), 0.005); % ADDED: Values Of 'pH' For 'logcNP = 0.005'
end
pHx
%This is plot of the final curve equation (black) and calculated crossection with magenta line
semilogy(pH,logcNP,'-', 'color', 'k', 'MarkerFaceColor','w', 'MarkerEdgeColor', 'm');
box off;
hold on;
plot([0 crossectionX],[cMg cMg],'-', 'color', 'm', 'LineWidth',0.5);
plot(pHx, ones(size(pHx))*0.005, 'sc', 'MarkerSize', 10) % ADDED: Plot Cyan Squares At Intersection Points
box off;
hold on;
axis([0 14 0.000000000000000001 1000000000000000000]);
xlabel('pH (-)');
ylabel('c {(mol/dm^{3})}');
axis square;
grid on;
axis on;
box off;
hold off;
.
2 Comments
Star Strider
on 15 May 2023
My pleasure!
If my Answer helped you solve your problem, please Accept it!
.
More Answers (0)
See Also
Categories
Find more on Lighting, Transparency, and Shading 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!