! 20111212 Alexander Heger ! 1/z**2 d/dz (z**2 d/dz theta(z)) + theta(z)**n = 0 ! for small z one can approximate ! theta(z) = 1. + (-1/6.)*z**2 + (n/120.)*z**4 + O(z**6) ! Therefore lim(z)-->0 d**2 theta(z)/d z**2 = -1/3 ! 20120423 Alexander Heger ! if we include a constant rotation rate Omega the equation becomes ! 1/z**2 d/dz (z**2 d/dz theta(z)) + theta(z)**n - w = 0 ! where ! w = W/rho_c ! W = 2 Omega**2 / 4 pi G ! for small z one can approximate ! theta(z) = 1. + (w - 1)/6. *z**2 + (1.-w)*(n/120.)*z**4 + O(z**6) ! Therefore lim(z)-->0 d**2 theta(z)/d z**2 = (w-1.)/3 module types implicit none INTEGER, PARAMETER :: int32 = SELECTED_INT_KIND(8) INTEGER, PARAMETER :: int64 = SELECTED_INT_KIND(16) INTEGER, PARAMETER :: real64 = SELECTED_REAL_KIND(15) end module types module lanedata use types, only: & real64, int32 implicit none save integer(kind=int32), parameter :: & maxdata = 2**20-1 real(kind=real64), dimension(:,:), allocatable :: & theta integer(kind=int32) :: & ndata end module lanedata subroutine freelanedata() use lanedata, only: & ndata, theta implicit none if (allocated(theta)) deallocate(theta) ndata = -1 end subroutine freelanedata subroutine lane(dx, n, w) ! lane emden intrgration, return data, ! including first invalid (rho < 0) point for interpolation. use types, only: & int32, real64 use lanedata, only: & maxdata, theta, ndata implicit none !f2py real(kind=real64), intent(in) :: dx, n, w real(kind=real64), intent(in) :: & dx, n, w real(kind=real64) :: & x integer(kind=int32) :: & i, j real(kind=real64), dimension(0:1) :: & y, z real(kind=real64), dimension(:,:), allocatable :: & theta_ if (allocated(theta)) deallocate(theta) j = maxdata allocate(theta_(0:j, 0:1)) x = 0.d0 y(0:1) = (/1.d0, 0.d0/) i = 0 theta_(i,:) = y(:) do while (.True.) call rk4(x,y(0),y(1),dx,n,w,z(0),z(1)) if (i == j) then j = j + maxdata allocate(theta(0:j,0:1)) theta(0:i,:) = theta_(0:i,:) deallocate(theta_) call move_alloc(theta, theta_) endif i = i+1 theta_(i, :) = z(:) x = x + dx y(:) = z(:) if (y(0) < 0.d0) exit enddo ndata = i allocate(theta(0:ndata, 0:1)) theta(0:ndata,0:1) = theta_(0:ndata,0:1) deallocate(theta_) end subroutine lane subroutine rk4(x0, y0, y1, dx, n, w, z0, z1) use types, only: & real64 implicit none real(kind=real64), parameter :: & p13 = 1.d0 / 3.d0, & p16 = 1.d0 / 6.d0 real(kind=real64), intent(in) :: & x0, y0, y1 real(kind=real64), intent(in) :: & dx, n, w real(kind=real64), intent(out) :: & z0, z1 real(kind=real64) :: & xh, dh real(kind=real64) :: & k10, k11, k20, k21, k30, k31, k40, k41 !f2py real(8), intent(in) :: x0, dx, n, w !f2py real(8), intent(in) :: y0, y1 !f2py real(8), intent(out) :: z0, z1 xh = x0 + 0.5d0 * dx dh = 0.5d0 * dx k10 = y1 if (x0 == 0) then k11 = (w - 1.d0) * p13 else k11 = -2.d0 / x0 * y1 - (max(y0, 0.d0))**n + w endif k20 = y1 + dh*k11 k21 = -2.d0 / xh * k20 - (max(y0 + dh*k10, 0.d0))**n k30 = y1 + dh*k21 k31 = -2.d0 / xh * k30 - (max(y0 + dh*k20, 0.d0))**n k40 = y1 + dx*k31 k41 = -2.d0 / (x0+dx) * k40 - (max(y0 + dx*k30,0.d0))**n z0 = y0 + dx*(k10 + 2.d0 * (k20 + k30) + k40) * p16 z1 = y1 + dx*(k11 + 2.d0 * (k21 + k31) + k41) * p16 end subroutine rk4