function [r,drhodr] = drho_dr(geop)
r = geop.epsilon.*geop.Rgeom;
drhodr = 1./(geop.Rgeom.*geop.epsilonp + geop.Rgeomp.*geop.epsilon); 
% drho/dr = 1/(dr/drho), r(rho)=epsilon(rho)*Rgeom(rho)

return
