| 1 | Module LE_module
|
|---|
| 2 |
|
|---|
| 3 | USE GlobalDeclarations
|
|---|
| 4 | USE PhysicalDeclarations
|
|---|
| 5 | USE Profiles
|
|---|
| 6 | Implicit None
|
|---|
| 7 |
|
|---|
| 8 | Private
|
|---|
| 9 |
|
|---|
| 10 | PUBLIC ConstructPolyTropeProfile
|
|---|
| 11 |
|
|---|
| 12 | CONTAINS
|
|---|
| 13 |
|
|---|
| 14 | SUBROUTINE ConstructPolyTropeProfile(Profile, Mass_Profile,P_profile,rho_profile,R_outer, poly_index, rho_cent, M_star, xi_last,a_scale,p_cent)
|
|---|
| 15 | TYPE(ProfileDef), POINTER :: Profile
|
|---|
| 16 | INTEGER :: nPoints
|
|---|
| 17 | REAL :: R_outer, poly_index, rho_cent, xi_last
|
|---|
| 18 | REAL :: h, a_scale, K, p_cent, dxi, xi, theta, dtheta, M_star
|
|---|
| 19 | PARAMETER(pi=3.14159265359)
|
|---|
| 20 |
|
|---|
| 21 | a_scale=(-M_star/(4*pi*rho_cent*(xi_last**2)*dtheta_last))**(1/3)
|
|---|
| 22 | R_outer=a_scale*xi_last
|
|---|
| 23 | nPoints=max(R_outer/levels(maxLevel)%dx, 1000)
|
|---|
| 24 | h=R_outer/nPoints !Computational step size
|
|---|
| 25 |
|
|---|
| 26 | IF (poly_index >= 5) THEN
|
|---|
| 27 | IF (xi_last <= 0) THEN
|
|---|
| 28 | IF (MPI_ID == 0) write(*,*) 'Error - must specify dimensionless truncation cutoff xi_last if the polytropic index > 5 in ConstructPolyTropeProfile'
|
|---|
| 29 | STOP
|
|---|
| 30 | END IF
|
|---|
| 31 | ELSE
|
|---|
| 32 | !Initilize solution at some small xi
|
|---|
| 33 | xi=1e-4
|
|---|
| 34 | 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)
|
|---|
| 35 | 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)
|
|---|
| 36 |
|
|---|
| 37 | !Initial step size
|
|---|
| 38 | dxi=1e-2
|
|---|
| 39 |
|
|---|
| 40 | DO
|
|---|
| 41 | theta_last=theta
|
|---|
| 42 | xi_last=xi
|
|---|
| 43 | dtheta_last=dtheta
|
|---|
| 44 | CALL LE_Advance(xi, theta, dtheta, dxi)
|
|---|
| 45 | IF (xi > 10d0) dxi=.01*theta/dtheta !adjust step size once were beyond xi=10 in case poly_index is just under 5...
|
|---|
| 46 | IF (theta < 0d0) EXIT
|
|---|
| 47 | END DO
|
|---|
| 48 | dtheta = dtheta - (xi-xi_last)*(dtheta-dtheta_last)/dxi
|
|---|
| 49 | xi = xi + (theta/(theta_last-theta))*dxi
|
|---|
| 50 | theta=0d0
|
|---|
| 51 | xi_last=xi
|
|---|
| 52 | theta_last=theta
|
|---|
| 53 | dtheta_last=dtheta
|
|---|
| 54 | END IF
|
|---|
| 55 |
|
|---|
| 56 | !Construct Profile Object
|
|---|
| 57 | CALL CreateProfile(Profile, nPoints, (/Mass_Field, Press_Field/), RADIAL)
|
|---|
| 58 |
|
|---|
| 59 |
|
|---|
| 60 | !Calculate various other scales
|
|---|
| 61 | a_scale=(-M_star/(4*pi*rho_cent*(xi_last**2)*dtheta_last))**(1/3)
|
|---|
| 62 | p_cent=4*pi*ScaleGrav*(a_scale**2)*(rho_cent**2)*(theta_last**(poly_index+1d0))/(poly_index+1)
|
|---|
| 63 |
|
|---|
| 64 | !Start solution at small positive xi
|
|---|
| 65 | xi=1e-4
|
|---|
| 66 | 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)
|
|---|
| 67 | 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)
|
|---|
| 68 |
|
|---|
| 69 | !Integrate to first profile zone
|
|---|
| 70 | dxi=.5*(h/a_scale)-xi
|
|---|
| 71 | nsteps=ceiling(dxi/.01)
|
|---|
| 72 | ddxi=dxi/nsteps
|
|---|
| 73 | DO j=1,nsteps
|
|---|
| 74 | CALL LE_Advance(xi, theta, dtheta, ddxi)
|
|---|
| 75 | END DO
|
|---|
| 76 |
|
|---|
| 77 |
|
|---|
| 78 | !Calculate step sizes for rest of profile
|
|---|
| 79 | dxi=h/a_scale
|
|---|
| 80 | nsteps=ceiling(h/a_scale/dxi)
|
|---|
| 81 | ddxi=dxi/nsteps
|
|---|
| 82 |
|
|---|
| 83 | data(:,1)=h*(/(i,i=1,npoints)/)-.5*h
|
|---|
| 84 | data(:,2:3:4)=0
|
|---|
| 85 |
|
|---|
| 86 | !Initialize profile
|
|---|
| 87 | DO i=1, size(data, 1)
|
|---|
| 88 | ! Convert from dimenionless to computational units
|
|---|
| 89 | data(i,2)=rho_cent*theta**poly_index
|
|---|
| 90 | data(i,3)=4*pi*ScaleGrav*(a_scale**2)*(rho_cent**2)*(theta**(poly_index+1d0))/(poly_index+1)
|
|---|
| 91 | data(i,4)=-4*pi*(a_scale**3)*(xi**2)*rho_cent*dtheta/(Msun)
|
|---|
| 92 | DO j=1,nsteps
|
|---|
| 93 | !backup values for theta, xi, dtheta in case we overshoot
|
|---|
| 94 | theta_last=theta
|
|---|
| 95 | xi_last=xi
|
|---|
| 96 | dtheta_last=dtheta
|
|---|
| 97 |
|
|---|
| 98 | !Advance to next position
|
|---|
| 99 | CALL LE_Advance(xi, theta, dtheta, ddxi)
|
|---|
| 100 | IF (theta < 0d0) EXIT
|
|---|
| 101 | END DO
|
|---|
| 102 | IF (theta < 0d0) EXIT
|
|---|
| 103 | END DO
|
|---|
| 104 |
|
|---|
| 105 | !
|
|---|
| 106 | dtheta = dtheta - (xi-xi_last)*(dtheta-dtheta_last)/ddxi
|
|---|
| 107 | xi = xi + (theta/(theta_last-theta))*ddxi
|
|---|
| 108 | theta=0d0
|
|---|
| 109 |
|
|---|
| 110 | M_star=4d0*pi*rho_cent*(R_star**3)*(-1d0/xi_last)*(dtheta)
|
|---|
| 111 |
|
|---|
| 112 | IF (MPI_ID == 0) THEN
|
|---|
| 113 | print*,'final value for xi',xi
|
|---|
| 114 | print*,'final value for theta',theta
|
|---|
| 115 | print*,'final value for dtheta/dphi',dtheta
|
|---|
| 116 |
|
|---|
| 117 | print*, 'first root of phi at xi=',xi
|
|---|
| 118 | print*, 'at this point dtheta/dt = ',dtheta
|
|---|
| 119 | END IF
|
|---|
| 120 | END SUBROUTINE ConstructPolyTropeProfile
|
|---|
| 121 |
|
|---|
| 122 | SUBROUTINE LE_Advance(xi,dtheta,poly_index,theta,h)
|
|---|
| 123 | Real:: xi, dtheta, poly_index, theta, h
|
|---|
| 124 | REAL :: k1(2), k2(2),k3(2),k4(2),g(2)
|
|---|
| 125 |
|
|---|
| 126 | g=(/theta, dtheta/)
|
|---|
| 127 | k1=h*f(poly_index,xi,g)
|
|---|
| 128 | k2=h*f(poly_index,xi+h/2.0,g+k1/2d0)
|
|---|
| 129 | k3=h*f(poly_index,xi+h/2.0,g+k2/2d0)
|
|---|
| 130 | k4=h*f(poly_index,xi+h,g+k3)
|
|---|
| 131 |
|
|---|
| 132 | g=g+(k1+2.0*k2+2.0*k3+k4)/6.0
|
|---|
| 133 | theta=g(1)
|
|---|
| 134 | dtheta=g(2)
|
|---|
| 135 | xi=xi+h
|
|---|
| 136 | end SUBROUTINE LE_Advance
|
|---|
| 137 |
|
|---|
| 138 |
|
|---|
| 139 | !function f calculate the second order dtheta
|
|---|
| 140 |
|
|---|
| 141 | function f(poly_index,xi,g)
|
|---|
| 142 | implicit none
|
|---|
| 143 | real poly_index,xi,g(2), f(2)
|
|---|
| 144 | f(1)=g(2)
|
|---|
| 145 | f(2)=-g(1)**poly_index-2.0*g(2)/xi
|
|---|
| 146 | end function f
|
|---|
