function V = sphericalChargePotential( Q, R, r )

% Calculates the electrostatic potential of a homogeneous spherical 
% distribution of charge Q and radius R at a distance r, given by
%   V(r) = K*Q*(3*R^2-r^2)/(2*R^3)  inside  (r<R)
%        = K*Q/r                    outside (r>R)
% where K = 9e9 J*m^2/C^2 is Coulomb's constant
% Input:
%   Q  : total charge (C)
%   R  : radius of spherical charge distribution (m)
%   r  : distance(s) to center of charge (m)
% Output:
%   V  : electrostatic potential at distance(s) r (V)
% Usage:
%   V = sphericalChargePotential( -3.e-9, 0.2, 0.1 )
% Author: J.M.Soler. Oct.2010, Sep.2015

K = 9.e9;   % Coulomb constant in J*m^2/C^2

nr = numel(r);
V = zeros(size(r));
% V = r*0;            % another way to do the same
for ir = 1:nr
    if r(ir)<R
       V(ir) = K * Q/2/R^3 * (3*R^2 - r(ir)^2);  % V inside R
    else
       V(ir) = K * Q / r(ir);                    % V outside R
    end
end

% A faster way. The +1.e-50 is to avoid a division by zero.
% inside = (r<R);
% V = inside .* K*Q/2/R^3*(3*R^2-r.^2) + ~inside .* K*Q./(r+1.e-50);

end % of function
