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')