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
|
/*@@
@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
DECLARE_CCTK_FUNCTIONS
logical firstcall
integer qtype
CCTK_REAL brillq,rho,z,phi
data firstcall /.true./
save firstcall,qtype
c Get parameters at first call.
if (firstcall) then
if (CCTK_EQUALS(q_function,"exp")) then
qtype = 0
else if (CCTK_EQUALS(q_function,"eppley")) then
qtype = 1
else if (CCTK_EQUALS(q_function,"gundlach")) then
qtype = 2
else
call CCTK_WARN(0,"Brill wave data type not recognised")
end if
firstcall = .false.
end if
if (rho < 0) then
rho = -rho
end if
c Calculate q(rho,z) from a choice of forms.
c "exp" brill data
if (qtype.eq.0) then
brillq = exp_a*dabs(rho)**(2 + exp_b)
. / dexp((rho - exp_rho0)**2)/(rho**2 + z**2)
if (exp_sigmaz.ne.0.0D0) then
brillq = brillq/dexp((z/exp_sigmaz)**2)
end if
else if (qtype.eq.1) then
c "Eppley" brill data. This includes Eppleys choice of q.
c But note that q(Eppley) = 2q(Cactus).
brillq = eppley_a*(rho/eppley_sigmarho)**eppley_b
& / ( 1.0D0 + ((rho**2 + z**2 - eppley_r0**2)/
& eppley_sigmar**2)**(eppley_c/2) )
else if (qtype.eq.2) then
c "Gundlach" brill data This includes my (Carstens) notion of what a
c smooth "pure quadrupole" q should be, plus the choice of
c Holz et al.
brillq = gundlach_a*(rho/gundlach_sigmarho)**gundlach_b
& / dexp(((rho**2 + z**2 - gundlach_r0**2)/
& gundlach_sigmar**2)**(gundlach_c/2))
else
c Unknown type for q function.
call CCTK_WARN(0,"brillq: Unknown type of Brill function q")
end if
c If desired, multiply with a phi dependent factor.
if (CCTK_EQUALS(initial_data,"brilldata")) then
brillq = brillq*(1.0D0 + brill3d_d*rho**brill3d_m*
& cos(brill3D_n*(phi-brill3d_phi0))**2
& / (1.0D0 + brill3d_e*rho**brill3d_m))
end if
return
end
|
