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
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
|
/*@@
@file brilldata.F
@date
@author Carsten Gundlach, Miguel Alcubierre.
@desc
Construct Brill wave q function.
@enddesc
@version $Header$
@@*/
#include "cctk.h"
#include "cctk_Parameters.h"
function brillq(rho,z,phi)
c Calculates the function q that appear in the conformal
c metric for Brill waves:
c
c ds^2 = psi^4 ( e^(2q) (drho^2 + dz^2) + rho^2 dphi^2 )
c
c There are three different choices for the form of q depending
c on the value of the parameter "brill_q":
c
c brill_q = 0:
c 2
c 2+b - (rho-rho0) 2 2
c q = a rho e (z/sz) / r
c
c
c brill_q = 1: (includes Eppleys form as special case)
c
c
c b 2 2 2 c/2
c q = a (rho/srho) / { 1 + [ ( r - r0 ) / sr ] }
c
c
c brill_q = 2: (includes Holz et al form as special case)
c
c 2 2 2 c/2
c b - [ ( r - r0 ) / sr ]
c q = a (rho/srho) e
c
c
c If we want 3D initial data, the function q is multiplied by an
c additional factor:
c
c m 2 m
c q -> q [ 1 + d rho cos (n (phi)+phi0 ) / ( 1 + e rho ) ]
implicit none
DECLARE_CCTK_PARAMETERS
logical firstcall
integer qtype
integer b,c
CCTK_REAL brillq,rho,z,phi
CCTK_REAL a,r0,sr,rho0,srho,sz
CCTK_REAL d,e,m,n,phi0
data firstcall /.true./
save firstcall,qtype
save a,b,c,r0,sr,rho0,srho,sz
save d,e,m,n,phi0
c Get parameters at first call.
if (firstcall) then
qtype = brill_q
a = brill_a
b = brill_b
c = brill_c
r0 = brill_r0
sr = brill_sr
rho0 = brill_rho0
srho = brill_srho
sz = brill_sz
if (axisym.eq.0) then
d = brill_d
e = brill_e
m = brill_m
n = brill_n
phi0 = brill_phi0
end if
firstcall = .false.
end if
c Calculate q(rho,z) from a choice of forms.
c brill_q = 0.
if (qtype.eq.0) then
brillq = a*dabs(rho)**(2 + b)
. / dexp((rho - rho0)**2)/(rho**2 + z**2)
if (sz.ne.0.0D0) then
brillq = brillq/dexp((z/sz)**2)
end if
else if (qtype.eq.1) then
c brill_q = 1. This includes Eppleys choice of q.
c But note that q(Eppley) = 2q(Cactus).
brillq = a*(rho/srho)**b
. / ( 1.0D0 + ((rho**2 + z**2 - r0**2)/sr**2)**(c/2) )
else if (qtype.eq.2) then
c brill_q = 2. This includes my (Carstens) notion of what a
c smooth "pure quadrupole" q should be, plus the choice of
c Holz et al.
brillq = a*(rho/srho)**b
. / dexp(((rho**2 + z**2 - r0**2)/sr**2)**(c/2))
else
c Unknown type for q function.
call CCTK_WARN(0,"Unknown type of Brill function q")
end if
c If desired, multiply with a phi dependent factor.
if (axisym.eq.0) then
brillq = brillq*(1.0D0 + d*rho**m*cos(n*(phi-phi0))**2
. / (1.0D0 + e*rho**m))
end if
return
end
|
