Self-gravitating equilibrium profile using density profile and external pressure

Mathematically, we are solving the equation of hydrostatic equilibrium

\nabla P(r) = -\rho(r) \nabla \phi(r)

subject to the constraint that P(R)=P_0.

Since \phi(r) = -\frac{G M(r)}{r} we can rewrite the equation for HSE as

begin{equation} \frac{dP}{dr}=-\frac{GM(r)\rho(r)}{r^2} \end{equation}

and we can express the enclosed mass as

setcounter{equation}{1} \begin{equation} M(r)=4\pi\displaystyle\int_0^r{\rho(r') r'^2dr'} \end{equation}

Now we can discretize equation 2 as

setcounter{equation}{2} \begin{equation} M_{i+1}=M_i + 4\pi \rho_i r_i^2\Delta r_i \end{equation}

where \Delta r_j = r_{j+1}-r_j

and equation 1 as

setcounter{equation}{3} \begin{equation} P_{i+1}=P_i-\Delta r_i \frac{G M_i \rho_i}{r_i^2} \end{equation}

but a taylor series expansion shows that the error at each step in the integration goes as

\frac{\Delta r_j}{r_j} + \frac{\Delta \rho_j}{\rho_j} + ...

so we need to limit out point spacing - or improve the accuracy of our discretization.

If we apply linear interpolation to our density profile, we can discretize equation 2 as

setcounter{equation}{5} \begin{equation} M_i=M(r_i)=4\pi\displaystyle\sum_{j=0}^i{\rho_j r_j^2\Delta r_j \left ( 1 + \frac{\Delta r_j}{r_j} + \frac{\Delta \rho}{2\rho} + \frac{\Delta r_j^2}{3 r_j^2} + \frac{2\Delta \rho \Delta r_j}{3 \rho r_j} + \frac{\Delta \rho \Delta r_j^2}{4 \rho r_j^2}\right )} \end{equation}

and equation 4 becomes a bit of a monster - and we have to worry about limits as r → 0.

if (r0 == 0d0) then
    dp=drho*((pi*drho*dr**4)/4d0 + (4d0*pi*rho*dr**3)/9d0) + (rho*(pi*drho*dr**3 + 2d0*pi*rho*dr**2))/3d0
else
    epsilon=log(r0/(r0+dr))
    dp=((-36d0*drho*r0)*epsilon*M0*dr + (-36d0*drho*r0**2)*epsilon*M0 + (- 12d0*pi*drho**2*r0**5 + 48d0*pi*rho*drho*r0**4)*epsilon*dr + (-12d0*pi*drho**2*r0**6 + 48d0*pi*rho*drho*r0**5)*epsilon + (36d0*rho - 36d0*drho*r0)*M0*dr + (9d0*pi*drho**2*r0)*dr**5 + (17d0*pi*drho**2*r0**2 + 28d0*pi*rho*drho*r0)*dr**4 + (2d0*pi*drho**2*r0**3 + 64d0*pi*drho*r0**2*rho + 24d0*pi*r0*rho**2)*dr**3 + (- 6d0*pi*drho**2*r0**4 + 24d0*pi*drho*r0**3*rho + 72d0*pi*r0**2*rho**2)*dr**2 +(- 12d0*pi*drho**2*r0**5 + 48d0*pi*rho*drho*r0**4)*dr)/(36d0*r0**2 + 36d0*dr*r0)
end if
data(i,i_Press)=data(i+1,i_Press)+ScaleGrav*dp