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

Timestamp:

01/09/14 14:48:30 (11 years ago)

Author:

Jonathan

Comment:

—

u/johannjc/scratchpad

@@ -65 +65 @@
 }}}
 
-and equation 4 becomes a bit of a monster...
+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)*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
+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
 }}}