1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
|
#ifdef HAVE_ODEPACK
c------------------------------------------------------------------------
c equation of state, converts p to rho
c currently using a polytrope
c------------------------------------------------------------------------
real*8 function p2rho(p)
implicit none
real*8 p
p2rho = p ** (3/4.0d0)
return
end
c------------------------------------------------------------------------
c equation of state, converts rho to p
c currently using a polytrope
c------------------------------------------------------------------------
real*8 function rho2p(rho)
implicit none
real*8 rho
rho2p = rho ** (4/3.0d0)
return
end
c-----------------------------------------------------------
c Integrates static equations of motion for spherically
c symmetric, general-relativistic TOV star:
c
c dp/dr = - (p + rho)*(m+4*pi*r**3*p)/ r / (r-2*m)
c dm/dr = 4*pi*r**2 * rho
c d alpha/dr = alpha*(m+4*pi*r**3*rho )/ r / (r-2*m)
c
c (see, e.g. Walds textbook pp126-127.)
c where g_tt = - alpha**2, and g_rr = 1/(1-2m/r)
c
c Note that (my rho) = (the-rest-of-GRHydros rho)*(1+epsilon)
c
c-----------------------------------------------------------
c neq = 3
c-----------------------------------------------------------
subroutine tov_fcn(neq,r,y,yprime)
implicit none
real*8 rho0
common / com_tov / rho0
real*8 p2rho
integer neq
real*8 r, y(neq), yprime(neq)
real*8 pi
parameter ( pi = 3.1415926 5358 9793 d0 )
real*8 p, m, alf, rho, fourpr2
p = y(1)
m = y(2)
alf = y(3)
if( r .eq. 0.0d0 ) then
yprime(1) = 0.0d0
yprime(2) = 0.0d0
yprime(3) = 0.0d0
else
fourpr2 = 4*pi*r*r
c BUG: should be using call to CactusEOS
rho = p2rho(p)
yprime(1) = - (p + rho)*(m+fourpr2*r*p)/ r / (r-2*m)
yprime(2) = fourpr2 * rho
yprime(3) = alf*(m+fourpr2*r*rho)/ r / (r-2*m)
endif
return
end
c-----------------------------------------------------------
c Implements Jacobian (optional) ...
c Does absolutely nothing
c-----------------------------------------------------------
subroutine tov_jac
implicit none
real*8 rho0
common / com_tov / rho0
return
end
#else
subroutine tov_fcn()
return
end
#endif
|
