Does have Mathematica a built-in function for write Dirac delta in spherical coordinate?

Question

It has been appeared $\delta(\mathbf r_1 -\mathbf r_2)$ in one of my integrals which I have to calculate. The differentials of the integral are $dr_1d\theta_1d\phi_1 dr_2d\theta_2d\phi_2$ So I need to write this Dirac delta in spherical coordinates. I just found $$ \delta(\mathbf r_1 -\mathbf r_2)=\frac{1}{r_1^2}\delta(r_1 -r_2)\delta(\cos\theta_1 -\cos\theta_2)\delta(\phi_1 -\phi_2) $$ by searching the web, but I'm not sure about that. How can I deal with this problem?

Addendum

This is my code and I try to calculate the result of ep. Note that by DirecDelta[re-rp] I mean $\delta(\mathbf r_e-\mathbf r_p)$ namely vectors of $r_e$ and $r_p$ beacaue I didn't know how to write them in my code (spirit of my question)!

mp = 1836;
\[Omega] = 10^-2;
me = 1;
M = mp + me;
{\[Alpha], \[Beta]} = \
{0.99930605221867707665950320006231777369976043701171875`50., 
   0.00014974969686706935865193324186606105286045931279659271240234375\
`50.};
\[Gamma] = 1/2*M*\[Omega];
norm = 2.12025232778742133718363166442892649864;
a = SetPrecision[2.35, 50];
b = SetPrecision[2.4, 50];
c = SetPrecision[3.2, 50];
intrE[r_, re_] := (
  E^(-2 r (\[Alpha] + r \[Beta])) (E^(-((
      2 (me re + mp (-r + re))^2 \[Gamma])/M^2)) - 
     E^(-((2 (me re + mp (r + re))^2 \[Gamma])/M^2))) M^2)/(
  4 mp (me + mp) \[Gamma]);
intE[re_?NumericQ] := (2*\[Pi])/re norm^2*
   NIntegrate[r intrE[r, re], {r, 0, \[Infinity]}];
intrp[r_, rp_] := (
  E^(-2 r (\[Alpha] + r \[Beta])) (E^(-((
      2 (mp rp + me (-r + rp))^2 \[Gamma])/M^2)) - 
     E^(-((2 (mp rp + me (r + rp))^2 \[Gamma])/M^2))) M^2)/(
  4 me (me + mp) \[Gamma]);
intP[rp_?NumericQ] := (2*\[Pi])/rp norm^2*
   NIntegrate[r intrp[r, rp], {r, 0, \[Infinity]}];



ep = -NIntegrate[(
   re^2*rp^2*intE[re]*intP[rp]*
    DiracDelta[re - rp] Sin[\[Theta]e] Sin[\[Theta]p])/(
   a - b*Sqrt[intE[re]*intP[rp]] + c *intE[re]*intP[rp]), {\[Theta]e, 
    0, \[Pi]}, {\[Theta]p, 0, \[Pi]}, {re, 0, \[Infinity]}, {rp, 
    0, \[Infinity]}, {\[Phi]e, 0, 2 \[Pi]}, {\[Phi]p, 0, 2 \[Pi]}]

Answer 1

The integration Integrate[\[Delta][...]dV about the infinite sphere must evaluate to 1!

Test of your formula doesn't evaluate to one, but for the modified approach

$\delta (r-\text{r0},\theta -\text{$\theta $0},\phi -\text{$\phi $0})=\frac{\csc (\theta ) \tan (\text{$\theta $0}) \delta (\theta -\text{$\theta $0}) \delta (r-\text{r0}) \delta (\phi -\text{$\phi $0})}{r^2} $

the integration gives

sol=Integrate[ DiracDelta[r - r0]/(r^2 Sin[\[Theta]])DiracDelta[\[Phi] - \[Phi]0] DiracDelta[ \[Theta] - \[Theta]0 ]/Cot[\[Theta]0]  r^2 Cos[\[Theta]] 
, {r, 0,Infinity}, {\[Theta], -Pi/2, Pi/2}, {\[Phi], 0, 2 Pi}, 
Assumptions -> {r0 > 0, 0 < \[Phi] < 2 Pi, -Pi/2 < \[Theta] < Pi/2,Element[{r0, \[Theta]0, \[Phi]0}, Reals]} ]

(*UnitStep[\[Pi] - 2 \[Theta]0, \[Pi]/2 + \[Theta]0, 2 \[Pi] - \[Phi]0, \[Phi]0]*)

Plot3D[sol, {\[Phi]0, 0, 2 Pi}, {\[Theta]0, -Pi/2, Pi/2},AxesLabel -> {\[Phi]0, \[Theta]0} ]

The result equals 1 for 0 < \[Phi]0 < 2 Pi, -Pi/2 < \[Theta]0 < Pi/2.

addendum/modified

The added integral ep==... contains DiracDelta["re" -"rp"] which stands for the vectorform of DiracDelta (see first part of the question)

It must be substituted by

dirac = DiracDelta[re - rp]/re^2 DiracDelta[\[Phi]e - \[Phi]p] DiracDelta[\[Theta]e - \\[Theta]p]/( Sin[\[Theta]e] Cot[\[Theta]p])

Knowing Integrate[...DiracDelta[x - y],{x,0,Infinity}, {y,0,Infinity}]= Integrate[... /.x->y, {y,0,Infinity}] the integration can be partially executed. The integral can be divided into three parts

ep = -NIntegrate[ (re^2*rp^2*intE[re]*intP[rp]*1/re^2)/(a - b*Sqrt[intE[re]*intP[rp]] + c*intE[re]*intP[rp])  /.re -> rp  , {rp, 0, \[Infinity]}]*
NIntegrate[1 , {\[Phi]p, 0, 2 \[Pi]}]*
NIntegrate[1/( Sin[\[Theta]e] Cot[\[Theta]p])Sin[\[Theta]e] Sin[\[Theta]p] /. \[Theta]e -> \[Theta]p, {\[Theta]p, - Pi/2, \[Pi]/2}]    
(*-2.09081 *)