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 denisty of the clump object
xi_last	The root of the solve lame-emden equation - this represents the polytrope radius
R_outer	The maximum domain size, needed to determine computational step size

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

K	K is the free parameter from the polytropic equation of state
p_cent	Central pressure of clump
M_star	Mass of clump
data(*,1:3)	profiles - data(*,2): density profile, data(*,3): pressure profile

Code

Module LE_module

USE GlobalDeclarations USE PhysicalDeclarations USE Profiles Implicit None

Private

PUBLIC ConstructPolyTropeProfile

CONTAINS

SUBROUTINE ConstructPolyTropeProfile(Profile, R_outer, poly_index, rho_cent, M_star, xi_last)

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 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-(xi2.0)/6.0+(poly_index/120.0)*(xi4.0)-(8.0*poly_index2.0-5.0*poly_index)/(15120.0)*(xi5.0) dtheta=-xi/3.0+(poly_index/30.0)*(xi3.0)-(8.0*poly_index2.0-5.0*poly_index)/(2520.0)*(xi5.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

END IF

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

!Calculate various other scales a_scale=xi_last/R_outer K=4d0*pi*ScaleGrav/((poly_index+1)*a_scale2)*rho_cent((poly_index-1)/poly_index) p_cent=K*rho_cent(1d0+1d0/poly_index)

!Start solution at small positive xi xi=1e-4 theta=1-(xi2.0)/6.0+(poly_index/120.0)*(xi4.0)-(8.0*poly_index2.0-5.0*poly_index)/(15120.0)*(xi5.0) dtheta=-xi/3.0+(poly_index/30.0)*(xi3.0)-(8.0*poly_index2.0-5.0*poly_index)/(2520.0)*(xi5.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)=0

!Initialize profile DO i=1, size(data, 1)

! Convert from dimenionless to computational units data(i,2)=rho_cent*thetapoly_index data(i,3)=k*data(i,2)(1d0+1d0/poly_index)

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_star3)*(-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