clear all

% Kraftkoeffizienten

ktc=-3956.878095;
krc=-1982.30869;
kac=-1231.365042;
kte=-54.39686283;
kre=-42.83233385;
kae=3.226134429;

n_Umdrehungen=1; % Anzahl Werkzeugumdrehungen

%%%% Systemgrößen Werkzeug %%%%

% Durchmesser [mm]
D=6;
R=D/2;
%Drallwinkel [°]
io=30;
io=io/180*pi; %rad
%Zähnezahl
Ztotal=2;

%%%% Stellgrößen %%%%

% Drehzahl in [U/min] für Stahl
N=3820;
%Schnitttiefe [mm]
ap=0.15;
% Zeilenbreite [mm]
ae=0.3;
%Zahnvorschub [mm]
fz=0.062;
% Anstellwinkel [°]
beta_fn=83;
% Teilungswinkel [°]
phi_p=2*pi/Ztotal;


%%%% Integrationsgrenzen %%%%

% Minimaler und maximaler axialer Eingriffswinkel [°]
kappa_min=asin(1-(ap/R))-(90-beta_fn)/180*pi;
kappa_max=beta_fn/180*pi+asin(ae/D);

%% Aufteilung der Z-Achse (im Eingriff)
z1=R*(1-cos(kappa_min));
z2=R*(1-cos(kappa_max));
% Anzahl Intervalle
L=100;
% Länge Intervalle
dz=(z2-z1)/(L-1);

% Radiale Eingriffswinkel (Ein- und Austritt)
% phi_aus=pi/2+asin(fz/(2*R))+pi/2;
% phi_ein=asin(1-(ap/R))+pi/2;

z=(z1:dz:z2)'; % Erzeuge Vektor z und befülle ihn
E=(R-z)./R; % Hilfsgröße
kappa=asin(sqrt(1-(E).^2)); % Berechne axialen Eingriffswinkel
R_z=R*(1-E.^2).^0.5; % Radius entlang der z-Achse
i_z=atan(R_z/R*tan(io)); % Drallwinkel entlang der z-Achse

% Radiale Eingriffswinkel (Ein- und Austritt)
phi_aus=pi/2+asin(fz./(2*R_z))+pi/2;
phi_ein=asin(1-(ap./R_z))+pi/2;

% Länge eines Schneidensegments
dS=dz*((tan(io))^2*(1-E.^2)+1./(1-E.^2)).^0.5;

%Spanungsbreite
db=dz./sin(kappa);

M=360; % Anzahl der Intervalle (Drehwinkel)
phi_end=2*pi*n_Umdrehungen; % Maximaler Drehwinkel
d_phi0=phi_end/(M-1); % Länge des Drehwinkelintervalls
phi0=(0:d_phi0:phi_end); % Erzeuge Vektor phi0 und befülle ihn
n=length(phi0); % Länge des Vektors phi0

%Testweise
h=zeros(n,L);
dFt=zeros(n,L);
dFr=zeros(n,L);
dFa=zeros(n,L);


for  i=1:n % Schleife für den Drehwinkel des Werkzeugs
     % Erzeuge Kraftkomponentenvektoren der Länge n und befülle sie mit 0
     % außerhalb des Eingriffs
     Fx(i)=0;
     Fy(i)=0;
     Fz(i)=0;     
     for k=1:Ztotal   % Summe über die Anzahl der Schneiden
            for j=1:L % Integration über z
                % Eingriffswinkel
                phi=phi0(i)+(k-1)*2*pi/Ztotal-z(j)*tan(io)/R;
 
            
                    if (mod(phi,2*pi)<=phi_aus(j))&(mod(phi,2*pi)>=phi_ein(j)) %Schneide im Eingriff?
                        % Berechnung der infinitesimalen Kräfte
%                         DFt=kte*dS(j)+ktc*fz*sin(phi)*dz;
%                         DFr=kre*dS(j)+krc*fz*sin(phi)*dz;
%                         DFa=kae*dS(j)+kac*fz*sin(phi)*dz;

                        %Spanungsdicke
                        h(i,j)=fz*sin(phi)*sin(kappa(j));
                        
                        DFt=kte*dS(j)+ktc*h(i,j)*db(j);
                        DFr=kre*dS(j)+krc*h(i,j)*db(j);
                        DFa=kae*dS(j)+kac*h(i,j)*db(j);
                        
               
                      
                        % Transformationsmatrix 
                       
                        T=[-cos(phi) -sin(kappa(j))*sin(phi) -cos(kappa(j))*sin(phi);...
                        sin(phi) -sin(kappa(j))*cos(phi) -cos(kappa(j))*cos(phi);...
                         0            cos(kappa(j))          -sin(kappa(j))];
                    
                 
                     dFt(i,j)=DFt;
                     dFr(i,j)=DFr;
                     dFa(i,j)=DFa;
                     
                     
                     DF=T*[DFt;DFr;DFa];
                         
                        DFx=DF(1);
                        DFy=DF(2);
                        DFz=DF(3);
                        
                        % Summation
                        Fx(i)=DFx+Fx(i);
                        Fy(i)=DFy+Fy(i);
                        Fz(i)=DFz+Fz(i);
                    end     
            end
     end
end



% Rücktransformation von TCS in MCS
F_TCS=[Fx; Fy; Fz];
for i=1:length(Fx)
    F_MCS(:,i)=TCStoMCS(F_TCS(:,i));
    
end

figure(2);plot(phi0*180/pi,F_MCS(1,:),'r',phi0*180/pi,F_MCS(2,:),'g',phi0*180/pi,F_MCS(3,:),'b',phi0*180/pi,F_R,'k');grid on;