Stuff

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);
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
fcmax=2*sqrt(fcprima);
mcr=fcmax*it/(cm);
curva=mcr/(EC*it);
mcr;
curva;

% TERMINA LA LINEALIDAD DEL HORMIGON


%si el centro de area de compresion esta en el ala

%c1=1;
%tol=100;
%while tol>10
%    C=L*c1^2/2;
%    T=n*as*(r-c1);
%    tol=abs(C-T);
%    c1=c1+0.2;
%end
%c1
 
%si el centro de area de compresion esta en el alma

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.0005;
end

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


ecm=0.4*fcprima/EC;
curv2=ecm/c1;

momentos = [0 mcr mcr ];
curvaturas = [0 curva curva1];
plot(curvaturas,momentos)