%  This script computes the sizes of a shaft that are necessary to for a
%  2.5 safety factor.  Fatigue, stress concentration, local yielding, and a
%  table of bearing, gear, and sheave sizes are considered.

%  Open output file

Out = fopen('ShaftSize.txt', 'w');

%  Input data for the problem

Sut = 400;          %  Ultimate tensile strength in MPa
Sy = 250;           %  Yield strength in MPa
Fillet = 0.5;       %  Step fillet size in mm
Nmin = 2.5;         %  Minimum Required safety factor
RPM = 1000;         %  Shaft speed in RPM
Power = 2000;       %  Average power transmitted in watts
Variance = 0.5;     %  Variance in power transmitted
DIAgear = 150;      %  Gear diameter in mm
DIAsheave = 125;    %  Sheave diameter in mm

DiaAC = 12;         %  Initial diameter at points A and C
DiaB = 15;          %  Initial diameter at point B
DiaD = 10;          %  Initial diameter at point D

fprintf(Out, '\n    Shaft Design    \n\n');

fprintf(Out, 'Ultimate tensile stress = %8.1f MPa\n', Sut);
fprintf(Out, 'Yield stress = %8.1f MPa\n', Sy);
fprintf(Out, 'Fillet Radius = %8.3f mm\n', Fillet);
fprintf(Out, 'Minimum Safety Factor = %8.2f\n', Nmin);
fprintf(Out, 'Rotational velocity = %8.2f\n', RPM);
fprintf(Out, 'Average Power Transmitted = %8.2f\n', Power);
fprintf(Out, 'Power variance = %8.3f\n', Variance);
fprintf(Out, 'Gear Diameter = %8.2f\n', DIAgear);
fprintf(Out, 'Sheave Diameter = %8.2f\n', DIAsheave);

Done = 0;           %  Iteration control

%  Shaft size Table

AllowableShaftSize = [10, 12, 15, 17, 20, 25, 30, 35, 40, 45, 50, 55];

%  Compute part of the Correction factor for Se

C = Cload('bending') * Csurf_Machined_m(Sut) * Crelib(99);

%  1.  Compute the torque created by the power transmitted from the gear to
%      the sheave.

W = RPM * 2 * pi / 60;          %  rotational velocity in radians / second
TorqueMean = 1000 * Power / W;  %  Mean torque in Newton-millimeters
TorqueAlt = Variance * 1000 * Power / W; %  Alternating torque in Newton-millimeters

%  2.  Compute the tangential forces on the sheave

FsheaveXmean = 1.5 * TorqueMean / (DIAsheave / 2);   %  Force in Newtons
FsheaveXalt  = 1.5 * TorqueAlt  / (DIAsheave / 2);   %  Force in Newtons

%  3.  Compute the tangential and normal forces on the gear.

FgearYmean = TorqueMean / (DIAgear / 2);
FgearXmean = FgearYmean * tand(20);

FgearYalt = TorqueAlt / (DIAgear / 2);
FgearXalt = FgearYalt * tand(20);

%  4.  Compute the moments at A, B, C, and D in both the X-Z and Y-Z planes

LoadsXmean = [85, FgearXmean; 195, FsheaveXmean];
LoadsXalt = [85, FgearXalt; 195, FsheaveXalt];
LoadsYmean = [85, FgearYmean];
LoadsYalt = [85, FgearYalt];

%           Compute the moments at A

[MXmean, R1, R2] = Shaft_Moment(15, LoadsXmean, 0, 160);
[MXalt, R1, R2] = Shaft_Moment(15, LoadsXalt, 0, 160);
[MYmean, R1, R2] = Shaft_Moment(15, LoadsYmean, 0, 160);
[MYalt, R1, R2] = Shaft_Moment(15, LoadsYalt, 0, 160);
MomentAmean = sqrt(MXmean^2 + MYmean^2);
MomentAalt  = sqrt(MXalt^2 + MYalt^2);

%           Compute the moments at B

[MXmean, R1, R2] = Shaft_Moment(56, LoadsXmean, 0, 160);
[MXalt, R1, R2] = Shaft_Moment(56, LoadsXalt, 0, 160);
[MYmean, R1, R2] = Shaft_Moment(56, LoadsYmean, 0, 160);
[MYalt, R1, R2] = Shaft_Moment(56, LoadsYalt, 0, 160);
MomentBmean = sqrt(MXmean^2 + MYmean^2);
MomentBalt = sqrt(MXalt^2 + MYalt^2);

%           Compute the moments at C

[MXmean, R1, R2] = Shaft_Moment(145, LoadsXmean, 0, 160);
[MXalt, R1, R2] = Shaft_Moment(145, LoadsXalt, 0, 160);
[MYmean, R1, R2] = Shaft_Moment(145, LoadsYmean, 0, 160);
[MYalt, R1, R2] = Shaft_Moment(145, LoadsYalt, 0, 160);
MomentCmean = sqrt(MXmean^2 + MYmean^2);
MomentCalt = sqrt(MXalt^2 + MYalt^2);

%           Compute the moments at D

[MXmean, R1, R2] = Shaft_Moment(185, LoadsXmean, 0, 160);
[MXalt, R1, R2] = Shaft_Moment(185, LoadsXalt, 0, 160);
[MYmean, R1, R2] = Shaft_Moment(185, LoadsYmean, 0, 160);
[MYalt, R1, R2] = Shaft_Moment(185, LoadsYalt, 0, 160);
MomentDmean = sqrt(MXmean^2 + MYmean^2);
MomentDalt = sqrt(MXalt^2 + MYalt^2);

%   ***** START THE LOOP LOOKING FOR CORRECT DIAMETERS *****

while (Done == 0)
    
    %  5.  Compute the alternating and mean nominal stresses due to bending at
    %      each location
    
    I = pi * DiaAC^4 / 64;
    SaltA = MomentAalt * (DiaAC / 2) / I;
    SmeanA = MomentAmean * (DiaAC / 2) / I;
        
    I = pi * DiaB^4 / 64;
    SaltB = MomentBalt * (DiaB / 2) / I;
    SmeanB = MomentBmean * (DiaB / 2) / I;
        
    I = pi * DiaAC^4 / 64;
    SaltC = MomentCalt * (DiaAC / 2) / I;
    SmeanC = MomentCmean * (DiaAC / 2) / I;
        
    I = pi * DiaD^4 / 64;
    SaltD = MomentDalt * (DiaD / 2) / I;
    SmeanD = MomentDmean * (DiaD / 2) / I;
        
    %  6.  Compute the alternating and mean shear stresses at each location due
    %      to the torque
    
    TaltA = 0.0;
    TmeanA = 0.0;
    
    TaltB = 0.0;
    TmeanB = 0.0;
    
    J = pi * DiaAC ^ 4 / 32;
    TaltC = TorqueAlt * (DiaAC / 2) / J;
    TmeanC = TorqueMean * (DiaAC / 2) / J;
    
    J = pi * DiaD ^ 4 / 32;
    TaltD = TorqueAlt * (DiaD / 2) / J;
    TmeanD = TorqueMean * (DiaD / 2) / J;
    
    %  7.  Compute the endurance stress at teach location
    
    SeAtAC = 0.5 * Sut * C * Csize_Round_Rotating_m(DiaAC);
    SeAtB  = 0.5 * Sut * C * Csize_Round_Rotating_m(DiaB);
    SeAtD  = 0.5 * Sut * C * Csize_Round_Rotating_m(DiaD);
    
    %  8.  Compute the static stress Concentration Factors for both bending and
    %  shear.
    
    KtA = Stress_Step_Bending(DiaB + 2, DiaAC, Fillet);
    KtsA = Stress_Step_Torsion(DiaB + 2, DiaAC, Fillet);
    
    KtB = Stress_Step_Bending(DiaB + 2, DiaB, Fillet);
    KtsB = Stress_Step_Torsion(DiaB + 2, DiaB, Fillet);
    
    KtC = Stress_Step_Bending(DiaB, DiaAC, Fillet);
    KtsC = Stress_Step_Torsion(DiaB, DiaAC, Fillet);
    
    KtD = Stress_Step_Bending(DiaAC, DiaD, Fillet);
    KtsD = Stress_Step_Torsion(DiaAC, DiaD, Fillet);
    
    %  9.  Reduce the stress concentration factors for fatigue
    
    KfA = Stress_Reduce_Bend_m(KtA, Fillet, Sut);
    KfsA = Stress_Reduce_Torque_m(KtsA, Fillet, Sut);
    
    KfB = Stress_Reduce_Bend_m(KtB, Fillet, Sut);
    KfsB = Stress_Reduce_Torque_m(KtsB, Fillet, Sut);
    
    KfC = Stress_Reduce_Bend_m(KtC, Fillet, Sut);
    KfsC = Stress_Reduce_Torque_m(KtsC, Fillet, Sut);
    
    KfD = Stress_Reduce_Bend_m(KtD, Fillet, Sut);
    KfsD = Stress_Reduce_Torque_m(KtsD, Fillet, Sut);
    
    %  9.  Compute the alternating and mean stresses at each location - we
    %  combind shear and torsion with von Mises stress
    
    VMaltA = VM_Stress_Alt(KfA, SaltA, 1.0, 0.0, 1.0, TaltA);
    VMmeanA = VM_Stress_Mean (KfA, SmeanA, 1.0, 0.0, 1.0, TmeanA);
    
    VMaltB = VM_Stress_Alt(KfB, SaltB, 1.0, 0.0, 1.0, TaltB);
    VMmeanB = VM_Stress_Mean (KfB, SmeanB, 1.0, 0.0, 1.0, TmeanB);
    
    VMaltC = VM_Stress_Alt(KfC, SaltC, 1.0, 0.0, 1.0, TaltC);
    VMmeanC = VM_Stress_Mean (KfA, SmeanC, 1.0, 0.0, KfsC, TmeanC);
    
    VMaltD = VM_Stress_Alt(KfD, SaltD, 1.0, 0.0, 1.0, TaltD);
    VMmeanD = VM_Stress_Mean (KfD, SmeanD, 1.0, 0.0, KfsD, TmeanD);
    
    %  10.  Compute the safety Factor using Langer's method
    
    %           Location A
    
    NlangerA = Langer(VMaltA, VMmeanA, Sy);
    NgoodmanA = Goodman(VMaltA, VMmeanA, SeAtAC, Sut);
    N = min(NlangerA, NgoodmanA);
    if (N < Nmin)
        DiaAC = NextSize(DiaAC, AllowableShaftSize);
        if (DiaB < DiaAC + 2)
            DiaB = NextSize(DiaB, AllowableShaftSize);
        end
        
        continue;
    end
    
    %           Location B
    
    NlangerB = Langer(VMaltB, VMmeanB, Sy);
    NgoodmanB = Goodman(VMaltB, VMmeanB, SeAtAC, Sut);
    N = min(NlangerB, NgoodmanB);
    if (N < Nmin)
        [DiaB] = NextSize(DiaB, AllowableShaftSize);
        
        continue;
    end
    
    %           Location C
    
    NlangerC = Langer(VMaltC, VMmeanC, Sy);
    NgoodmanC = Goodman(VMaltC, VMmeanC, SeAtAC, Sut);
    N = min(NlangerC, NgoodmanC);
    if (N < Nmin)
        [DiaAC] = NextSize(DiaAC, AllowableShaftSize);
        if (DiaB < DiaAC + 2)
            [DiaB] = NextSize(DiaB, AllowableShaftSize);
        end
        
        continue;
    end
    
    %           Location D
    
    NlangerD = Langer(VMaltD, VMmeanD, Sy);
    NgoodmanD = Goodman(VMaltD, VMmeanD, SeAtAC, Sut);
    N = min(NlangerD, NgoodmanD);
    if (N < Nmin)
        [DiaD] = NextSize(DiaD, AllowableShaftSize);
        if (DiaAC < DiaD + 2)
            [DiaAC] = NextSize(DiaB, AllowableShaftSize);
        end
        
        if (DiaB < DiaAC + 2)
            [DiaB] = NextSize(DiaB, AllowableShaftSize);
        end
        
        continue;
        
    else
        Done = 1;    
        
    end
    
    %  We have the Langer and modified Goodman safety factors but we need
    %  the Gerber and ASME
    
    NgerberA = Gerber(VMaltA, VMmeanA, SeAtAC, Sut);
    NgerberB = Gerber(VMaltB, VMmeanB, SeAtB, Sut);
    NgerberC = Gerber(VMaltC, VMmeanC, SeAtAC, Sut);
    NgerberD = Gerber(VMaltD, VMmeanD, SeAtD, Sut);
    
    NasmeA = ASME (VMaltA, VMmeanA, SeAtAC, Sut);
    NasmeB = ASME (VMaltB, VMmeanB, SeAtB, Sut);
    NasmeC = ASME (VMaltC, VMmeanC, SeAtAC, Sut);
    NasmeD = ASME (VMaltD, VMmeanD, SeAtD, Sut);
    
    %  Print the results
    
    fprintf(Out, '                         A        B        C        D\n');
    fprintf(Out, 'Mean Moments (N-mm)%8.0f %8.0f %8.0f %8.0f\n',MomentAmean, MomentBmean, MomentCmean, MomentDmean);
    fprintf(Out, 'Alt Moments  (N-mm)%8.0f %8.0f %8.0f %8.0f\n',MomentAalt, MomentBalt, MomentCalt, MomentDalt);
    fprintf(Out, 'Alt Stress (MPa)   %8.1f %8.1f %8.1f %8.1f\n',SaltA, SaltB, SaltC, SaltD);
    fprintf(Out, 'Mean Stress (MPa   %8.1f %8.1f %8.1f %8.1f\n',SmeanA, SmeanB, SmeanC, SmeanD);
    fprintf(Out, 'Alt Shear (MPa)    %8.1f %8.1f %8.1f %8.1f\n',TaltA, TaltB, TaltC, TaltD);
    fprintf(Out, 'Mean Shear (MPa)   %8.1f %8.1f %8.1f %8.1f\n',TmeanA, TmeanB, TmeanC, TmeanD);
    fprintf(Out, 'Se (MPa)           %8.0f %8.0f %8.0f %8.0f\n',SeAtAC, SeAtB, SeAtAC, SeAtD);
    fprintf(Out, 'Kt                 %8.3f %8.3f %8.3f %8.3f\n',KtA, KtB, KtC, KtD);
    fprintf(Out, 'Kts                %8.3f %8.3f %8.3f %8.3f\n',KtsA, KtsB, KtsC, KtsD);
    fprintf(Out, 'Kf                 %8.3f %8.3f %8.3f %8.3f\n',KfA, KfB, KfC, KfD);
    fprintf(Out, 'Kfs                %8.3f %8.3f %8.3f %8.3f\n'      ,KfsA, KfsB, KfsC, KfsD);
    fprintf(Out, 'VM alt (MPa)       %8.1f %8.1f %8.1f %8.1f\n',VMaltA, VMaltB, VMaltC, VMaltD);
    fprintf(Out, 'VM mean (MPa)      %8.1f %8.1f %8.1f %8.1f\n',VMmeanA, VMmeanB, VMmeanC, VMmeanD);
    fprintf(Out, 'Langer N           %8.3f %8.3f %8.3f %8.3f\n',NlangerA, NlangerB, NlangerC, NlangerD);
    fprintf(Out, 'Goodman N          %8.3f %8.3f %8.3f %8.3f\n',NgoodmanA, NgoodmanB, NgoodmanC, NgoodmanD);
    fprintf(Out, 'Gerber N           %8.3f %8.3f %8.3f %8.3f\n',NgerberA, NgerberB, NgerberC, NgerberD);
    fprintf(Out, 'ASME N             %8.3f %8.3f %8.3f %8.3f\n',NasmeA, NasmeB, NasmeC, NasmeD);
    fprintf(Out, 'Shaft Dia (mm)     %8.1f %8.1f %8.1f %8.1f\n',DiaAC, DiaB, DiaAC, DiaD);
    
    
end


    
    