function d = findTubeD(t,dmax,P,Leff,E,Sy)
% findTubeD  diameter of round tube that does not buckle under an applied load
%
% Synopsis:  d = findTubeD(t,dmax,P,Leff,E,Sy)
%
% Input: t = desired wall thickness
%        dmax = maximum allowable pipe diameter
%        P = applied load
%        Leff = effective length of the tubing -- real length modified
%               by a factor to account for boundary conditions
%        E = modulus of elasticity of tubing material
%        Sy = yield stress of the tubing material

opts = optimset('Display','off');
d = fzero(@tubeDpD,dmax*[0.1 1],opts,t,P,Leff,E,Sy);

  % --- Nested function to define error in f(d) when trying to find d
  %     tubeDpD is visible only to code in this (findTubeD) mfile
  function dP = tubeDpD(d,t,P,Leff,E,Sy)

    ro = d/2;
    ri = ro - t;
    I = 0.25*pi*(ro^4 - ri^4);
    A = pi*(ro^2 - ri^2);
    k = 0.25*sqrt(ro^2 + ri^2);

    Sr = Leff/sqrt(I/A);     %  slenderness ratio
    if Sr > pi*sqrt(2*E/Sy)
      Pcr = pi^2 * A*E/Sr^2;
    else
      Pcr = A*(Sy - (1/E)*(Sy*Sr/(2*pi))^2 );
    end
    dP = Pcr - P;
  end

end