function [Sigma] = SurfaceStressSpurGear(Qv, Vt, FW, DrivingShock, ...
    DrivenShock, size, I, Cp, surf, Wt, d)
% Calculates the surface stresses setup between two spur gears. English Units.
%
% Synopsis: [Sigma] = SurfaceStressSpurGear(Qv, Vt, FW, DrivingShock, ...
%    DrivenShock, size, I, Cp, surf, Wt, d)
%
% input:    Qv = (required) Gear quality. Constant between 3 and 11.
%           Vt = (required) pitch-line velocity. r cross omega.
%           FW = (required) Face width factor. Number in inches.
%           DrivingShock = (required) 'Uniform', 'Light Shock', or
%               'Medium Shock'.
%           DrivenShock = (required) 'Uniform', 'Moderate Shock', or
%               'Heavy Shock'.
%           size = (required) Must be a structured array. size must contain
%                   size.shape: 'rectangle', 'circle', or 'number'
%                   size.rotation: 'nonrotating', 'rotating', NA
%                   size.dim: 1x1 (radius of circle) or 1x2 (height (i.e. the
%                   cubed value) and width); must be in inches.
%                   fur number just put a number
%           I = (required) Surface geometry factor. See page 673 for values.
%           surf = (required) should equal 1. If surface is very rough then
%                   the value can be increased.
%           Wt = (required) The tangential force on the tooth
%           d = (required) diameter of smaller gear
%
% output:   Sigma = surface stress.
%
% Warning: Error most likely due to the interpolation of Km.
%
% All values must be constants, no matricies allowed. Use the following to
% get around this rule though:
%
% use matix for Qv, Vt, or other variable, for example:
%
% Quality = [5, 6, 7, 8];
% len = length(Quality);
% Sigma = ones(1,len);
% for coun = 1:len
%   Qv = Quality(coun);
%   Sigma(coun) = SurfaceStressSpurGear(Qv, Vt, FW, DrivingShock, ...
%    DrivenShock, size, I, Cp, surf, Wt, d);
% end



%%%%%% Dynamic Factor Kv %%%%%%%
if (Qv >= 6) && (Qv <= 11)
    B = (12 - Qv).^(2/3)./4;
    A = 50 + 56.*(1 - B);
    Kv = (A/(A + (Vt).^(1/2))).^B;
elseif (Qv <= 5)
    Kv = 50/(50 + sqrt(Vt));
else
    Kv = 1;
    fprintf('Error: Qv must be a number between 1 and 11.')
end


%%%%%% Load distribution factor %%%%%%%
if FW <= 2
    Km = 1.6;
elseif FW < 20 % Error most likely due to this interpolation.
    % Use two point formula for the two points of the form (FW, Km):
    % Points: (6, 1.7) & (9, 1.8).
    Km = ((1.8-1.7)/(9-6)).*((FW)-6) + 1.7;
else
    Km = 2.0;
end


%%%%%% Application factor %%%%%%%
switch DrivingShock
    case 'Uniform'
        switch DrivenShock
            case 'Uniform'
                Ka = 1;
            case 'Moderate Shock'
                Ka = 1.25;
            case 'Heavy Shock'
                Ka = 1.75;
            otherwise,
                fprintf('Error: must specify acceptable word for DrivenShock.')
        end
    case 'Light Shock'
        switch DrivenShock
            case 'Uniform'
                Ka = 1.25;
            case 'Moderate Shock'
                Ka = 1.5;
            case 'Heavy Shock'
                Ka = 2;
            otherwise,
                fprintf('Error: must specify acceptable word for DrivenShock.')
        end
    case 'Medium Shock'
        switch DrivenShock
            case 'Uniform'
                Ka = 1.5;
            case 'Moderate Shock'
                Ka = 1.75;
            case 'Heavy Shock'
                Ka = 2.25;
            otherwise,
                fprintf('Error: must specify acceptable word for DrivenShock.')
        end
    otherwise,
        fprintf('Error: must specify acceptable word for DrivingShock.')
        Ka = 1;
end


%%%%%% Size %%%%%%%
if  strcmp(size.shape, 'number')
    Ks = size.dim;
elseif strcmp(size.shape,'rectangle')
    A_95 = 0.05.*size.dim(1,1).*size.dim(1,2);

    d_equi = (A_95/0.0766)^(1/2);
    if d_equi <= 0.3
        Ks = 1;
    elseif (d_equi > 0.3) && (d_equi <= 10)
        Ks = 0.869.*d_equi.^(-0.097);
    else
        fprintf('Error: Equivalent diameter is too large.')
        Ks = 1;
    end

elseif strcmp(size.shape,'circle')
    if strcmp(size.rotation,'rotating')
        A_95 = 0.0766.*(size.dim).^2;
    elseif strcmp(size.rotation,'nonrotating')
        A_95 = 0.010462.*(size.dim).^2;
    else
        error('CorrectedStrength:WrongWord', 'Must specify either "rotating" or "nonrotating" only for size.rotation.')
    end

    d_equi = (A_95/0.0766)^(1/2);
    if d_equi <= 0.3
        Ks = 1;
    elseif (d_equi > 0.3) && (d_equi <= 10)
        Ks = 0.869.*d_equi.^(-0.097);
    else
        fprintf('Error: Equivalent diameter is too large.')
        Ks = 1;
    end

else
    error('Ks:WrongWord', 'must specify acceptable word for size.shape.')
end


%%%%%% Surface %%%%%%%
Cf = surf;



Cv = Kv; Cm = Km; Ca = Ka; Cs = Ks;

%%%%%% output %%%%%%%
Sigma = Cp.*sqrt((Wt)./((FW).*I.*d).*(Ca.*Cm)./(Cv).*Cs.*Cf);