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
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
|
c The metric given here corresponds to a constant
c density star, also known as a "Schwarzschild" star.
c There is corresponding code in
c include/Scalar_CalcTmunu.inc
c to set up the matter variables.
C
c Author: Mitica Vulcanov
C see ../../README for copyright & licensing info
C
c $Header$
c
c The metric is given as a conformally flat metric.
c Turns out that in the original areal radius, the
c metric variables have a kink at the surface of the
c star, but they are smooth in the conformal form.
c
c Thanks to Philippos Papadopoulos for suggesting
c the use of this metric.
#include "cctk.h"
#include "cctk_Parameters.h"
subroutine Exact__constant_density_star(
$ x, y, z, t,
$ gdtt, gdtx, gdty, gdtz,
$ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
$ gutt, gutx, guty, gutz,
$ guxx, guyy, guzz, guxy, guyz, guzx,
$ psi, Tmunu_flag)
implicit none
DECLARE_CCTK_PARAMETERS
c input arguments
CCTK_REAL x, y, z, t
c output arguments
CCTK_REAL gdtt, gdtx, gdty, gdtz,
$ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
$ gutt, gutx, guty, gutz,
$ guxx, guyy, guzz, guxy, guyz, guzx
CCTK_REAL psi
LOGICAL Tmunu_flag
c local static variables
logical firstcall
CCTK_REAL mass,radius
data firstcall /.true./
save firstcall, mass, radius
c local variables
CCTK_REAL r,c,psi4
c constants
CCTK_REAL zero,one,two
parameter (zero=0.0d0, one=1.0d0, two=2.0d0)
c This model has matter
c ==> it sets the stress-energy tensor in the "CalcTmunu" code
Tmunu_flag = .true.
c Get parameters of the metric.
if (firstcall) then
mass = constant_density_star__mass
radius = constant_density_star__radius
firstcall = .false.
end if
c Find r.
r = sqrt(x**2 + y**2 + z**2)
c Find conformal factor.
if (r.le.radius) then
c = mass/(two*radius)
psi4 = (one + c)**6/(one + c*(r/radius)**2)**2
else
c = mass/(two*r)
psi4 = (one + c)**4
end if
c Find metric components.
gdtt = - psi4
gdtx = zero
gdty = zero
gdtz = zero
gdxx = psi4
gdyy = psi4
gdzz = psi4
gdxy = zero
gdyz = zero
gdzx = zero
c Find inverse metric.
gutt = -one/psi4
gutx = zero
guty = zero
gutz = zero
guxx = one/psi4
guyy = one/psi4
guzz = one/psi4
guxy = zero
guyz = zero
guzx = zero
return
end
|
