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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
|
/*@@
@file planewaves.f
@date Wed Jan 24 16:51:03 1996
@author Malcolm Tobias + Joan Masso
@desc
Routine to initialize linear plane wave spacetime
@enddesc
@hdate Tue Mar 16 12:38:04 1999 @hauthor Gerd Lanfermann
@hdesc Converted to Cactus4.0. the code and its comments of the
original authors are kept. There is NO BM suport at this point.
Id like to know how Cactus handles this w/o #ifdef !
@@*/
/*@@
@routine planewaves
@date Mon Jan 29 11:57:00 1996
@author Malcolm Tobias + Joan Masso
@desc
Routine to initialize a linear plane wave travelling in an
arbitrary direction. The form of the packet is: <p>
$ A*exp\left[-(kp_ix^i-\omega_p (time-ra))^2\right]
cos(k_ix^i-\omega \ time) $
<p>
where: <p>
A = amplitude of the wave <br>
k = the wave # of the sine wave <br>
$\omega$ = the frequency of the sine wave <br>
kp = the wave # of the gaussian modulating the sine wave<br>
$\omega_p$ = the frequency of the gaussian <br>
ra = the initial position of the packet(s)
@enddesc
@calledby linearWaves
@calls
@@*/
#include "cctk.h"
#include "declare_arguments.h"
#include "declare_parameters.h"
subroutine planewaves(CCTK_FARGUMENTS)
implicit none
DECLARE_CCTK_FARGUMENTS
DECLARE_PARAMETERS
INTEGER iin,iout
REAL pi
REAL amplitude,ra,the,phi
REAL wave,wavep
REAL kx,ky,kz,w
REAL kxp,kyp,kzp,wp
c local arrays (scalars! Man, what a sweat!)
REAL plus,minus,plusp,minusp,ain,aout
INTEGERE i,j,k
pi = 3.14159265358979d0
c Wave characteristics
c amplitude is declared by DECLARE_PARAMETER, for all other parameters
c we stick to the naming in the cactus3.2.0 version of cactus
c wavecentering
ra = wavecenter
c wavelength
wave = wavelength
c wavepulse
wavep= wavepulse
c determine whether the wave is ingoing, outgoing, or both
iin = 0
iout = 0
if(CCTK_Equals(wavesgoing,'in').ne.0) then
iin = 1
elseif(CCTK_Equals(wavesgoing,'out').ne.0) then
iout = 1
elseif(CCTK_Equals(wavesgoing,'both').ne.0) then
iin = 1
iout = 1
endif
c determine the direction for the plane wave to travel
c and convert it from degrees to radians
the = pi*wavetheta/180.d0
phi = pi*wavephi/180.d0
write(*,*)'Plane waves'
write(*,*)'amplitude = ',amplitude
write(*,*)'wavecenter = ',ra
write(*,*)'wavelength = ',wave
write(*,*)'wavepulse = ',wavep
c precalc
kx = 2*pi*sin(the)*cos(phi)/wave
ky = 2*pi*sin(the)*sin(phi)/wave
kz = 2*pi*cos(the)/wave
w = (kx*kx+ky*ky+kz*kz)**0.5
kxp = 2*pi*sin(the)*cos(phi)/wavep
kyp = 2*pi*sin(the)*sin(phi)/wavep
kzp = 2*pi*cos(the)/wavep
wp = (kxp*kxp+kyp*kyp+kzp*kzp)**0.5
c *************** plane waves ********************
do k=1,sh(3)
do j=1,sh(2)
do i=1,sh(1)
plus = (kx*x(i,j,k)+ky*y(i,j,k)+kz*z(i,j,k)+w*time)
minus = (kx*x(i,j,k)+ky*y(i,j,k)+kz*z(i,j,k)-w*time)
plusp =(kxp*x(i,j,k)+kyp*y(i,j,k)+kzp*z(i,j,k)+wp*(time-ra))
minusp =(kxp*x(i,j,k)+kyp*y(i,j,k)+kzp*z(i,j,k)-wp*(time-ra))
ain = iin*amplitude*cos(plus)*exp(-(plusp)**2)
aout = iout*amplitude*cos(minus)*exp(-(minusp)**2)
c the metric functions
gxx(i,j,k) = cos(the)**2*cos(phi)**2*(1+ain+aout) +
& sin(phi)**2*(1-ain-aout) + sin(the)**2*cos(phi)**2
gxy(i,j,k) = cos(the)**2*sin(phi)*cos(phi)*(1+ain+aout) +
& sin(phi)*cos(phi)*(1-ain-aout) - sin(the)**2*sin(phi)*cos(phi)
gxz(i,j,k) = sin(the)*cos(the)*cos(phi)*(1+ain+aout) -
& sin(the)*cos(the)*cos(phi)
gyy(i,j,k) = cos(the)**2*sin(phi)**2*(1+ain+aout) +
& cos(phi)**2*(1-ain-aout) + sin(the)**2*sin(phi)**2
gyz(i,j,k) = -sin(the)*cos(the)*sin(phi)*(1+ain+aout) +
& sin(the)*cos(the)*sin(phi)
gzz(i,j,k) = sin(the)**2*(1+ain+aout) + cos(the)**2
c At this pint, the CAC3.2.0 version has some code for the dxgxx,....
c defined, how does CAC4.0 handles this without #ifdef THORN_BONAMASSO ?
c and finally the extrinsic curvatures
c Joan: I removed the division by alp as the initialization
c of the lapse may take place later. This data is consistent
c for initial lapse equal one.
hxx(i,j,k) = amplitude*
& (cos(the)**2*cos(phi)**2*
& (iin*(-w*sin(plus)*exp(-(plusp)**2) -
& 2.d0*wp*plusp*cos(plus)*exp(-(plusp)**2)) +
& iout*(w*sin(minus)*exp(-(minusp)**2) +
& 2.d0*wp*minusp*cos(minus)*exp(-(minusp)**2))) -
& sin(phi)**2*
& (iin*(-w*sin(plus)*exp(-(plusp)**2) -
& 2.d0*wp*plusp*cos(plus)*exp(-(plusp)**2)) +
& iout*(w*sin(minus)*exp(-(minusp)**2) +
& 2.d0*wp*minusp*cos(minus)*exp(-(minusp)**2))))/(-2.d0)
hxy(i,j,k) = amplitude*
& (cos(the)**2*sin(phi)*cos(phi)*
& (iin*(-w*sin(plus)*exp(-(plusp)**2) -
& 2.d0*wp*plusp*cos(plus)*exp(-(plusp)**2)) +
& iout*(w*sin(minus)*exp(-(minusp)**2) +
& 2.d0*wp*minusp*cos(minus)*exp(-(minusp)**2))) -
& sin(phi)*cos(phi)*
& (iin*(-w*sin(plus)*exp(-(plusp)**2) -
& 2.d0*wp*plusp*cos(plus)*exp(-(plusp)**2)) +
& iout*(w*sin(minus)*exp(-(minusp)**2) +
& 2.d0*wp*minusp*cos(minus)*exp(-(minusp)**2))))/(-2.d0)
hxz(i,j,k) = amplitude*
& (sin(the)*cos(the)*cos(phi)*
& (iin*(-w*sin(plus)*exp(-(plusp)**2) -
& 2.d0*wp*plusp*cos(plus)*exp(-(plusp)**2)) +
& iout*(w*sin(minus)*exp(-(minusp)**2) +
& 2.d0*wp*minusp*cos(minus)*exp(-(minusp)**2))))/(-2.d0)
hyy(i,j,k) = amplitude*
& (cos(the)**2*sin(phi)**2*
& (iin*(-w*sin(plus)*exp(-(plusp)**2) -
& 2.d0*wp*plusp*cos(plus)*exp(-(plusp)**2)) +
& iout*(w*sin(minus)*exp(-(minusp)**2) +
& 2.d0*wp*minusp*cos(minus)*exp(-(minusp)**2))) -
& cos(phi)**2*
& (iin*(-w*sin(plus)*exp(-(plusp)**2) -
& 2.d0*wp*plusp*cos(plus)*exp(-(plusp)**2)) +
& iout*(w*sin(minus)*exp(-(minusp)**2) +
& 2.d0*wp*minusp*cos(minus)*exp(-(minusp)**2))))/(-2.d0)
hyz(i,j,k) = amplitude*
& (-sin(the)*cos(the)*sin(phi)*
& (iin*(-w*sin(plus)*exp(-(plusp)**2) -
& 2.d0*wp*plusp*cos(plus)*exp(-(plusp)**2)) +
& iout*(w*sin(minus)*exp(-(minusp)**2) +
& 2.d0*wp*minusp*cos(minus)*exp(-(minusp)**2))))/(-2.d0)
hzz(i,j,k) = amplitude*
& (sin(the)**2*
& (iin*(-w*sin(plus)*exp(-(plusp)**2) -
& 2.d0*wp*plusp*cos(plus)*exp(-(plusp)**2)) +
& iout*(w*sin(minus)*exp(-(minusp)**2) +
& 2.d0*wp*minusp*cos(minus)*exp(-(minusp)**2))))/(-2.d0)
c loop over sh ends here:
enddo
enddo
enddo
return
end
|
