Blog: Erini's Update 4/16/2013: LE_module.txt

Module LE_module

USE DataDeclarations
USE GlobalDeclarations
USE PhysicalDeclarations

Implicit None

Private

real xi, dtheta, poly_index, theta, xlast, dxdtlast,k1,k2,k3,k4,h,f,xi_last,dtheta_last,theta_last,dtheta1

PUBLIC CalcParams,LE_params,

CONTAINS

Subroutine CalcParams(R_star,M_star,poly_index,h,rho_cent,xi_last,dtheta1,K,a_scale,p_cent,rho_profile)

!requires a user input star/clump radius, a polytropic index, and a central density and integration length h 
!and calculated xi_last (root), and corresponding dtheta/dxi
!outputs scaling constant a, mass of star/clump,gas constant K, and central Pressue

Real(Kind=qPrec) :: R_star,M_star,rho_new,h,poly_index,rho_cent,xi_last,dtheta1,K,a_scale,rho_profile
Real:: xi_last

a_scale=xi_last/R_star

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

K=4*pi*G/((poly_index+1)(a_scale**2))*rho_cent**((poly_index-1)/poly_index)

p_cent=K*rho_cent**(1+1/poly_index)

rho_profile=rho_cent*theta**poly_index

end subroutine calcparams

!initial conditions
numsteps=0
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)

Function LE_params(xi,dtheta,poly_index,theta,xi_last,dtheta1,dtheta_last,theta_last)
Real:: LE_rho,xi,xi, dtheta, poly_index, theta,k1,k2,k3,k4,h,f,xi_last,dtheta_last,theta_last,dtheta1
do while(theta.ge.0)
xi=xi+h
theta=theta+h*dtheta
!calculating the runge kutta

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

theta_last=theta
dtheta_last=dtheta
enddo   

print*,'final value for xi',xi
print*,'final value for theta',theta
print*,'final value for dtheta/dphi',dtheta

!final conditions (for root)
xi = xi + (theta/(theta_last-theta))*h
dtheta1 = dtheta - (xi-xi_last)*(dtheta-dtheta_last)/h

        write(unit=6,FMT=*) 'first root of phi at xi=',xi_last 
        write(unit=6,FMT=*) 'at this point dtheta/dt = ',dtheta1
end     


!function f calculate the second order dtheta

real function f(poly_index,xi,theta,dtheta)
implicit none
real poly_index,xi,theta,dtheta
f=-theta**poly_index-2.0*dtheta/xi
return
end