% APPENDIX B 

%%%% MATT STOFLETH | 10/21/15 %%%%



clear all

% INPUTS 
R2 = 2;  % CRANK ARM (in)
R3 = 8;  % CON ROD (in)
OM_2 = (50*2*pi)/60; % RAD/SEC
G = 386.09;  % IN/SEC^2
M_slider = 10;  % LB
M_shock = 10;   % LB 

% LINK PARAMETERS (in)
Dshaft = 1.2 ;  
FRshaft = .08 ;
Ashaft = (.25 * pi * (Dshaft^2));

Ratio = 3;
Warm = 2.5 ;   % b
THarm = Warm / Ratio ;  % c
Aarm = THarm * Warm;
% ALPHA AND BETA FOR b/c for ratio of 3
Alp_1 = 0.267;
Alp_2 = 0.355;
Beta = 0.263;
I_X = (Warm*(THarm^3))/12;
I_Y = ((Warm^3)*THarm)/12;

Dpin = 1 ; 
FRpin = .08 ;
Apin = (.25 * pi * (Dpin^2));

THrod = 1.5 ;
Wrod = 1 ; 
Arod = THrod * Wrod;
FRrod = .08 ;
ETrod = .6 ;  
ArodE = FRrod * ETrod; 

Lshaft = .5 - (.5  * THarm) ;   
Larm = 2 - (0.5 * Dshaft) + (0.5 * Dpin);
Lpin = FRpin + Wrod ;  
HDrod = Dpin ; 
Lrod = 8 + HDrod + (2 * ETrod) ; 

% Material Properties (psi)
Sut_S = 85000;
Sut_A = Sut_S;
Sut_P = Sut_S;
Sut_R = 85000;
Sut_RE = Sut_R;

Sy_S = 71000;
Sy_A = Sut_S;
Sy_P = Sy_S;
Sy_R = 71000;
Sy_RE = Sy_R;

% Se Calculation

Se_prime_S = 0.5 * Sut_S;
Se_prime_A = 0.5 * Sut_A;
Se_prime_P = 0.5 * Sut_P;
Se_prime_R = 0.5 * Sut_R;

ka_S = 0.8319;
ka_A = 0.8319;
ka_P = 0.8319;
ka_R = 0.8319;

kb_S = 0.879*(Dshaft^-0.107);
kb_A = 1;
kb_P = 0.879*(Dpin^-0.107);
kb_R = 1;

kc = 1; % Accounted for in Von Mises stress calc. 
kd = 1; % Assume room temp
ke = 0.897; % 90 percent reliability



%%%%%%%%%%%%%%%%%%%%%%%%%%%% LINK KINEMATICS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

i = 1; 

for A2 = 0:1:360;
    
    % KINEMATICS
    TH_2 = A2 * (pi/180);
    TH_3 = asin((-R2*sin(TH_2))/R3);
    R4 = R2*cos(TH_2) + R3*cos(TH_3);
    OM_3 = -((R2*cos(TH_2)*OM_2)/(R3*cos(TH_3))); 
    R4_D = -(R2*sin(TH_2)*OM_2) - (R3*sin(TH_3)*OM_3);
    ALPHA_3 = ((R2*sin(TH_2)*(OM_2^2)) + (R3*sin(TH_3)*(OM_3^2))) / (R3*cos(TH_3));
    R4_DD = (-R2*cos(TH_2)*(OM_2^2)) - (R3*cos(TH_3)*(OM_3^2)) - (R3*sin(TH_3)*ALPHA_3);
    R4_DDG = R4_DD/G;
    
    R2x = R2*cos(TH_2);
    R2y = R2*sin(TH_2);
    R3x = R3*cos(TH_3);
    R3y = R3*sin(TH_3);
    
    
    % FORCE ANALYSIS
    if R4_D <= -0.3
        F_SH = (32.5 * R4_D) - 105;
    elseif R4_D >= 0.3
        F_SH = (19.8 * R4_D) + 77;
    else 
        F_SH = 285.7 * R4_D;
    end
    
    % X Components
    F34x = (M_slider * R4_DDG) + F_SH;
    F43x = -1 * F34x; 
    F23x = -1 * F43x;
    F32x = F43x; 
    F12x = F23x;
    
    % Y Components
    F43y = (R3y / R3x) * F43x;
    F34y = -1 * F43y;
    F14y = M_slider - F34y;
    F23y = -1 * F43y;
    F32y = F43y;
    F12y = F23y;
    
    Tor = (R2y * F32x) - (R2x * (-1 * F23y));
    
    
    % CONVERSION TO ROTATING SYSTEM
    F12X = (F12x * cos(TH_2)) + (F12y * sin(TH_2));
    F12Y = (-1 * F12x * sin(TH_2)) + (F12y * cos(TH_2));
    
    F32X = (F32x * cos(TH_2)) + (F32y * sin(TH_2));
    F32Y = (-1 * F32x * sin(TH_2)) + (F32y * cos(TH_2));
    
    F23X = (F23x * cos(TH_3)) + (F23y * sin(TH_3));
    F23Y = (-1 * F23x * sin(TH_3)) + (F23y * cos(TH_3));
    
    F43X = (F43x * cos(TH_3)) + (F43y * sin(TH_3));
    F43Y = (-1 * F43x * sin(TH_3)) + (F43y * cos(TH_3));
    
    
    % MAGNITUDES
    F12 = sqrt((F12x^2) + (F12y^2));
    F23 = sqrt((F23x^2) + (F23y^2));
    F32 = sqrt((F32x^2) + (F32y^2));
    F34 = sqrt((F34x^2) + (F34y^2));
    

    
    % VECTOR FORMATION
    vector2(1,i) = TH_2 * (180/pi);
    vector3(1,i) = TH_3 * (180/pi);
    vector4(1,i) = R4;
    vector5(1,i) = OM_3;
    vector6(1,i) = R4_D;
    vector7(1,i) = ALPHA_3;
    vector8(1,i) = R4_DD;
    vector9(1,i) = F_SH;
    vector10(1,i) = Tor;
    vector11(1,i) = F43x;
    vector12(1,i) = F23x;
    vector13(1,i) = F43y;
    vector14(1,i) = F23y;
    vector15(1,i) = R2x;
    vector16(1,i) = R2y;
    vector17(1,i) = F12;
    vector18(1,i) = F32;
    vector19(1,i) = F32x;
    vector20(1,i) = F12X;
    vector21(1,i) = F12Y;
    vector22(1,i) = F32X;
    vector23(1,i) = F32Y;
    vector24(1,i) = F23X;
    vector25(1,i) = F23Y;
    vector26(1,i) = F12x;
    vector27(1,i) = F12y;
    vectorF43X(1,i) = F43X;
    vectorF43Y(1,i) = F43Y;
    
    i = i+1;
    
end



%%%%%%%%%%%%%%%%%%%%%%%% CRITICAL PLANE ANALYSIS %%%%%%%%%%%%%%%%%%%%%%%%%%


% SHAFT %

INC = .001; % increment of link iteration (in)
L = 0; % starting position on link (in)

q = 1;

while L <= Lshaft
        
    % SHEAR FORCE
    VshaftX = vector20(1,46);
    vectorVsX(1,q) = VshaftX;
    
    VshaftY = vector21(1,46);
    vectorVsY(1,q) = VshaftY;
      
    % MOMENT
    MshaftX = VshaftY * ((Lshaft + THarm + Lpin) - L) ; % About X
    vectorMsX(1,q) = -1 * MshaftX;
    
    MshaftY = VshaftX * ((Lshaft + THarm + Lpin) - L);
    vectorMsY(1,q) = -1 * MshaftY;
    
    % TORQUE
    Tshaft = (2 + (.5 * Dpin)) * VshaftY;  % At 45 degrees
    vectorTs(1,q) = Tshaft;
    
    vectorLs(1,q) = L;
    
    % ITERATION
    q = q + 1;
    L = L + INC;
       
end


% PIN %

INC = .001; % increment of link iteration (in)
L = 0; % starting position on link (in)

q = 1;

while L <= Lpin
        
    % SHEAR FORCE
    VpinX = vector20(1,46);
    vectorVpX(1,q) = VpinX;
    
    VpinY = vector21(1,46);
    vectorVpY(1,q) = VpinY;
      
    % MOMENT
    MpinX = VpinY * (Lpin - L); % About X
    vectorMpX(1,q) = -1 * MpinX;
    
    MpinY = VpinX * (Lpin - L); % About Y
    vectorMpY(1,q) = -1 * MpinY;
   
    vectorLp(1,q) = L;
    
    % ITERATION
    q = q + 1;
    L = L + INC;
         
end


% ARM ALONG Z %

INC = .001; % increment of link iteration (in)
L = 0; % starting position on link (in)

q = 1;

while L <= THarm
    
    % SHEAR FORCE 
    VarmX = vector20(1,46);
    vectorVaX1(1,q) = VarmX;
    
    VarmY = vector21(1,46);
    vectorVaY1(1,q) = VarmY;

    % MOMENT
    MarmY = VarmX * ((THarm + Lpin) - L); % About Y
    vectorMaY(1,q) = -1 * MarmY;
    
    % TORSION 
    Tarm = VarmY * (Lpin + (.5 * THarm));
    vectorTa(1,q) = -1 * Tarm;
    
    vectorLaz(1,q) = L;
    
    % ITERATION
    q = q + 1;
    L = L + INC;
    
end


% ARM ALONG X %

INC = .001; % increment of link iteration (in)
L = 0; % starting position on link (in)

q = 1;

while L <= Larm
    
    % SHEAR FORCE (USE ABOVE)

    % SHEAR FORCE 
    VarmX2 = vector20(1,46);
    vectorVaX2(1,q) = VarmX2;
    
    VarmY2 = vector21(1,46);
    vectorVaY2(1,q) = VarmY2;
    
    %MOMENT
    MarmX =  -1 * VarmY2 * (Larm - L); % About X
    vectorMaX(1,q) = MarmX;
    
    vectorLaX(1,q) = L;
    
    % ITERATION
    q = q + 1;
    L = L + INC;
    
end


%%%%%%%%%%%%%%%%%%%%%%%% CRITICAL POINT ANALYSIS %%%%%%%%%%%%%%%%%%%%%%%%%%


% SHAFT CP %

g = 1;

for PHI = 0:1:360; % PHI
     

   w = 1;
    for THETA_2 = 0:1:360 % THETA 2
        
        % NORMAL STRESS
    
        c = (Dshaft/2)*sind(180 - PHI);
        
        ShaftN_MY = ((vector20(1,(THETA_2+1))*(THarm + Lpin)) * c) / ((pi*(Dshaft^4))/ 64);   % About Y
        ShaftN_MX = ((vector21(1,(THETA_2+1))*(THarm + Lpin)) * c) / ((pi*(Dshaft^4))/ 64);   % About X
        ShaftN_T = ShaftN_MX + ShaftN_MY; 
        
        ShaftS_Tz = (16 * vectorTs(1,10)) / (pi * (Dshaft^3));
        ShaftS_VY = (4/3) * (vector21(1,(THETA_2+1))/Ashaft) * abs(sind(PHI));
        ShaftS_VX = (4/3) * (vector20(1,(THETA_2+1))/Ashaft) * abs(cosd(PHI));
        ShaftS_T = ShaftS_Tz + ShaftS_VY + ShaftS_VX; 
        
        
        % INTO VECTORS
        vector28(1,w) = ShaftN_T;
        vector29(1,w) = ShaftS_T;
        
        Shaft_N_T_Max = max(vector28);
        Shaft_N_T_Min = min(vector28);
        Shaft_S_T_Max = max(vector29); 
        Shaft_S_T_Min = min(vector29);
        
        w = w + 1;
         
    end
    
    vector30(1,g) = Shaft_N_T_Max;
    vector31(1,g) = Shaft_S_T_Max;
    vector32(1,g) = Shaft_N_T_Min;
    vector33(1,g) = Shaft_S_T_Min;
     
    g = g + 1;
    
end


% PIN CP %

g = 1;

for PHI = 0:1:360; % PHI
     

   w = 1;
    for THETA_2 = 0:1:360 % THETA 2
        
        % NORMAL STRESS
    
        c = (Dpin/2)*sind(180 - PHI);
        
        PinN_MY = ((vector20(1,(THETA_2+1))*Lpin) * c) / ((pi*(Dpin^4))/ 64);   % About Y
        PinN_MX = ((vector21(1,(THETA_2+1))*Lpin) * c) / ((pi*(Dpin^4))/ 64);   % About X
        PinN_T = PinN_MX + PinN_MY; 
        
        PinS_VY = (4/3) * (vector21(1,(THETA_2+1))/Apin) * abs(sind(PHI));
        PinS_VX = (4/3) * (vector20(1,(THETA_2+1))/Apin) * abs(cosd(PHI));
        PinS_T = PinS_VY + PinS_VX; 
        
        
        % INTO VECTORS
        vector34(1,w) = PinN_T;
        vector35(1,w) = PinS_T;
        
        Pin_N_T_Max = max(vector34);
        Pin_N_T_Min = min(vector34);
        Pin_S_T_Max = max(vector35);  
        Pin_S_T_Min = min(vector35);
        
        w = w + 1;
         
    end
 
    vector36(1,g) = Pin_N_T_Max;
    vector37(1,g) = Pin_S_T_Max;
    vector38(1,g) = Pin_N_T_Min;
    vector39(1,g) = Pin_S_T_Min;
     
    g = g + 1;
    
end 


% ARM CP %

g = 1;

for POS = 1:1:4 % Position in arm cross section looking from slider. 1 in positive z and progressing CCW.                     
    
    w = 1;
    for THETA_2 = 0:1:360 % THETA 2

        if POS == 1 
            ArmS_Tor = ((vector21(1,(THETA_2+1)))*(Lpin + (.5 * THarm))) / (Alp_1*Warm*(THarm^2));
            ArmN_MY = ((vector20(1,(THETA_2+1))*(Lpin))*(Warm/2))/I_X; % At point 1
        end
        
        if POS == 2
            ArmS_Tor = ((vector21(1,(THETA_2+1)))*(Lpin + (.5 * THarm))) / (Alp_2*Warm*(THarm^2));
            ArmN_MY = ((vector20(1,(THETA_2+1))*(Lpin + .5  * THarm))*(Warm/2))/I_X; % At point 2
        end
        
        if POS == 3
            ArmS_Tor = ((vector21(1,(THETA_2+1)))*(Lpin + (.5 * THarm))) / (Alp_1*Warm*(THarm^2));
            ArmN_MY = ((vector20(1,(THETA_2+1))*(Lpin + THarm))*(Warm/2))/I_X; % At point 4
        end
        
        if POS == 4
            ArmS_Tor = ((vector21(1,(THETA_2+1)))*(Lpin + (.5 * THarm))) / (Alp_2*Warm*(THarm^2));
            ArmN_MY = ((vector20(1,(THETA_2+1))*(Lpin + .5 * THarm))*(Warm/2))/I_X; % At point 4
        end
        
             
        ArmN_MZ = ((vector21(1,(THETA_2+1))*(Larm))*(Warm/2))/I_Y;
        ArmN_AX = (vector20(1, (THETA_2+1)))/Aarm; 
        ArmN_T = ArmN_MY + ArmN_MZ + ArmN_AX;
        ArmN_Bend = ArmN_MY + ArmN_MZ;  % sigma a bending
        
        ArmS_VY = (4/3) * (vector21(1,(THETA_2+1))/Aarm);
        ArmS_VX = (4/3) * (vector20(1,(THETA_2+1))/Aarm);
      
        ArmS_T = ArmS_Tor + ArmS_VY + ArmS_VX;
        
        % INTO VECTORS
        vector40(1,w) = ArmN_T;
        vector41(1,w) = ArmS_T;
        vector48(1,w) = ArmN_Bend;
        vector49(1,w) = ArmN_AX;
        
        Arm_N_T_Max = max(vector40);
        Arm_N_T_Min = min(vector40);
        Arm_S_T_Max = max(vector41);
        Arm_S_T_Min = min(vector41);
        Arm_N_Bend_Max = max(vector48);
        Arm_N_Bend_Min = min(vector48);
        Arm_N_AX_Max = max(vector49);
        Arm_N_AX_Min = min(vector49);
        
        w = w + 1;
        
    end
  
    vector42(1,g) = Arm_N_T_Max;
    vector43(1,g) = Arm_S_T_Max;
    vector44(1,g) = Arm_N_T_Min;
    vector45(1,g) = Arm_S_T_Min;
    vector50(1,g) = Arm_N_Bend_Max;
    vector51(1,g) = Arm_N_Bend_Min; 
    vector52(1,g) = Arm_N_AX_Max;
    vector53(1,g) = Arm_N_AX_Min;
    vectorg(1,g) = g;
    
    g = g + 1;
    
    
end


% CONNECTING ROD CP %

w = 1;

for THETA_2 = 0:1:360 % THETA 2
    
    RodN_AX = sqrt((((vectorF43X(1, (THETA_2+1)))^2) + ((vectorF43Y(1, (THETA_2+1)))^2))) / Arod; 
    RodEN_AX = sqrt((((-1*vectorF43X(1, (THETA_2+1)))^2) + ((vectorF43Y(1, (THETA_2+1)))^2))) / ArodE;
    
    % INTO VECTORS
    vector46(1,w) = RodN_AX;
    vector47(1,w) = RodEN_AX;
    
    Rod_N_AX_Max = max(vector46);
    Rod_N_AX_Min = min(vector46);
    RodE_N_AX_Max = max(vector47);
    RodE_N_AX_Min = min(vector47);
    
    w = w+1;
end


%%%%%%%%%%%%%%%%%%%%%%%% SHAFT STRESS COMPONENTS %%%%%%%%%%%%%%%%%%%%%%%%%%

altsigma_S = (vector30(1,91) - vector32(1,91))/2; 
meansigma_S = (vector30(1,91) + vector32(1,91))/2;
alttau_S = (vector31(1,91) - vector33(1,91))/2;
meantau_S = (vector31(1,91) + vector33(1,91))/2;

%%check table A-15-6
Dd_S = Warm/Dshaft;
rd_S = FRshaft/Dshaft; % ratio of fillet radius to shaft diameter
KT_S = 1.95; % SEE TABLE A-15-9 FOR VALUE - Stress concentration (bending)
KTS_S = 1.7; % SEE TABLE A-15-8 FOR VALUE - Stress Concentration (torsion)
sqrta_bendS = .246 - (3.08*(10^(-3))*Sut_S) + (1.51*(10^(-5))*Sut_S^2) - (2.67*(10^(-8))*(Sut_S^3));
sqrta_torsS = .190 - (2.51*(10^(-3))*Sut_S + (1.35*(10^(-5))*Sut_S^2) - (2.67*(10^(-8))*Sut_S^3));
Kf_S = 1 + (KT_S - 1)/(1 + (sqrta_bendS/sqrt(Dshaft/2)));
KfS_S = 1 + (KTS_S - 1)/(1 + (sqrta_torsS/sqrt(Dshaft/2)));
altsigmaprime_S = sqrt((Kf_S*altsigma_S)^2 + 3*(KfS_S*alttau_S)^2);
meansigmaprime_S = sqrt((Kf_S*meansigma_S)^2 + 3*(KfS_S*meantau_S)^2);


%%%%%%%%%%%%%%%%%%%%%%%%%% PIN STRESS COMPONENTS %%%%%%%%%%%%%%%%%%%%%%%%%%

altsigma_pin = (vector36(1,91) - vector38(1,91))/2;
meansigma_pin = (vector36(1,91) + vector38(1,91))/2;
alttau_pin = (vector37(1,91) - vector39(1,91))/2;
meantau_pin = (vector37(1,91) + vector39(1,91))/2;

Dd_P = Warm/Dpin; 
rd_P = FRpin/Dpin;%ratio of fillet radius to shaft diameter
KT_P = 1.8; %SEE TABLE A-15-12 FOR VALUE
sqrta_bendP = .246 - 3.08*10^(-3)*Sut_P + 1.51*10^(-5)*Sut_P^2 - 2.67*10^(-8)*Sut_P^3;
sqrta_torsP = .190 - 2.51*10^(-3)*Sut_P + 1.35*10^(-5)*Sut_P^2 - 2.67*10^(-8)*Sut_P^3;
Kf_P = 1 + (KT_P - 1)/(1 + (sqrta_bendP/sqrt(Dpin/2)));
altsigmaprime_P = sqrt((Kf_P*altsigma_pin)^2);
meansigmaprime_P = sqrt((Kf_P*meansigma_pin)^2);


%%%%%%%%%%%%%%%%%%%%%%%% ARM STRESS COMPONENTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%
altsigma_Arm_Bend = (vector50(1,3) - vector51(1,3))/2;
altsigma_Arm_AX = (vector52(1,3) - vector52(1,3))/2;
meansigma_Arm_Bend = (vector50(1,3) + vector51(1,3))/2;
meansigma_Arm_AX = (vector52(1,3) + vector52(1,3))/2;

alttau_Arm = (Arm_S_T_Max - Arm_S_T_Min)/2;
meantau_Arm = (Arm_S_T_Max + Arm_S_T_Min)/2;
altsigmaprime_A = sqrt((altsigma_Arm_Bend + (altsigma_Arm_AX/0.85))^2 + (3 * (alttau_Arm^2)));
meansigmaprime_A = sqrt((meansigma_Arm_Bend + (meansigma_Arm_AX))^2 + (3 * (meantau_Arm^2)));


%%%%%%%%%%%%%%%%%%%%%%%% ROD STRESS COMPONENTS %%%%%%%%%%%%%%%%%%%%%%%%%%%% 

altsigma_Rod_AX = (Rod_N_AX_Max - Rod_N_AX_Min)/2;
meansigma_Rod_AX = (Rod_N_AX_Max + Rod_N_AX_Min)/2;
alttau_Rod = 0;
meantau_Rod = 0;

Dd_R = ((2*ETrod)+Dpin)/THarm; 
rd_R = FRrod/THarm;  %ratio of fillet radius to shaft diameter
KT_R = 2.7; %SEE TABLE A-15-5 FOR VALUE

sqrta_bendR = .246 - 3.08*10^(-3)*Sut_R + 1.51*10^(-5)*Sut_R^2 - 2.67*10^(-8)*Sut_R^3;
sqrta_torsR = .190 - 2.51*10^(-3)*Sut_R + 1.35*10^(-5)*Sut_R^2 - 2.67*10^(-8)*Sut_R^3;
Kf_R = 1 + (KT_R - 1)/(1 + (sqrta_bendR/sqrt(THarm/2)));

altsigmaprime_R = sqrt((0 + (Kf_R*(altsigma_Rod_AX/0.85)))^2); 
meansigmaprime_R = sqrt((Kf_R*meansigma_Rod_AX)^2);


%%%%%%%%%%%%%%%%%%%%%%% ROD END STRESS COMPONENTS %%%%%%%%%%%%%%%%%%%%%%%%% 

altsigma_RodE_AX = (RodE_N_AX_Max - RodE_N_AX_Min)/2;
meansigma_RodE_AX = (RodE_N_AX_Max + RodE_N_AX_Min)/2;
alttau_RodE = 0;
meantau_RodE = 0;

dw_RE = Dpin/((2*ETrod)+Dpin);
hw_RE = (ETrod + (Dpin/2)) / ((2*ETrod)+Dpin);
KT_RE = 2.5; %SEE TABLE A-15-12 FOR VALUE

sqrta_axialRE = .246 - 3.08*10^(-3)*Sut_RE + 1.51*10^(-5)*Sut_RE^2 - 2.67*10^(-8)*Sut_RE^3;
Kf_RE = 1 + (KT_RE - 1)/(1 + (sqrta_axialRE/sqrt(ETrod/2)));
altsigmaprime_RE = sqrt((Kf_RE*(altsigma_RodE_AX/0.85))^2);
meansigmaprime_RE = sqrt((Kf_RE*meansigma_RodE_AX)^2);


%%%%%%%%%%%%%%%%%%%%%%%%%%% FAILURE CRITERION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Se_S = Se_prime_S * ka_S * kb_S * kc * kd * ke * Kf_S;
Se_A = Se_prime_A * ka_A * kb_A * kc * kd * ke;
Se_P = Se_prime_P * ka_P * kb_P * kc * kd * ke * Kf_P;
Se_R = Se_prime_R * ka_R * kb_R * kc * kd * ke * Kf_R;
Se_RE = Se_prime_R * ka_R * kb_R * kc * kd * ke * Kf_R;



%%% FACTORS OF SAFETY %%%

% SHAFT %

% Langer 
NSy_S1 = Sy_S / (altsigmaprime_S + meansigmaprime_S) % More conservative approach than the square root
% von Mises
NSy_S2 = Sy_S / sqrt((meansigma_S + altsigma_S)^2 + (3*((alttau_S + meantau_S)^2)))
% Mod-Good
N_S = 1 / ((altsigmaprime_S / Se_S) + (meansigmaprime_S / Sut_S))


% ARM %

% Langer
NSy_A = Sy_A / (altsigmaprime_A + meansigmaprime_A)
% Mod-Good
N_A = 1 / ((altsigmaprime_A / Se_A) + (meansigmaprime_A / Sut_A))

 
% PIN %

% Langer
NSy_P = Sy_P / (altsigmaprime_P + meansigmaprime_P)
% Mod-Good
N_P = 1 / ((altsigmaprime_P / Se_P) + (meansigmaprime_P / Sut_P))


% CONNECTING ROD MIDDLE

% Langer %
NSy_R = Sy_R / (altsigmaprime_R + meansigmaprime_R)
% Mod-Good
N_R = 1 / ((altsigmaprime_R / Se_R) + (meansigmaprime_R / Sut_R))


% CONNECTING ROD END

% Langer %
NSy_RE = Sy_RE / (altsigmaprime_RE + meansigmaprime_RE)
% Mod-Good
N_RE = 1 / ((altsigmaprime_RE / Se_RE) + (meansigmaprime_RE / Sut_RE))


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% KINEMATICS AND FORCES
% figure(9)
% plot(vector2, vector10);
% title('Crank Position vs. Slider Velocity');
% xlabel('Theta 2 (degrees)');
% ylabel('R4 Dot (in/sec)');
% axis([0 360 0 1000]);
% 
% figure(10)
% plot(vector2, vector11);
% title('Forces On Connecting Rod');
% xlabel('Crank Angle (degrees)');
% ylabel('Force (lb)');
% axis([0 360 -500 500]);
% hold on;
% plot(vector2, vector12);
% hold on;
% plot(vector2, vector13);
% hold on
% plot(vector2, vector14);
% hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% figure(11) % SHAFT
% subplot(4,1,1)
% plot(vectorLs, vectorVsY);
% title('Shaft Shear In Y');
% xlabel('Length of Shaft (in)');
% ylabel('Shaft Shear in Y (lbs)');
% axis([0 0.085 324 327]);
% subplot(4,1,2)
% plot(vectorLs, vectorMsX);
% title('Shaft Moment About X');
% xlabel('Length of Shaft (in)');
% ylabel('Shaft Moment (in-lbs)');
% axis([0 0.085 -650 -620]);
% subplot(4,1,3)
% plot(vectorLs, vectorVsX);
% title('Shaft Shear In X');
% xlabel('Length of Shaft (in)');
% ylabel('Shaft Shear in X (lbs)');
% axis([0 0.085 -228 -225]);
% subplot(4,1,4)
% plot(vectorLs, vectorMsY);
% title('Shaft Moment about Y');
% xlabel('Length of Shaft (in)');
% ylabel('Shaft Moment (in-lbs)');
% axis([0 0.085 430 460]);

% figure(12) % ARM IN Z 
% subplot(4,1,1)
% plot(vectorLaz, vectorVaY1);
% title('Arm Shear In Y')
% xlabel('Length of Arm in Z (in)');
% ylabel('Arm Shear in Y (lbs)');
% axis([0 0.85 324 327]);
% subplot(4,1,2)
% plot(vectorLaz, vectorTa);
% title('Arm Torque About X')
% xlabel('Length of Arm in Z (in)');
% ylabel('Arm Torque (in-lbs)');
% axis([0 0.85 -488 -485]);
% subplot(4,1,3)
% plot(vectorLaz, vectorVaX1);
% title('Arm Shear In X')
% xlabel('Length of Arm in Z (in)');
% ylabel('Arm Shear in X (lbs)');
% axis([0 0.85 -228 -225]);
% subplot(4,1,4)
% plot(vectorLaz, vectorMaY);
% title('Arm Moment About Y')
% xlabel('Length of Arm in Z (in)');
% ylabel('Arm Moment About Y (in-lbs)');
% axis([0 0.85 200 500]);

% figure(13) % ARM IN X
% subplot(2,1,1)
% plot(vectorLaX, vectorVaY2);
% title('Arm Shear In Y');
% xlabel('Length of Arm in X (in)');
% ylabel('Arm Shear in Y (lbs)');
% axis([0 1.9 324 327]);
% subplot(2,1,2)
% plot(vectorLaX, vectorMaX);
% title('Arm Moment About Z')
% xlabel('Length of Arm in X (in)');
% ylabel('Arm Moment About Z (in-lbs)');
% axis([0 1.9 -800 0]);

% figure(14) % PIN
% subplot(4,1,1)
% plot(vectorLp, vectorVpY);
% title('Pin Shear In Y');
% xlabel('Length of Pin (in)');
% ylabel('Pin Shear in Y (lbs)');
% axis([0 1.1 324 327]);
% subplot(4,1,2)
% plot(vectorLp, vectorMpX);
% title('Pin Moment About X');
% xlabel('Length of Pin (in)');
% ylabel('Pin Moment About X (in-lbs)');
% axis([0 1.1 -400 0]);
% subplot(4,1,3)
% plot(vectorLp, vectorVpX);
% title('Pin Shear In X');
% xlabel('Length of Pin (in)');
% ylabel('Pin Shear in X (lbs)');
% axis([0 1.1 -228 -225]);
% subplot(4,1,4)
% plot(vectorLp, vectorMpY);
% title('Pin Moment About Y');
% xlabel('Length of Pin (in)');
% ylabel('Pin Moment About Y (in-lbs)');
% axis([0 1.1 0 300]);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% figure(15)
% subplot(2,1,1)
% plot(vector2, vector30);
% title('Shaft Max Normal Stress');
% xlabel('PHI (degrees)');
% ylabel('Normal Stress (psi)');
% axis([0 360 0 8000]);
% subplot(2,1,2)
% plot(vector2, vector31);
% title('Shaft Max Shear Stress');
% xlabel('PHI (degrees)');
% ylabel('Shear Stress (psi)');
% axis([0 360 2600 3000]);

% figure(16)
% subplot(2,1,1)
% plot(vector2, vector36);
% title('Pin Max Normal Stress');
% xlabel('PHI (degrees)');
% ylabel('Normal Stress (psi)');
% axis([0 360 0 8000]);
% subplot(2,1,2)
% plot(vector2, vector37);
% title('Pin Max Shear Stress');
% xlabel('PHI (degrees)');
% ylabel('Shear Stress (psi)');
% axis([0 360 400 800]);
% 
% figure(17)
% subplot(2,1,1)
% scatter(vectorg, vector42,'k','x');
% title('Arm Max Normal Stress');
% xlabel('Point On Cross-Section');
% ylabel('Normal Stress (psi)');
% 
% subplot(2,1,2)
% scatter(vectorg, vector43,'k','x');
% title('Arm Max Shear Stress');
% xlabel('Point On Cross-Section');
% ylabel('Shear Stress (psi)');

% figure(18)
% subplot(2,1,1)
% plot(vector2, vector46);
% title('Rod Max Normal Stress');
% subplot(2,1,2)
% plot(vector2, vector47);
% title('Rod End Max Normal Stress');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(4)
subplot(3,1,1)
plot(vector2, vector4);
title('Crank Position vs. Slider Position');
xlabel('Theta 2 (degrees)');
ylabel('R4 (in)');
axis([0 360 6 10]);

subplot(3,1,2)
plot(vector2, vector6);
title('Crank Position vs. Slider Velocity');
xlabel('Theta 2 (degrees)');
ylabel('R4 Dot (in/sec)');
axis([0 360 -20 20]);

subplot(3,1,3)
plot(vector2, vector8);
title('Crank Position vs. Slider Acceleration');
xlabel('Theta 2 (degrees)');
ylabel('R4 Double Dot (in/sec^2)');
axis([0 360 -100 50]);


% ALL PLOTS TOGETHER%%%%%%%%%%%%%%%%%%%%%%%%%%

% figure (5)
% plot(vector2, vector4);
% hold on;
% plot(vector2, vector6);
% hold on;
% plot(vector2, vector8);
% hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% figure(1)
% 
% subplot(2,1,1)
% plot(vector2, vector3);
% title('Input Influence On TH_3');
% xlabel('Theta 2 (degrees)');
% ylabel('Theta 3 (degrees)');   
% 
% subplot(2,1,2)
% plot(vector2, vector4);
% title('Input Influence On R4');
% xlabel('Theta 2 (degrees)');
% ylabel('R4 (in)');
% 
% 
% figure(2)
% 
% subplot(2,1,1)
% plot(vector2, vector5);
% title('Input Influence On OM_3');
% xlabel('Theta 2 (degrees)');
% ylabel('OM_3 (rad/s)');   
% 
% subplot(2,1,2)
% plot(vector2, vector6);
% title('Input Influence On R4_D');
% xlabel('Theta 2 (degrees)');
% ylabel('R4 dot (in/sec)');
% 
% 
% figure(3)
% 
% subplot(2,1,1)
% plot(vector2, vector7);
% title('Input Influence On ALPHA_3');
% xlabel('Theta 2 (degrees)');
% ylabel('ALPHA 3 (rad/s^2)');   
% 
% subplot(2,1,2)
% plot(vector2, vector8);
% title('Input Influence On R4 Double Dot');
% xlabel('Theta 2 (degrees)');
% ylabel('R4_DD (in/sec^2)');