function y = truncatedNormalDraw(lbound,ubound,mu,sigma)
%{

mean is mu
variance is sigma

Call the function trandn.m to draw from a truncated normal
very carefully

if mu or sigma are vectors, then expand lengths of bounds
to do a bunch of draws efficiently in trandn

%}
muLength = length(mu);
sigmaLength = length(sigma);
if muLength>1 || sigmaLength > 1
    if muLength > 1 && sigmaLength > 1 && muLength ~= sigmaLength
        error('mean and variance different lengths');
    end
    lbound = repmat(lbound,max(muLength,sigmaLength,1));
    ubound = repmat(ubound,max(muLength,sigmaLength,1));
end

X = trandn((lbound-mu)./sqrt(sigma),(ubound-mu)./sqrt(sigma));
y = mu + sqrt(sigma).*X;

end

