function [r1,alpha] = h3d_rho2index(s, rho)
	% r1 is the closest cell
	% alpha between -0.5 and 0.5 specifies the distance from the center of the cell
	% If r1 == 0, out of the buffer.

	if (rho < s.rho_min) || (rho > s.rho_max) 
		r1 = 0;
		alpha = 0;
	else
		% x belongs to [0, n)
		x = s.rho_ncells * (rho-s.rho_min) / (s.rho_max-s.rho_min);
		ra = floor( x );	
		alpha = (ra+0.5)-x;
		
		if abs(alpha) > 0.5
			error('rho = %f  x = %f  ra = %f  alpha = %f \n',rho,x,ra,alpha);
		end

		% stupid matlab counts from 1
		r1 = ra + 1;
		if (r1 > s.rho_ncells) || (r1<1)
			fprintf('ops ');
		end
	end