matlab死循环怎么调试,一直是死循环，不知道原因

% 程序说明：适用于矩形柱，截面尺寸600mm*550mm，

clc;clear all;close all;

%-------------------尺寸属性---------------------------------------------

h=600;          %柱截面的高度mm

d=550;          %柱截面的宽度mm

H=1200;         %柱的高度

N_0=3000000;    %为轴力

c=30;           %为保护层厚度

D=24;           %纵筋直径

d_1=10;         %箍筋直径

m=4;            %箍筋肢数

nn=2;            %箍筋肢数

s_01=75;        %箍筋间距

p=4;            %高一侧纵筋根数

q=3;            %宽一侧纵筋根数

%%%---------------材料参数---------------------

Es=2e5;             %钢筋的弹性模量Es(N/mm2)

fy=375;             %纵筋屈服强度

fu=635.6;           %纵筋峰值强度

ey=fy/Es;           %纵筋屈服应变

esm=0.0012;          %箍筋最大应变 (usually ~0.10-0.15)*    %%%有问题

e1=0.002;       %无约束混凝土峰值应力对应的应变

fco=21.4;       %无约束混凝土抗压强度

ecu=0.0738;             %约束混凝土极限压应变

% ====================              结束输入数据       ==========================

%---------------------------------------主程序-----------------------------------

%将截面分成三部分，中间为约束混凝土和普通混凝土，其余为普通混凝土本构，分层就为柱高。

Ms=zeros(1,h);                   %初始化弯矩

Ns=zeros(1,h);                   %初始化轴力

eall=zeros(1,h);                 %调试时加的

angle=linspace(0,0.0001,2000);  %设0.0001为最小曲率

e_0=0;                          %最初应变

angle_02=0;

for l=1:2000

angle_01=angle(l);

n=1;

wucha=0.006;                %0.006为规定的误差值

e_0=e_0;

e_1=0.000001;            %预定截面中心应变调整增量

while (wucha>0.005)&&(n<300)   %迭代次数过多或误差超过容限则终止

for u=1:h

eall(u)=e_0+angle_01*(h/2-u);   %总应变

end

for s=1:(c+d_1)   %保护层混凝土

Ns(s)=unconfined(eall(s),fco,e1)*h;

Ms(s)=Ns(s)*(h/2-s);

end

for s=(c+d_1+1):1:(h-c-d_1)    %约束混凝土和保护层%%%%这里有问题。

Ns(s)=unconfined(eall(s),fco,e1)*(2*c+2*d_1)*(h-2*c-2*d_1)*2+confined(eall(c),fco,nn,d,h,c,d_1,fy,m,D,p,q,s_01)*(h-2*2*c-2*d_1)*(d-2*c-2*d_1);

Ms(s)=Ns(s)*(h/2-s);

end

for s=(h-c-d_1+1):1:h     %保护层混凝土

Ns(s)=unconfined(eall(s),fco,e1)*h;

Ms(s)=Ns(s)*(h/2-s);

end

ss=c+d_1;    %钢筋

while round(ss)<=(h-c-d_1+1)    %大于h-c-d_1+1就跳出循环

f(round(ss))=steel(eall(round(ss)),fu,Es,ey,fy);

if ss==(c+d_1)

area=q*d^2*pi/4;

elseif  (round(ss)<=(h-c-d_1+1))&&(round(ss)>=(h-c-d_1-1))    %四舍五入是因为除法可能会出现小数

area=q*d^2*pi/4;

else

area=2*d^2*pi/4;      %多数试验柱仅有边柱，故仅考虑边侧钢筋

end

Ns(round(ss))=f(round(ss))*area;

Ms(round(ss))=f(round(ss))*area*(h/2-ss);

ss=ss+(h-2*c-2*d_1)/(p-1);

end

N=sum(Ns);

M=sum(Ms);

wucha=abs((N_0-N)/N_0);

if wucha>0.001&&(N>=N_0)

e_0=e_0-e_1;

elseif wucha>0.001&&(N

e_0=e_0+e_1;

end

n=n+1;

end

MM(l)=M;

if eall(c+d_1)>ecu

angle=angle(:,1:(l));

break         %分析终止

end

%  if abs(eall(h-c-d_1))>=0.98*ey&&abs(eall(h-c-d_1))<=1.02*ey

if abs(eall(h-c-d_1+1))>=ey&&abs(eall(h-c-d_1-1))<=ey

angle_02=angle(l);      %截面最外侧纵向受拉钢筋首次屈服时的曲率

end

end

lp=0.08*H+0.022*fy*D;

MM=MM/10^6;

ff=zeros(1,l);                   %初始化位移

PP=zeros(1,l);                   %初始化荷载

for r=1:l

if angle(l)<=angle_02

ff(r)=angle(l)*H^2/3 ;    %ff为柱顶弯曲变形

elseif angle(l)>angle_02

ff(r)=angle_01*H^2/3+(angle(l)-angle_01)*lp*(H-lp/2);

end

PP(r)=(MM(r)-N_0*ff(r)/10^6)/H/1000;

end

ff;

PP;

figure; plot(ff(r),PP(r),'ro'); grid on;

xlabel('Displacement(mm)','FontSize',16); ylabel('load(KN)','FontSize',16); grid on

title('荷载位移曲线','FontSize',16);