Matlab - Ayuda con Algoritmo

% PROGRAMA PARA EL CALCULO Y DISEÑO DE COLUMNAS RECTANGULARES Y CUADRADAS DE CONCRETO ARMADO

disp('CÓDIGO PARA DIBUJAR EL DIAGRAMA DE INTERACCIÓN DE COLUMNAS RECTANGULARES Y CUADRADAS')

disp('UNIVERSIDAD EXPERIMENTAL DE LOS LLANOS OCCIDENTALES EZEQUIEL ZAMORA')

disp('INGRESO DE DATOS')

fc=input('Resistencia del Concreto Kg/cm2:');

if fc<=280

B1=0.85; %Valor del factor B1

else B1=(1.05-fc)/1400; %Valor del factor B1

end

rec=input('Recubrimieto cm:');

Fy=input('Resistencia del Acero Kg/cm2:');

Es=input ('Módulo de elasticidad del acero Es:');

B=input('Base de Columna cm:');

if B>=30;

display('Base minima de 30 cm');

end

H=input('Altura de Columna cm:');

if B>=30;

disp('Base minima de 30 cm');

end

Nlas=input('Numero de Lineas de Acero:');

Ast=0;

for i=1:1:Nlas

 %Diametros de Barras Segun Lineas de Acero

 Db(i)=input('Ingrese el Diametro de la Barra por Fila en Pulg:');

 Nb(i)=input('Ingrese el Numero de Barras por Fila:');

 d(i)=input('Ingrese la Separación entre Barras desde la Fibra Comprimida:');

 As(i)=Nb(i)*((Db(i)*2.54)^(2))*pi/4;

 Ast=sum(As);

end

% CALCULO DE LOS MOMENTOS Y CARGAS NOMINALES

disp('CÁLCULOS')

disp('1) Fuerza a compresion pura en Ton');

% PARA EL PUNTO P0

P0=(0.75*0.65*(0.85*fc*((B*H)-Ast)+Fy*Ast))/1000

M0=0;

eb=M0/P0;

dt=H-rec;

Ecu=0.003;

C=0;

fs=0;

Es1=0;

Est=0;

fact=0;

% Condicionante para ver que aceros están en tracción y compresión

for i=1:1:H;

  C(i)=i/B1;

  display('2)Esfuerzos para cada Linea de Acero');

  % Aceros a Tensión Negativos , Aceros a Compresión Positivos

  for j=1:Nlas

  Es1(j)=Ecu*((C(i)-d(j))/C(i))

  fs(j)=Es*Es1(j)

  if abs(fs(j))>Fy

  if fs(j)>=0;

  fs(j)=Fy

  else

  fs(j)=-Fy

  end

  end

  Fi(j)=fs(j)*As(j)

  display('3)Calculo del factor de Minoracion de Carga')

  Est(j)=Ecu*((C(i)-dt)/C(i))

  if abs(Est(j))<=0.002

  fact(j)=0.65;

  else fact(j)=0.90;

  end

  display('4)Calculo de Cargas y Momentos')

  P(i)=((0.85*fc*i*B+sum(Fi))/1000)*fact(j)

  M(i)=((((0.85*fc*i*B)/1000)*((H/2-i/2)/100))+(Fi(j)/1000)*((d(j)-H/2)/100))*fact(j)

  end

end