function s = lhsu(xmin,xmax,nsample)
% Latin Hypercube sampling 
% modified from Jasper Vrugt's SCEM-UA code

nvar=length(xmin);
ran=rand(nsample,nvar);
s=zeros(nsample,nvar);
for j=1: nvar
   idx=randperm(nsample);
   P =(idx'-ran(:,j))/nsample;
   s(:,j) = xmin(j) + P.* (xmax(j)-xmin(j));
end
