Someone

clear;
clc;
format long;
disp('     L                ___');
disp(' -----------           |') ;
disp(' |          |  A       |');
disp(' ----    ---           |   '); 
disp('    |   |              |  D');
disp('    |   |              |');
disp('    |   |              |');
disp('    |   |              |');
disp('    |___|             _|_');
disp('      B                ');
disp(' ');


fcprima=200;%=input('INGRESE "fc prima" en megapascales: ');
ES=2100000;
EC=15100*sqrt(fcprima);
FY=2800;
n=ES/EC;
fr=2*sqrt(fcprima);
% Basta con eliminar el valor y el % para poder ingresar datos. Estan
% puestos los datos del problema, con 3 barras de diametro 26mm.

L=50;%input('INGRESE L (largo del ala): ');
A=14;%input('INGRESE A (espesor del ala): ');
B=25;%input('INGRESE B (espesor del alma): ');
D=100;%input('INGRESE D (alto de la viga): ');
b=3;%input('INGRESE CUANTAS BARRAS DE ACERO: ');
d=26;%input('INGRESE DIAMETRO DE LAS BARRAS EN MM: ');
r=97.5;%input('INGRESE DISTANCIA DESDE TOPE DEL ALA HASTA EL CENTRO DE LAS BARRAS: ');
as=b*(pi*d^2/400); % en cm^2

% SECCION TRANSFORMADA
cm=[L*A*A/2+B*(D-A)*(A+(D-A)/2)+as*(n-1)*(r)]/[L*A+B*(D-A)+(n-1)*as]; % Centro de area de seccion transf.
it=L*A^3/12+L*A*(cm-A/2)^2+(B*(D-A)^3)/12+B*(D-A)*((D-A)/2-cm)^2+(n-1)*as*(r-cm)^2; % Inercia seccion transf.


% Primera grieta, hormigon comportandose elasticamente.
fcmax=2*sqrt(fcprima);
mcr=fcmax*it/(cm);
curva1=mcr/(EC*it);

% Cambia la inercia
c1=1;
tol=10000;
while tol>1000
    C=(L*A*(c1-A/2)+B*(c1-A)^2/2);
    T=n*as*(r-c1);
    tol=abs(C-T);
    c1=c1+0.01;
end

i1=L*A^3/12+L*A*(c1-A/2)^2+B*(c1-A)^3/12+B*(c1-A)*(c1-A)/2;
curva2=mcr/(EC*i1);

% Acero fluyendo y EC=E0
ecmax=0.002;
Ty=as*FY; %
c3=0.1;
tol=10000;
while tol>1000
    C3=(2/3)*fcprima*B*c3+2*(2/3)*fcprima*((L-B)/2)*c3; % Se divide en tres partes la seccion, dos son iguales.
    tol=abs(C3-Ty);
    c3=c3+0.01;
end
m3=Ty*((5/8)*c3+r-c3);
curva3=ecmax/c3;


% Deformacion última del hormigón.
ecmax=0.003;
Ty=as*FY;
tol=10000;
c4=0.1;
while tol>1000
    C41=fcprima*B*c4/3+2*fcprima*((L-B)/2)*c4/3;
    C42=(2/3)*fcprima*B*(2/3)*c4+2*(2/3)*fcprima*((L-B)/2)*(2/3)*c4;
    C4=C42+C41;    
    tol=abs(C4-Ty);
    c4=c4+0.01;
end
m4=C41*(r-c4+5*c4/6)+C42*(r-c4+(10/24)*c4);
curva4=ecmax/c4;


% Gráfico
momentos = [0 mcr mcr m3 m4];
curvaturas = [0 curva1 curva2 curva3 curva4 ];
plot(curvaturas,momentos)