%col_stab

clear

%user input----------------------------
yield_strength = 170e3; %psi
modulus = 28.5e6; %psi
max_dia = 1; %in
min_dia = 0.25; %in
delta_dia = .125; %in
max_l = 20; %in
min_l = 6; %in
delta_l = 0.25; %in
load = 1000; %lbs
%--------------------------------------

dia = min_dia:delta_dia:max_dia;
l = min_l:delta_l:max_l;


max_compression_force = load/(pi()*(min_dia/2)^2);
if(max_compression_force > yield_strength)
    fprintf('WILL FAIL IN COMPRESSION')
end

%from machinery handbook v26 p262 
denominator = 2*pi()^2*modulus; %eq3
for i=1:size(dia,2)
    for j = 1:size(l,2)
        A = pi() * ( dia(i)/2 )^2; %in^2
        r = dia(i)/4; %in
        Q = (yield_strength*l(j)^2)/denominator; %eq3
        P(i,j) = A*yield_strength*(1 - (Q/(4*r^2))); %lbs
    end
end

s = surf(l,dia,P);
s.EdgeColor = 'none';
view(2);
c = colorbar;
c.Label.String = 'Critical load (lbs)';
xlabel('length (in)');
ylabel('dia (in)');


