MATLAB: Finding a value and marking it in a graph

helpmarkingvalues

Hi everyone,
I am trying to find the value of M when epss=epssh=0.009 for every graph (there are 4 of them) and to mark it on the graph unsucssesfully.
Can anyone help me please with this please?
b=400;  %mm

d=500;  %mm
Asv = [3000 5000 3000 5000]; %mm^2
fckv = [30  30  90  90]; %Mpa





num_fckv = numel(fckv);
epscmv = linspace(0.1, 100, 5000)*1E-3;
num_epscvm = numel(epscmv);
fsolve_opts = optimoptions('fsolve','Display','off');
M = zeros(num_epscvm, num_fckv);
phi = zeros(num_epscvm, num_fckv);
idx = zeros(1, num_fckv);
Mmax = zeros(1, num_fckv);
phiAtMmax = zeros(1, num_fckv);
epsAtMmax = zeros(1, num_fckv);
tic
for k = 1:num_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;
    As = Asv(k);
    c0 = 1000;
    c = zeros(1, num_epscvm);
    for i=1:num_epscvm
        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), c0, fsolve_opts);
        c0 = c(i);
        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);
        epss1(k,i)= (d-c(i))/c(i)*epscm;
        sigmaSteel1(k,i)= Es*epss1(k,i) .* (epss1(k,i)<=epsy) + fy .* (epss1(k,i)>epsy & epss1(k,i)<=epssh) + (fsu+(fy-fsu)*abs((epssu-epss1(k,i))./(epssu-epssh)).^(1/P)) .* (epss1(k,i)>epssh & epss1(k,i)<=epssu) + 0 .* (epss1(k,i)>epssu);
        c_mtx(i,k) = c(i);
    end
    idx(k) = find(diff(M(:,k)) < 0, 1, 'first')
    [Mmax(k),idx(k)]=max(M(:,k)) %[kNm]
    phiAtMmax(k)=phi(idx(k),k)  %[1/mm]
    epsAtMmax(k)=epscmv(idx(k))*1000 %[promil]
end
toc
% epss_c = fsolve(@(c) epss(c)-epssh, 100, fsolve_opts)
% 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
hp = zeros(1, num_fckv);
lgdstr = cell(1, num_fckv);
for k = 1:num_fckv
    hp(k) = plot(phi(1:idx(k)+50,k), M(1:idx(k)+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')

Best Answer

The only way I can see to calculate ‘M’ where ‘epss’ is equal to 0.009 is to duplicate the first nested loop, defining: 
epss = @(c) epssh;
If you do that: 
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) epssh;
%         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);
        M009(i,k)=b*fck*c(i)/epscm*integral(funM,0,epscm)/1000000;
        phi(i,k)=epscm/c(i);
        c_mtx(i,k) = c(i);
    end
end
toc
figure
mesh(M009)
grid on
and also do the surface plot, you can see the result.
I have no idea what to do with this result, so I leave that to you.
Related Question