the general equation pairs are:

i[R_]=Simplify[2 Integrate[rho[r] r /Sqrt[r^2-R^2],{r,R,Infinity}],R>0]
rho[r_]=-1/Pi Integrate[D[i[R],R]/Sqrt[R^2-r^2],{R,r,Infinity}]

rho[x_]=Exp[-x^2]
i[R_]=Simplify[2 Integrate[rho[r] r /Sqrt[r^2-R^2],{r,R,Infinity}],R>0]

In[4]:= Sqrt[Pi]/E^R^2

-1/Pi Integrate[D[i[R],R]/Sqrt[R^2-r^2],{R,r,Infinity}]

                                 2
                               -r
Out[5]= ConditionalExpression[E   , Re[r] > 0 && Im[r] == 0]

since r is >0, this is the same gaussian back

---r=1, rho=1 sphere--

i[R_]=2 Integrate[ r/Sqrt[r^2-R^2],{r,R,1},GenerateConditions->False]

                     2
Out[10]= 2 Sqrt[1 - R ]


-1/Pi Integrate[D[i[R],R]/Sqrt[R^2-r^2],{R,r,1}]


Out[11]= ConditionalExpression[1, 0 < Re[r] < 1 && Im[r] == 0]
yep; same as started (r is real and when r<1 we are inside the ball and rho=1)

---starting shell---

shell[R_]=If[R<1,i[R],0]
inside[R_]=If[R<.9,Evaluate[-2 Integrate[ r /Sqrt[r^2-R^2],{r,R,.9},GenerateConditions->False]],0]
                                     2
Out[13]= If[R < 0.9, -2 Sqrt[0.81 - R ], 0]

total[R_]=If[R<.9, inside[R]+shell[R],shell[R]]

Plot[total[R],{R,0,1.1},PlotRange->All]  

---spiral disk---
rho[r_]=Exp[-r]
i[R_]=Simplify[2 Integrate[rho[r] r /Sqrt[r^2-R^2],{r,R,Infinity}],R>0]
Out[19]= 2 R BesselK[1, R]
Plot[{i[R],rho[R]},{R,0,3},PlotRange->All]

---elliptical galaxy---
i[R_]=Exp[-7.669(R^(1/4)-1)]
Plot[{i[x]},{x,.5,4},AspectRatio->1] 
ParametricPlot[{x,2.5 Log10[i[x^4]]},{x,0,Sqrt[2]},AspectRatio->1] 

rho[r_]:=-1/Pi NIntegrate[D[i[R],R]/Sqrt[R^2-r^2],{R,r,16}]
Plot[{rho[x],i[x]},{x,.5,4},AspectRatio->1] 

ParametricPlot[{x,2.5 Log10[rho[x^4]]},{x,0,Sqrt[2]},AspectRatio->1]