%%  Boyd, Kreinin, Ting, and Webster 2019

%%  GEAR PROPERTIES
%%  THROUGH HARDENED STEEL - GRADE 2 - SPUR GEARS
%   A,C = PINIONS;      B,D = GEARS;

Functions; %Include functions file
clear


%% Givens
%%Gear bending stress variables
reqd_speed = 6 / 3.6 * 39.3701; %kph->m/s->in/s mower driving speed
life_hours = 1000; %hours
min_gear_wheel_rad_clearance = 0.05 * 39.3701; % 5 cm off the ground ->inches
eng_net_T = 2.0; %[lb.ft] at 4500rpm
eng_rpm = 4500; %Motor RPM
eng_w = eng_rpm * (2*pi/60); %Motor [rad/s]
Ce = 1; %Eqn 14-35 (for all other conditions)
Cf = 1; %surface condition factor
Ch = 1; %hardness ratio factors for pitting resistance, Section 14-12, gear only
Cp = 2300; %elastic coefficient, sqrt(psi), Table 14-8
Cmc = 1; %uncrowned teeth
Cpm = 1; % For calculating Km
H = 1.7; %horsepower (hp) at 4500rpm
HB = 400; %Brinell Hardness of gear material
Kb = 1; %rim-thickness factor
Ko = 1; %overload factor
Kr = 1;  %reliabilty factor
Kt = 1; %Temperature Factor (1 if T< 250 F)
Kv_a = 1.416; %dynamic factor for gears A and B
Kv_b = 1.1404; %dynamic factor for gears C and D
Sh = 1.5; %Wear factor of Safety
Sf = 2; %Bending safety factor
Sc = 349*HB + 34300;  %[lbf/inch^2], allowable contact stress for grade 2
St = 102*HB + 16400; %allowable bending stress for grade 2 through hardened steel, Fig. 14-2


%% Angular Speed of Req'd and Gearing Ratio
wheel_rad = 0.2 /2*39.3701; %Wheel radius m->in
w_o = reqd_speed / wheel_rad; %Wheel angular velocity in rad/s
w_o_rpm = w_o * 60 / (2*pi); %Wheel angular velocity in rpm
g_ratio = eng_w / w_o;% Total gear ratio to step down from motor to wheel
optimal_ratio = sqrt(g_ratio);%Optimal gear ratio
redxn_weight1 = 0.666;%Gear reduction weight for gears C and D
redxn_weight2 = 1.2; %Gear reduction weight for gears A and B
Gear(4).r = wheel_rad - min_gear_wheel_rad_clearance; %Gear D radius
Gear(3).r = Gear(4).r / (optimal_ratio*redxn_weight1);%Gear C radius
g_ratio1 = Gear(4).r/Gear(3).r; %Gear ratio between gears C and D
g_ratio2 = g_ratio / g_ratio1; %Gear ratio between gears A and B
Gear(1).r = Gear(3).r*redxn_weight2; %gear A > C to compensate for faster speed
Gear(2).r = Gear(1).r * g_ratio2; %similar ground clearance to gear D
valid = (((wheel_rad + Gear(4).r + Gear(3).r) - Gear(2).r) >= (5/2.54));
%Check that there is enough clearance with the gear radii

%% Gear RPMs
Gear(1).rpm = eng_rpm;
Gear(2).rpm = Gear(1).rpm / g_ratio2;
Gear(3).rpm = Gear(2).rpm;
Gear(4).rpm = w_o_rpm;

%% Gear life cycles
N = "Cycles";
for i = 1:4
    Gear(i).N = life_hours * 60 * Gear(i).rpm; %Cycles- Hours * 60 mins/hr * rpm
endfor


Gear(1).Pd = 10; %teeth/in, diametral pitch of A
Gear(2).Pd = Gear(1).Pd;
Gear(3).Pd = 6;
Gear(4).Pd = Gear(3).Pd;

for i = 1:4
    Gear(i).mg = 1/Gear(i).Pd; %[in/tooth]
    Gear(i).a = Gear(i).mg; % [in/tooth], gear addendum
    Gear(i).d = Gear(i).r * 2; % [in], pinion pitch diameter
endfor

%Number of teeth
Gear(1).Np = round(Gear(1).Pd * Gear(1).d);
Gear(2).Np = round(Gear(1).Np * g_ratio2);
Gear(4).Np = round(Gear(4).Pd * Gear(4).d);
Gear(3).Np = round(Gear(4).Np / g_ratio1);

for i = 1:4
    Gear(i).Y = GetY(Gear(i).Np); %lewis form factor, Table 14-2
endfor

%Surface strength geometry factor, Project Report Figure 3 & 4
Gear(1).I = 0.1063;
Gear(2).I = 0.1063;
Gear(3).I = 0.0909;
Gear(4).I = 0.0909;

%[in], face width
Gear(1).F = 1;
Gear(2).F = 1;
Gear(3).F = 2.5;
Gear(4).F = 2;

%Calculate Cma and Cpf for each gear
for i = 1:4
    Gear(i).Cma = 0.127+0.0158*Gear(i).F-0.930*10^-4*Gear(i).F^2; %Table 14-9
    Gear(i).Cpf = GetCpf(Gear(i).F, Gear(i).d);
endfor

%spur-gear geometry factor, Fig. 14-6
Gear(1).J = 0.22;
Gear(2).J = 0.42;
Gear(3).J = 0.21;
Gear(4).J = 0.32;

for i = 1:4
    %load-distribution factor
    Gear(i).Km = 1+Cmc*(Gear(i).Cpf*Cpm+Gear(i).Cma*Ce);
    %size factor
    Gear(i).Ks = 1.192*(((Gear(i).F*sqrt(Gear(i).Y))/Gear(i).Pd).^0.0535);
    %stress-cycle factor for bending stress, Fig 14-14
    Gear(i).Yn = GetYn(Gear(i).N);
    %stress-cycle factor, Fig. 14-15
    Gear(i).Zn = GetZn(Gear(i).N);

    %%Gear bending stress equation, Eqn 14-15
    Gear(i).V = (pi*Gear(i).d*Gear(i).rpm) / 12;
    Gear(i).Wt = (33000*H)/Gear(i).V; %[lbf], tangential transmitted load
endfor


%%Stress Calculations%%
for i = 1:4
    if (i == 1 || i == 2)
        Kv = Kv_a;
    else
        Kv = Kv_b;
    endif
%%Gear bending stress equation, Eqn(14-15)
    Gear(i).sigma = GetSigma( Gear(i).Wt, Ko, Kv, Gear(i).Ks, Gear(i).Pd, Gear(i).F, Gear(i).Km, Kb, Gear(i).J);

%%Gear bending endurance strength equation, Eqn (14-17)
    Gear(i).sigma_all = GetSigmaAll(St, Sf, Gear(i).Yn, Kt, Kr);

%%Gear contact stress equation, Eqn 14-16
    Gear(i).sigma_c = GetSigmaC (Cp, Gear(i).Wt, Ko, Kv, Gear(i).Ks, Gear(i).Km, Gear(i).d, Gear(i).F, Cf, Gear(i).I);

%%Gear contact endurance strength, Eqn 14-18
    Gear(i).sigma_c_all = GetSigmaCAll ( Sc, Gear(i).Zn, Ch, Sh, Kt, Kr );
endfor

%Print all gears
for i = 1:4
    Gear(i)
endfor

%Print safety ratios
for i = 1:4
    bendsafety.i = Gear(i).sigma_all / Gear(i).sigma;
    wear_safety.i = Gear(i).sigma_c_all / Gear(i).sigma_c;

    printf ("Gear %d safety: %d \n wear safety: %d \n", i, bendsafety.i, wear_safety.i );
endfor

% Calculate final speed of wheel using tooth ratios
FinalSpeed = 4500 * (Gear(1).Np / Gear(2).Np) * ( Gear(3).Np / Gear(4).Np)
%Print error from the ideal 159.15 rpm
Error = (159.15 - FinalSpeed) * 100 / 159.1