! 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