Module LE_module

   USE GlobalDeclarations
   USE PhysicalDeclarations
   USE Profiles
   Implicit None

   Private

   PUBLIC ConstructPolyTropeProfile

CONTAINS

   SUBROUTINE ConstructPolyTropeProfile(Profile, Mass_Profile,P_profile,rho_profile,R_outer, poly_index, rho_cent, M_star, xi_last,a_scale,p_cent)
      TYPE(ProfileDef), POINTER :: Profile
      INTEGER :: nPoints
      REAL :: R_outer, poly_index, rho_cent, xi_last
      REAL :: h, a_scale, K, p_cent, dxi, xi, theta, dtheta, M_star
      PARAMETER(pi=3.14159265359)

      a_scale=(-M_star/(4*pi*rho_cent*(xi_last**2)*dtheta_last))**(1/3)
      R_outer=a_scale*xi_last
      nPoints=max(R_outer/levels(maxLevel)%dx, 1000)
      h=R_outer/nPoints !Computational step size

      IF (poly_index >= 5) THEN
         IF (xi_last <= 0) THEN
            IF (MPI_ID == 0) write(*,*) 'Error - must specify dimensionless truncation cutoff xi_last if the polytropic index > 5 in ConstructPolyTropeProfile'
            STOP
         END IF
      ELSE
         !Initilize solution at some small xi
         xi=1e-4
         theta=1-(xi**2.0)/6.0+(poly_index/120.0)*(xi**4.0)-(8.0*poly_index**2.0-5.0*poly_index)/(15120.0)*(xi**5.0)
         dtheta=-xi/3.0+(poly_index/30.0)*(xi**3.0)-(8.0*poly_index**2.0-5.0*poly_index)/(2520.0)*(xi**5.0)

         !Initial step size
         dxi=1e-2

         DO
            theta_last=theta
            xi_last=xi
            dtheta_last=dtheta
            CALL LE_Advance(xi, theta, dtheta, dxi)
            IF (xi > 10d0) dxi=.01*theta/dtheta !adjust step size once were beyond xi=10 in case poly_index is just under 5...
            IF (theta < 0d0) EXIT
         END DO
         dtheta = dtheta - (xi-xi_last)*(dtheta-dtheta_last)/dxi
         xi = xi + (theta/(theta_last-theta))*dxi
         theta=0d0
         xi_last=xi
         theta_last=theta
         dtheta_last=dtheta
      END IF

      !Construct Profile Object
      CALL CreateProfile(Profile, nPoints, (/Mass_Field, Press_Field/), RADIAL)
      

      !Calculate various other scales
      a_scale=(-M_star/(4*pi*rho_cent*(xi_last**2)*dtheta_last))**(1/3)
      p_cent=4*pi*ScaleGrav*(a_scale**2)*(rho_cent**2)*(theta_last**(poly_index+1d0))/(poly_index+1)

      !Start solution at small positive xi
      xi=1e-4
      theta=1-(xi**2.0)/6.0+(poly_index/120.0)*(xi**4.0)-(8.0*poly_index**2.0-5.0*poly_index)/(15120.0)*(xi**5.0)
      dtheta=-xi/3.0+(poly_index/30.0)*(xi**3.0)-(8.0*poly_index**2.0-5.0*poly_index)/(2520.0)*(xi**5.0)
      
      !Integrate to first profile zone
      dxi=.5*(h/a_scale)-xi
      nsteps=ceiling(dxi/.01)
      ddxi=dxi/nsteps
      DO j=1,nsteps
         CALL LE_Advance(xi, theta, dtheta, ddxi)
      END DO

     
      !Calculate step sizes for rest of profile
      dxi=h/a_scale
      nsteps=ceiling(h/a_scale/dxi)
      ddxi=dxi/nsteps

      data(:,1)=h*(/(i,i=1,npoints)/)-.5*h
      data(:,2:3:4)=0

      !Initialize profile
      DO i=1, size(data, 1)
         ! Convert from dimenionless to computational units
         data(i,2)=rho_cent*theta**poly_index
         data(i,3)=4*pi*ScaleGrav*(a_scale**2)*(rho_cent**2)*(theta**(poly_index+1d0))/(poly_index+1)
         data(i,4)=-4*pi*(a_scale**3)*(xi**2)*rho_cent*dtheta/(Msun)
         DO j=1,nsteps
            !backup values for theta, xi, dtheta in case we overshoot
            theta_last=theta
            xi_last=xi
            dtheta_last=dtheta
           
            !Advance to next position
            CALL LE_Advance(xi, theta, dtheta, ddxi)
            IF (theta < 0d0) EXIT
         END DO
         IF (theta < 0d0) EXIT
      END DO

      !
      dtheta = dtheta - (xi-xi_last)*(dtheta-dtheta_last)/ddxi
      xi = xi + (theta/(theta_last-theta))*ddxi
      theta=0d0

      M_star=4d0*pi*rho_cent*(R_star**3)*(-1d0/xi_last)*(dtheta)

      IF (MPI_ID == 0) THEN
         print*,'final value for xi',xi
         print*,'final value for theta',theta
         print*,'final value for dtheta/dphi',dtheta
         
         print*, 'first root of phi at xi=',xi
         print*, 'at this point dtheta/dt = ',dtheta
      END IF
   END SUBROUTINE ConstructPolyTropeProfile

   SUBROUTINE LE_Advance(xi,dtheta,poly_index,theta,h)
      Real:: xi, dtheta, poly_index, theta, h
      REAL :: k1(2), k2(2),k3(2),k4(2),g(2)

      g=(/theta, dtheta/)
      k1=h*f(poly_index,xi,g)
      k2=h*f(poly_index,xi+h/2.0,g+k1/2d0)
      k3=h*f(poly_index,xi+h/2.0,g+k2/2d0)
      k4=h*f(poly_index,xi+h,g+k3)

      g=g+(k1+2.0*k2+2.0*k3+k4)/6.0
      theta=g(1)
      dtheta=g(2)
      xi=xi+h
   end SUBROUTINE LE_Advance


   !function f calculate the second order dtheta

   function f(poly_index,xi,g)
      implicit none
      real poly_index,xi,g(2), f(2)
      f(1)=g(2)
      f(2)=-g(1)**poly_index-2.0*g(2)/xi
   end function f