Changes between Version 12 and Version 13 of u/johannjc/scratchpad

Timestamp:

01/09/14 10:48:34 (12 years ago)

Author:

Jonathan

Comment:

—

u/johannjc/scratchpad

@@ -1 +1 @@
 == Self-gravitating equilibrium profile using density profile and external pressure ==
 
-Mathematically, we are solving the equation of hydrostatic equilibrium [[latex($\nabla P(r) = -\rho(r) \nabla \phi(r)$)]] subject to the constraint that  [[latex($P(R)=P_0$)]].  Since [[latex($\phi(r) = \frac{G M(r)}{r}$)]] we can rewrite the equation for HSE as [[latex($\frac{dP}{dr}=\frac{GM(r)\rho(r)}{r^2}$)]] and we can express the enclosed mass as [[latex($M(r)=4\pi\displaystyle\int_0^r{\rho(r') r'^2dr'}$)]] 
+Mathematically, we are solving the equation of hydrostatic equilibrium 
+
+[[latex($\nabla P(r) = -\rho(r) \nabla \phi(r)$)]] 
+
+subject to the constraint that  [[latex($P(R)=P_0$)]].  
+
+Since [[latex($\phi(r) = -\frac{G M(r)}{r}$)]] we can rewrite the equation for HSE as 
+{{{
+#!latex
+\begin{equation}
+\frac{dP}{dr}=-\frac{GM(r)\rho(r)}{r^2}
+\end{equation}
+}}}
 
 
+ and we can express the enclosed mass as 
+
+{{{
+#!latex
+\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
+
+{{{
+#!latex
+\setcounter{equation}{2}
+\begin{equation}
+M_{i+1}=M_i + 4\pi \rho_i r_i^2\Delta r_i
+\end{equation}
+}}}
+
+where [[latex($\Delta r_j = r_{j+1}-r_j$)]]
+
+and equation 1 as
+
+{{{
+#!latex
+\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 
+
+[[latex($\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
+
+{{{
+#!latex
+\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...
+{{{
+ 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)*epsi\
+lon + (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
+
+
+         if (r0 > 0d0) then
+            dp=dp
+         end if
+         data(i,i_Press)=data(i+1,i_Press)+ScaleGrav*dp
+}}}
+
+