function xx=griddef(x0,xf,npoints,xdensity,vdensity)
% function xx=griddef(x0,xf,npoints,xdensity,vdensity)
% x0: initial point
% xf: final point
% npoints: number of points
% [xdensity,vdensity]: function of point density (will be normalized
% to npoints)

% denser grid
xdens = linspace(x0,xf,1000*npoints);
vdens = interp1(xdensity,vdensity,xdens);

xequi = linspace(x0,xf,npoints);
% integrate
integ = cumtrapz(vdens)/trapz(vdens);
% reproject
iia=iround(integ,xequi); 
xx = xdens(iia);


% Previous version, discarded
% x0: initial point
% xf: final point
% npoints: number of points
% rhodens: location of densification
% wdens: width of densification
% frac: fraction of equidistant points
% % dense grid for gaussian
% xdens = linspace(x0,xf,1000*npoints);
% gauss=exp(-(xdens-rhodens).^2/wdens*4);
% 
% % iterate until right fraction obtained
% for c = 0:0.01:1
%     cgauss = max(gauss,c);
%     igau=find(gauss>c);
%     intcgauss = trapz(xdens(igau),cgauss(igau));
%     intgauss = cumtrapz(xdens,cgauss); % integrate
%     fracgauss = intcgauss/intgauss(end);
%     if fracgauss<(1-frac); break; end
% end
% xequi=linspace(x0,xf,npoints);
% % reproject
% iia=iround(intgauss/intgauss(end),xequi); 
% xx = xdens(iia);
