# Marking a specific point in a graph made out of a loop

2 views (last 30 days)
Shimon Katzman on 4 Dec 2019
Commented: Star Strider on 5 Dec 2019
Hi everyone.
Trying to mark 2 points at the 4 curves in the graph:
1) Marking the peak in each graph (the Mmax point)
2) Marking the point where the value of M belongs to epss=epssh in each graph.
Thank You very much.
b=400; %mm
d=500; %mm
Asv = [3000 5000 3000 5000]; %mm^2
fckv = [30 30 90 90]; %Mpa
for k = 1:numel(fckv)
fck = fckv(k);
Ecshah=57000/145*(fck*145)^0.5; %Mpa
Es=200000; %Mpa
Esh=8500; %Mpa
fy=500; %Mpa
fsu=750; %Mpa
epssh=0.009;
epssu=0.075;
eps0=1.027*10^-7*fck*145+0.00195;
kshah=0.025*fck*10^3;
A=Ecshah*eps0/fck;
P=Esh*((epssu-epssh)/(fsu-fy));
epsy=fy/Es;
epscmv = linspace(0.1, 100, 5000)*1E-3;
As = Asv(k);
for i=1:numel(epscmv);
epscm = epscmv(i);
epss=@(c) (d-c)/c*epscm;
funCshah=@(epsc) (1-(1-epsc./eps0).^A) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15) .* (epsc>eps0);
compression=@(c) b*fck*c/epscm*integral(funCshah,0,epscm)/1000;
sigmaSteel=@(c) Es*epss(c) .* (epss(c)<=epsy) + fy .* (epss(c)>epsy & epss(c)<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss(c))./(epssu-epssh)).^(1/P)) .* (epss(c)>epssh & epss(c)<=epssu) + 0 .* (epss(c)>epssu);
tension=@(c) sigmaSteel(c).*As/1000;
c(i)=fsolve(@(c) compression(c)-tension(c),1000);
funM=@(epsc) (1-(1-epsc./eps0).^A).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc>eps0);
M(i,k)=b*fck*c(i)/epscm*integral(funM,0,epscm)/1000000;
phi(i,k)=epscm/c(i);
end
end
[Mmax,idx]=max(M) %[kNm]
phiAtMmax=phi(idx) %[1/mm]
epsAtMmax=epscmv(idx)*1000 %[promil]
figure
hold all
for k = 1:numel(fckv)
plot(phi(1:idx+50,k), M(1:idx+50,k))
lgdstr{k} = sprintf('fck = %2d [Mpa], As = %4d [mm^2]',fckv(k), Asv(k));
plot(phi(idx,k),M(idx,k),'r*') % marking the Mmax
end
hold off
grid on
xlabel('phi [1/mm]')
ylabel('Moment [kNm]')
legend(lgdstr)

Shimon Katzman on 4 Dec 2019
Star Strider on 4 Dec 2019
Working on it. Windows 10 gave me a BSOD a few minutes ago for the second time in as many days (I have no kind works for Micro\$oft just now), and it took a while to put the pieces of everything back together.

Star Strider on 4 Dec 2019
Try this:
b=400; %mm
d=500; %mm
Asv = [3000 5000 3000 5000]; %mm^2
fckv = [30 30 90 90]; %Mpa
tic
for k = 1:numel(fckv)
fck = fckv(k);
Ecshah=57000/145*(fck*145)^0.5; %Mpa
Es=200000; %Mpa
Esh=8500; %Mpa
fy=500; %Mpa
fsu=750; %Mpa
epssh=0.009;
epssu=0.075;
eps0=1.027*10^-7*fck*145+0.00195;
kshah=0.025*fck*10^3;
A=Ecshah*eps0/fck;
P=Esh*((epssu-epssh)/(fsu-fy));
epsy=fy/Es;
epscmv = linspace(0.1, 100, 5000)*1E-3;
As = Asv(k);
for i=1:numel(epscmv);
epscm = epscmv(i);
epss=@(c) (d-c)/c*epscm;
funCshah=@(epsc) (1-(1-epsc./eps0).^A) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15) .* (epsc>eps0);
compression=@(c) b*fck*c/epscm*integral(funCshah,0,epscm)/1000;
sigmaSteel=@(c) Es*epss(c) .* (epss(c)<=epsy) + fy .* (epss(c)>epsy & epss(c)<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss(c))./(epssu-epssh)).^(1/P)) .* (epss(c)>epssh & epss(c)<=epssu) + 0 .* (epss(c)>epssu);
tension=@(c) sigmaSteel(c).*As/1000;
c(i)=fsolve(@(c) compression(c)-tension(c),1000);
funM=@(epsc) (1-(1-epsc./eps0).^A).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc<=eps0) + exp(-kshah*(epsc-eps0).^1.15).*(d-c(i)+(c(i)./epscm).*epsc) .* (epsc>eps0);
M(i,k)=b*fck*c(i)/epscm*integral(funM,0,epscm)/1000000;
phi(i,k)=epscm/c(i);
end
idx(k) = find(diff(M(:,k)) < 0, 1, 'first')%[kNm]
Mmax(k) = M(idx(k),k) %[kNm]
phiAtMmax(k)=phi(idx(k)) %[1/mm]
epsAtMmax(k)=epscmv(idx(k))*1000 %[promil]
end
toc
% [Mmax,idx]=max(M) %[kNm]
% phiAtMmax=phi(idx) %[1/mm]
% epsAtMmax=epscmv(idx)*1000 %[promil]
epss_c = fsolve(@(c) epss(c)-epssh, 100)
c_idx = find(c >= epss_c, 1, 'first')
epscm_exact = interp1(c([-5 5]+c_idx), epscmv([-5 5]+c_idx), epss_c)
figure
hold all
for k = 1:numel(fckv)
hp(k) = plot(phi(1:idx+50,k), M(1:idx+50,k));
% plot([1 1]*epss_c, ylim, ':k')
lgdstr{k} = sprintf('fck = %2d [Mpa], As = %4d [mm^2]',fckv(k), Asv(k));
plot(phi(idx(k),k),M(idx(k),k),'r*') % marking the Mmax
end
hold off
grid on
xlabel('phi [1/mm]')
ylabel('Moment [kNm]')
hl = legend(hp,lgdstr);
hl.FontSize = 7;
% text(epss_c, 400, sprintf('\\leftarrow epss_c = %.1f', epss_c), 'HorizontalAlignment','left')
This correctly plots the maxima of every curve (although the green curve maximum is beyond the green curve plotted limit (30 Mpa, 3000 mm^3), and I corrected the legend problem as well, so that it only refers to the lines, not the markers.
I do not understand about ‘epss_c’ with respect to every curve. It appears to be beyond every curve. In any event, I have no idea how to calculate it with respect to plotting it as a function of ‘phi’.
The Windows 10 crash (that I did not immediately recognise) put me abouit two hours behind.

Star Strider on 5 Dec 2019
As always, my pleasure!
In your code, ‘c’ is a vector, not a matrix. To create a matrrix for it, the easiest way would be to add:
c_mtx(i,k) = c(i);
just before the end of the ‘i’ loop, then change the plot call to:
hp(k) = plot(epscmv(1:idx(k)), c_mtx(1:idx(k),k));
EDIT (5 Dec 2019 at 14:03)
That appears to plot correctly. You will need to change the related calls to show the correct xlabel and ylabel strings, and remove or change the second plot call.
Shimon Katzman on 5 Dec 2019
Thank you very much!!!!
Star Strider on 5 Dec 2019
My pleasure!

Walter Roberson on 5 Dec 2019
Repaired code is attached. Some optimizations were made as well.

#### 1 Comment

Shimon Katzman on 5 Dec 2019
Wow Thank you very much.
Thank you.