BackLinksMenu()

LE_module.f90

Constructs a Polytrope Profile

Requirements

poly_index	This is the index,n, of the polytrope equation of state
rho_cent	The central density of the clump object (cgs)
M_star	The mass of the object (cgs)

If the polytrope index n that is given is greater than or equal to 5, we will get a solution with infinite radius - if the user explicitly provides a xi_last then the profile can be constructed, otherwise it is stopped.

Output

a_scale	scaling factor from dimensionless to physical values
p_cent	Central pressure of clump
t_cent	Central temp
R_outer	Radius of object
data(*,1:3)	profiles - data(*,2): density profile, data(*,3): pressure profile, data(*,4): temp profile

Code

Module LE_module
   USE GlobalDeclarations
   USE PhysicsDeclarations
   USE Profiles
   Implicit None

   Private

   PUBLIC ConstructPolyTropeProfile

CONTAINS

   SUBROUTINE ConstructPolyTropeProfile(Profile,R_outer, poly_index, rho_cent, M_star, xi_last,a_scale,p_cent,temp_cent,rho_outer)
      TYPE(ProfileDef), POINTER :: Profile
      INTEGER(KIND=8) :: nPoints
      REAL(KIND=QPREC) :: nPointsf
      REAL(KIND=QPREC) :: h, poly_index,press_profile,rho_profile
      REAL(KIND=QPREC) :: M_star,rho_cent  ! INPUT PARAMETERS
      REAL(KIND=QPREC) ::  xi, xi_last, dxi, dtheta, dtheta_last, theta, theta_last,ddxi,nsteps !calculated dimensionless variables
      REAL(KIND=QPREC) :: a_scale, R_outer, p_cent,temp_cent !calculated scaling values and other useful relations
      REAL(KIND=qPREC) :: theta_target, rho_outer
      INTEGER :: i,j
      !Profile data stored in Profile object 
      
      IF (poly_index >= 5) THEN
         IF (MPI_ID == 0) THEN
            write(*,*) 'not yet fully implemented for poly_index >= 5'
            STOP
         END IF
         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
         
         theta_target=(rho_outer/rho_cent)**(1d0/poly_index)

         !Initilize solution at some small xi
         xi=1d-4
         theta=1d0-(xi**2d0)/6d0+(poly_index/120d0)*(xi**4d0)-(8d0*poly_index**2d0-5d0*poly_index)/(15120d0)*(xi**5d0)
         dtheta=-xi/3d0+(poly_index/30d0)*(xi**3d0)-(8d0*poly_index**2d0-5d0*poly_index)/(2520d0)*(xi**5d0)

         !Initial step size
         dxi=1d-2
         DO
            CALL LE_Advance(xi,dtheta,poly_index,theta,dxi,1,0d0)
            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
!            write(*,*) xi, theta
         END DO
         xi_last=xi
         theta_last=theta
         dtheta_last=dtheta
         write(*,*) "xi_last :" ,xi_last   
         write(*,*) "dtheta_last: ", dtheta_last
       END IF

       !Calculate length scale
       a_scale=(-M_star/(4d0*Pi*rho_cent*(xi_last**2)*dtheta_last))**(1d0/3d0)
       write(*,*) "a_scale", a_scale, M_star, rho_cent, xi_last, dtheta_last

       !Calculate polytrope physical radius
       R_outer=a_scale*xi_last

       !Determine number of points based on resolution with a minimum of 1000
       nPoints=max(ceiling(R_outer/(levels(maxLevel)%dx*lScale)),999)+1

       !Physical step size between profile points
       h=R_outer/real(nPoints-1) 
       write(*,*) "step size h: ",h

      !Construct Profile Object and initialize default values for no polytrope
      CALL CreateProfile(Profile, nPoints, (/Mass_Field, P_Field, Temp_Field/), RADIAL)          
      Profile%data(:,1)=h*(/(i-1,i=1,nPoints)/)
      Profile%data(:,2:4)=0d0

      !Calculate central pressure
      p_cent=4d0*Pi*G*(a_scale**2)*(rho_cent**2)/(poly_index+1d0)
      write(*,*) "p_cent: ", p_cent

      !Initialize center cell
      profile%data(1,1:3)=(/0d0, rho_cent, p_cent/)
     

      !Start solution at small positive xi - but don't go past 1st point
      xi=min(h, 1d-4)
      theta=1d0-(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)
      dxi=h/a_scale-xi

      !Advance to second profile point
      write(*,*), "dxi:  ",dxi
      CALL LE_Advance(xi,dtheta,poly_index,theta,dxi,100,theta_target)
      write(*,*),"xi...WORK",xi  
      profile%data(2,1:3)=(/xi*a_scale,rho_cent*theta**poly_index, p_cent*theta**(poly_index+1d0)/)
     
      !Calculate step sizes for rest of profile
      dxi=h/a_scale
      write(*,*) "dxi for rest of profile: ", dxi     

      !Initialize profile
      DO i=2, size(Profile%data, 1)
         CALL LE_ADvance(xi,dtheta,poly_index,theta,dxi,100,theta_target)
         Profile%data(i,1:3)=(/xi*a_scale,rho_cent*theta**poly_index,p_cent*theta**(poly_index+1d0)/)
!         Print *, "Density: ", Profile%data(i,2)
!         IF (theta <= theta_target) EXIT
      END DO
      theta=theta_target
      Print *, xi
      profile%data(:,4)=profile%data(:,3)/profile%data(:,2) * Xmu*Hmass/Boltzmann
      
      IF (MPI_ID == 0) THEN
         print*,'final value for xi',xi
         print*,'final value for theta',theta
         print*,'final value for dtheta/dphi',dtheta
      END IF


   END SUBROUTINE ConstructPolyTropeProfile

   SUBROUTINE LE_Advance(xi,dtheta,poly_index,theta,step,nsteps,theta_target)
      REAL(KIND=QPREC):: xi, dtheta, poly_index, theta, h, step, theta_target
      REAL(KIND=QPREC) :: k1(2), k2(2),k3(2),k4(2)
      REAL(KIND=QPREC), Dimension(2) :: g
      REAL(KIND=QPREC) :: stepsize
      INTEGER :: i,nsteps
      IF (theta <= theta_target) THEN !use zero value for outside of polytrope
         xi=xi+step
         theta=theta_target
         RETURN
      END IF
      h=step/real(nsteps)
      DO i=1,nsteps
         g=(/theta, dtheta/)
         k1=h*f(poly_index,xi,g)
         k2=h*f(poly_index,xi+h/2.0,g+k1/2.0)
         k3=h*f(poly_index,xi+h/2.0,g+k2/2.0)
         k4=h*f(poly_index,xi+h,g+k3)

         g=g+(k1+2.0*k2+2.0*k3+k4)/6.0

         IF (g(1) <= theta_target .OR. ISNAN(g(1))) THEN !use previous values to estimate root
            xi = xi - ((theta-theta_target)/dtheta)
            theta=theta_target
            EXIT
         END IF
         theta=g(1)
         dtheta=g(2)
         xi=xi+h
      END DO
   end SUBROUTINE LE_Advance


   !function f calculate the second order dtheta

   function f(poly_index,xi,g)
      implicit none
      REAL(KIND=qPREC) poly_index,xi, f(2)
      REAL(KIND=qPREC), dimension(2) :: g
      f(1)=g(2)
      f(2)=-g(1)**poly_index-2.0*g(2)/xi
   end function f
END MODULE LE_module }}}