Matlab分段积分问题函数调用出现运行错误

[勇敢的搬砖人]

R=5;        Nt=4;    alpha=pi/6;    % Parameters about the mill 
ar=7.5;     ap=5;    ft=0.1;        % Cutting parameters
kts=500;    krs=500;    kas=500;    % Sheer cutting force coefficient in tangential, radial, axial direction, respectively
ktp=500;    krp=500;    kap=500;    % Plough cutting force coefficient in tangential, radial, axial direction, respectively
delta=2*pi/Nt;                      % The tooth space angle
beta0=ap*tan(alpha)/R;
phi0=acos((R-ar)/R);                % Milling angle, which equals phi_ex minus phi_st
phi_st=pi-phi0;   phi_ex=pi;        % Start and exit cutting angle
 syms t;                     

theta=0:0.01:30*pi;                 % Set the step length and the number of cycle periods


% Define the cutting force equations
gx=R/tan(alpha)*(-ft*kts*sin(t)^2/2-ft*krs*(-sin(t)*cos(t)/2+t/2)-ktp*sin(t)+krp*cos(t));
gy=R/tan(alpha)*(ft*kts*(-sin(t)*cos(t)/2+t/2)-ft*krs*sin(t)^2/2-ktp*cos(t)-krp*sin(t));
gz=R/tan(alpha)*(-ft*kas*cos(t)+kap*t);

% Initialization
Fx1=zeros(size(theta));Fx2=zeros(size(theta));Fx3=zeros(size(theta)); 
Fy1=zeros(size(theta));Fy2=zeros(size(theta));Fy3=zeros(size(theta)); 
Fz1=zeros(size(theta));Fz2=zeros(size(theta));Fz3=zeros(size(theta)); 
F1=zeros(size(theta));F2=zeros(size(theta));F3=zeros(size(theta)); 
hm1=zeros(size(theta));hm2=zeros(size(theta));hm3=zeros(size(theta));
b=zeros(size(theta));a=zeros(size(theta));
N=length(theta); 
for k=1:N
    if mod(theta(k),2*pi)>= phi_st && mod(theta(k),2*pi)< phi_st + beta0
        b(k)=mod(theta(k),2*pi);
        a(k)=phi_st;
        Fx1(k)=gx(b(k))-gx(a(k));
        Fy1(k)=gy(b(k))-gy(a(k));
        Fz1(k)=gz(b(k))-gz(a(k));        
        hm1(k)=(b(k)-a(k))*R/tan(alpha);
        
    elseif mod(theta(k),2*pi)>= phi_st + beta0 && mod(theta(k),2*pi)< phi_ex
        b(k)=mod(theta(k),2*pi);
        a(k)=mod(theta(k),2*pi)-beta0;
        Fx1(k)=gx(b(k))-gx(a(k));
        Fy1(k)=gy(b(k))-gy(a(k));
        Fz1(k)=gz(b(k))-gz(a(k));            
        hm2(k)=ap;
        
    elseif mod(theta(k),2*pi)>= phi_ex && mod(theta(k),2*pi)<= phi_ex + beta0
        b(k)=phi_ex;
        a(k)=mod(theta(k),2*pi)-beta0;
        Fx1(k)=gx(b(k))-gx(a(k));
        Fy1(k)=gy(b(k))-gy(a(k));
        Fz1(k)=gz(b(k))-gz(a(k));        
        hm1(k)=(b(k)-a(k))*R/tan(alpha);
    end
end
Fx=Fx1+Fx2+Fx3;
Fy=Fy1+Fy2+Fy3;
Fz=Fz1+Fz2+Fz3;
F=F1+F2+F3;
hm=hm1+hm2+hm3;

subplot(211)
plot(theta,Fx,'r')
hold on
plot(theta,Fy,'g')
hold on
plot(theta,Fz,'b')
hold on
plot(theta,F,'k')
grid on
legend('Fx','Fy','Fz','\SigmaF')
% title(Cutting Force)
subplot(212)
plot(theta,hm)
grid on
legend('Cutting Height')