问fortran (中子星振荡)打靶法

EN

   PROGRAM ROSCILLATION2
USE eos_parameters
IMPLICIT NONE
INTEGER ::i, j, k, l
INTEGER, PARAMETER :: N_ode = 5
REAL, DIMENSION(N_ode) :: y
REAL(8) :: rho0_cgs, rho0, P0, r0, phi0, pi
REAL(8) :: r, rend, mass, P, phi, delta, xi, eta
REAL(8) :: step, omega, omegastep, tiny, rho_print, Radius, B, a2, s0, lamda, E0, E
EXTERNAL :: fcn

!!!! User input 

rho0_cgs = 2.D+15           !central density in cgs unit
step = 1.D-4                ! step size dr
omegastep = 1.D-2           ! step size d(omega)
tiny = 1.D-8                ! small number P(R)/P(0) to define star surface
!!!!!!!!!
open(unit=15, file="data.dat", status="new")
pi = ACOS(-1.D0)

a2 =((((1.6022D-13)**4)*(6.674D-11)*((2.997D8)**-7)*((1.0546D-34)**-3)*(1.D6))**(0.5D0))*a2_MeV !convert to code unit (km^-1)

B = ((1.6022D-13)**4)*(6.674D-11)*((2.997D8)**-7)*((1.0546D-34)**-3)*(1.D6)*B_MeV !convert to code unit (km^-2)

s0 = (1.D0/3.D0) - (1/(6*pi**2))*a2*((1/(16*pi**2)*a2**2 + (pi**-2)*a4*(rho0 - B))**-0.5) !square of the spped of sound at r=0

lamda = -0.5D0*log(1-2*y(1)/r)

E0 = (r0**-2)*s0*exp(lamda + 3*phi0)



rho0 = rho0_cgs*6.67D-18 / 9.D0  !convert rho0 to code unit (km^-2)

!! Calculate central pressure  P0

P0 = (1.D0/3.D0)*rho0 - (4.D0/3.D0)*B - (1.D0/(a4*(12.D0)*(pi**2)))*a2**2 - &
&(a2/((3.D0)*a4))*(((1.D0/(16.D0*pi**4))*a2**2+(1.D0/(pi**2))*a4*(rho0-B))**0.5D0)


!! initial value for metric function phi 
phi0 = 0.1D0     ! arbitrary (needed to be adjusted later)
r0 = 1.D-30      ! integration starting point

!! Set initial conditions
!!!!!!!!!!!!!!!!!
!!Start integration loop
!!!!!!!!!!!!!!!!!
r = r0
y(1) = 0.D0
y(2) = P0
y(3) = phi0
y(4) = 1/(3*E0)
y(5) = 1
omega = 2*pi*1000/(2.997D5) !omega of 1kHz in code unit 
DO l = 1, 1000
    omega = omega + omegastep !shooting method part
DO i = 1, 1000000000

   rend = r0 + REAL(i)*step
   call oderk(r,rend,y,N_ode,fcn)

   r = rend
   mass = y(1)
   P = y(2)
   phi = y(3)
   xi = y(4)
   eta = y(5)

   IF (P < tiny*P0) THEN
      WRITE(*,*) "Central density (10^14 cgs) = ",  rho0_cgs/1.D14
      WRITE(*,*) " Mass (solar mass) = ", mass/1.477D0
      WRITE(*,*) " Radius (km) = ", r
      WRITE(*,*) " Compactness M/R ", mass/r
      WRITE(15,*) (omega*2.997D5/(2*pi)), y(5)

      GOTO 21
   ENDIF

ENDDO
ENDDO
21 CONTINUE

END PROGRAM roscillation2

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE fcn(r,y,yprime)
USE eos_parameters
IMPLICIT NONE 
REAL(8), DIMENSION(5) :: y, yprime
REAL(8) :: r, m, P, phi, rho, pi, B, a2, xi, eta, W, Q, E, s, lamda, omega
INTEGER :: j 

pi = ACOS(-1.D0)
a2 =((((1.6022D-13)**4)*(6.674D-11)*((2.997D8)**-7)*((1.0546D-34)**-3)*(1.D6))**(0.5D0))*a2_MeV !convert to code unit (km^-1)
B = ((1.6022D-13)**4)*(6.674D-11)*((2.997D8)**-7)*((1.0546D-34)**-3)*(1.D6)*B_MeV !convert to code unit (km^-2)
m = y(1)
P = y(2)
phi = y(3)
xi = y(4)
eta = y(5)

rho = 3.D0*P + 4.D0*B +((3.D0)/(4.D0*a4*(pi**2)))*a2**2+(a2/a4)*&
&(((9.D0/((16.D0)*(pi**4)))*a2**2+((3.D0/(pi**2))*a4*(P+B)))**0.5D0)

s = (1.D0/3.D0) - (1/(6*pi**2))*a2*((1/(16*pi**2)*a2**2 + (pi**-2)*a4*(rho - B))**-0.5) !square of speed of sound

W = (r**-2)*(rho + P)*exp(3*lamda + phi)

E = (r**-2)*s*exp(lamda + 3*phi)

Q = (r**-2)*exp(lamda + 3*phi)*(rho + P)*((yprime(3)**2) + 4*(r**-1)*yprime(3)- 8*pi*P*exp(2*lamda))

yprime(1) = 4.D0*pi*rho*r**2

yprime(2) = - (rho + P)*(m + 4.D0*pi*P*r**3)/(r*(r-2.D0*m))

yprime(3) = (m + 4.D0*pi*P*r**3)/(r*(r-2.D0*m)) 

yprime(4) = y(5)/(3*E)

yprime(5) = -(W*omega**2 + Q)*y(4)

END SUBROUTINE fcn

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!
!! Runge-Kutta method (from Numerical Recipes)
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine oderk(ri,re,y,n,derivs) 
INTEGER, PARAMETER :: NMAX=16
REAL(8) :: ri, re, step
REAL(8), DIMENSION(NMAX) :: y, dydx, yout
EXTERNAL :: derivs,rk4 

call derivs(ri,y,dydx) 
step=re-ri 
CALL rk4(y,dydx,n,ri,step,yout,derivs) 
do i=1,n 
   y(i)=yout(i) 
enddo
return 
end subroutine oderk

SUBROUTINE RK4(Y,DYDX,N,X,H,YOUT,DERIVS) 
INTEGER, PARAMETER :: NMAX=16 
REAL(8) :: H,HH,XH,X,H6 
REAL(8), DIMENSION(N) :: Y, DYDX, YOUT 
REAL(8), DIMENSION(NMAX) :: YT, DYT, DYM 
EXTERNAL :: derivs


HH=H*0.5D0 
H6=H/6D0 
XH=X+HH

DO I=1,N 
   YT(I)=Y(I)+HH*DYDX(I) 
ENDDO

CALL DERIVS(XH,YT,DYT) 

DO I=1,N 
   YT(I)=Y(I)+HH*DYT(I) 
ENDDO

CALL DERIVS(XH,YT,DYM) 

DO I=1,N 
   YT(I)=Y(I)+H*DYM(I) 
   DYM(I)=DYT(I)+DYM(I) 
ENDDO

CALL DERIVS(X+H,YT,DYT) 

DO I=1,N 
   YOUT(I)=Y(I)+H6*(DYDX(I)+DYT(I)+2*DYM(I)) 
ENDDO


END SUBROUTINE RK4