1;

clear all;
%%%%%%%%CALCULATED FORCES%%%%%%%%%%%
F_az = 588.60 ;

Rbz= 7.951587e+02;
Rcz= -4.486667e-01;
Rez= -9.968817e+01 ;
Rgz= -1.064218e+02;
Rby= -4.081000e+00;
Rcy= 1.632400e+01;
Rex= -6.588500e+01;
Rgx= -6.588500e+01;

%Old Values
%Rbz= 1.591013e+03  ;
%Rcz= 2.295387e+03  ;
%Rez= -1.559891e+03 ;
%Rgz= -1.737909e+03 ;
%Rby= -1.977733e+01 ;
%Rcy= 2.316773e+02  ;
%Rex= -1.017800e+03 ;
%Rgx= -1.090500e+03 ;
W1 =  206.11;
W2 =  131.77;
W3 =  13.243;
Wtg =  W3;
Wtw =  W2;
Wrg =  W1;
Wrw =  Wrg;
Wxg = Wtw;
Wxw = Wtg;
%W1 =  3297.8;
%W2 =  2108.3;
%W3 =  211.90;


%%%%%%%%parameters%%%%%%%%%%%
b = 0.05;
c = 0.20;
d = 0.25; %from 0.300;
rG = 0.090;
f = 0.04;
g = 0.08;
rW = 0.022;

%Winch shaft
x = linspace (0, d, 1000);
y1 = -F_az*H(x, 0) + Rbz*H(x, b) + Rcz*H(x, c) - Wrg*H(x, d);
y2 = -F_az*M(x, 0) + Rbz*M(x, b) + Rcz*M(x, c) - Wrg*M(x, d) - Wxg*rG*H(x, d);
figure(1)
plot (x, y1)
title ("Z Shear in Winch shaft");
xlabel ("Length");
ylabel ("Shear");
xlim([0, 0.26]);
grid on
figure(2)
plot (x, y2)
title ("Y Moment in Winch shaft");
xlabel ("Length");
ylabel ("Moment");
xlim([0, 0.26]);
grid on

y3 = Rby*H(x, b) + Rcy*H(x, c) - Wtg*H(x, d);
y4 = Rby*M(x, b) + Rcy*M(x, c) - Wtg*M(x, d);

figure(3)
plot (x, y3)
title ("Y Shear in Winch shaft");
xlabel ("Length");
ylabel ("Shear");
xlim([0, 0.26]);
grid on
figure(4)
plot (x, y4)
title ("Z Moment in Winch shaft");
xlabel ("Length");
ylabel ("Moment");
xlim([0, 0.26]);
grid on

%Motor shaft
x = linspace (0, g, 1000);

y5 = Rez*H(x, 0) + Wrw*H(x, f) + Rgz*H(x, g);
y6 = Rez*M(x, 0) + Wrw*M(x, f) + Rgz*M(x, g);

figure(5)
plot (x, y5)
title ("Z Shear in Motor shaft");
xlabel ("Length");
ylabel ("Shear");
xlim([0, 0.09]);
grid on
figure(6)
plot (x, y6)
title ("X Moment in Motor shaft");
xlabel ("Length");
ylabel ("Moment");
xlim([0, 0.09]);
grid on

y7 = Rex*H(x, 0) + Wtw*H(x, f) + Rgx*H(x, g);
y8 = Rex*M(x, 0) + Wtw*M(x, f) + Rgx*M(x, g);

figure(7)
plot (x, y7)
title ("X Shear in Motor shaft");
xlabel ("Length");
ylabel ("Shear");
xlim([0, 0.09]);
grid on
figure(8)
plot (x, y8)
title ("Z Moment in Motor shaft");
xlabel ("Length");
ylabel ("Moment");
xlim([0, 0.09]);
grid on

%%%% Max moment %%%%%%%%

MaxM_winch = 0;
MaxM_winch_x = 0;
for i = 1:1000
    if sqrt(y2(i)^2 + y4(i)^2) > MaxM_winch
        MaxM_winch = sqrt(y2(i)^2 + y4(i)^2);
        MaxM_winch_x = i*d/1000;
    end
end


MaxM_motor = 0;
MaxM_motor_x = 0;
for i = 1:1000
    if sqrt(y6(i)^2 + y8(i)^2) > MaxM_motor
        MaxM_motor = sqrt(y6(i)^2 + y8(i)^2);
        MaxM_motor_x = i*g/1000;
    end
end


%%%%% WINCH SHAFT DESIGN %%%%%%%%

rWinchShaft1 = 0.01; %mm to in
Ma = MaxM_winch; %Nm to lb.in
T_winch =  147.15; %Nm to lb.in
Tm = T_winch; %lb.in

%%%%%%%%%%%%% GOODMAN %%%%%%%%%%
conv = 0.1450377;
Sut = 686 ; %MPA--> kpsi
Sy = 490; %MPa --> kpsi
Nf = 2.5;
Lh = 6000; %hrs
%Ka = 2.7*Sut^-0.265; %machined surface finish,

%Kb = 0.879 * Sut^-0.107; %girth factor for 2in MAX shaft (51mm)
Ka = 4.51*Sut^-0.265 %machined surface finish,  Sut in MPa

Kb = 1.51 * (rWinchShaft1*2)^-0.107 %girth factor for shaft (>51mm)
Ke = 1;
Seprime = 0.5 * Sut;
Se = Seprime*Ka*Kb*Ke;

RTa= 0.246 - (3.08e-3)*Sut + (1.51e-5)*Sut^2 - 2.67e-8*Sut^3;
RTas= 0.190 - (2.51e-3) *Sut + (1.35e-5)*Sut^2 - 2.67e-8*Sut^3;

q = (1 + (RTa)/sqrt(rWinchShaft1))^(-1);
qs = (1 + (RTas)/sqrt(rWinchShaft1))^(-1);

%Kt, Kts from graph%%%%%%%%%%
Kt = 1.1; %guessed, need r/d and D/d
Kts =1.1; %guessed, need r/d and D/d

%Kf = 1 + q*(Kt - 1); %imperial
%Kfs = 1 + qs*(Kts - 1); %imperial

Kf = 1;
Kfs = 1;
dWinchShaft = ( (16*Nf/pi)*( 2*Kf*Ma/Se + sqrt(3)*Kfs*Tm/Sut ) )^(1/3)


function y = H(x, a)
    y = (x-a >= 0);
end

function y = M(x, a)
    y = ((x-a).*H(x,a));
end