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
|
c/*@@
c @file AHFinder_fun.F
c @date April 1998
c @author Miguel Alcubierre
c @desc
c Find horizon function.
c @enddesc
c@@*/
c Note that including cactus.h will also include AHFinder.h
!#include "cactus.h"
#include "cctk.h"
#include "cctk_parameters.h"
#include "cctk_arguments.h"
subroutine AHFinder_fun(CCTK_FARGUMENTS)
use AHFinder_dat
implicit none
DECLARE_CCTK_FARGUMENTS
DECLARE_CCTK_PARAMETERS
integer i,j,k
integer l,m
CCTK_REAL LEGEN
CCTK_REAL xp,yp,zp,rp
CCTK_REAL phi,cost,cosa,sina
CCTK_REAL zero,half,one,two
CCTK_REAL pi,halfpi,twopi
CCTK_REAL aux1,aux2
! **************************
! *** DEFINE NUMBERS ***
! **************************
zero = 0.0D0
half = 0.5D0
one = 1.0D0
two = 2.0D0
nx = cctk_lsh(1)
ny = cctk_lsh(2)
nz = cctk_lsh(3)
! *************************************
! *** DEFINE {pi,halfpi,twopi} ***
! *************************************
pi = 3.141592654D0
halfpi = half*pi
twopi = two*pi
! *********************************
! *** FIND HORIZON FUNCTION ***
! *********************************
do k=1,nz
do j=1,ny
do i=1,nx
xp = x(i,j,k) - xc
yp = y(i,j,k) - yc
zp = z(i,j,k) - zc
rp = dsqrt(xp**2 + yp**2 + zp**2)
! Find phi.
if (dabs(xp).gt.dabs(yp)) then
phi = atan(dabs(yp/xp))
else if (dabs(xp).lt.dabs(yp)) then
phi = halfpi - atan(dabs(xp/yp))
else
phi = 0.25D0*pi
end if
if ((xp.eq.zero).and.(yp.eq.zero)) then
phi = zero
else if ((xp.le.zero).and.(yp.ge.zero)) then
phi = pi - phi
else if ((xp.le.zero).and.(yp.le.zero)) then
phi = pi + phi
else if ((xp.ge.zero).and.(yp.le.zero)) then
phi = twopi - phi
end if
! Monopole term.
aux1 = c0(0)
! Axisymmetric terms.
if (rp.ne.zero) then
cost = zp/rp
else
cost = one
end if
do l=1+stepz,lmax,1+stepz
aux1 = aux1 + c0(l)*LEGEN(l,0,cost)
end do
! Non-axisymmetric terms. Notice how the sum over M is first.
! This will allow me to use the recursion relations to avoid
! having to start from scratch every time. Also, I sum over
! all l's even if I don't want some terms. This is because
! in order to use the recursion relations I need all polynomials.
if (nonaxi) then
do m=1,lmax
aux2 = dble(m)*phi
sina = sin(aux2)
cosa = cos(aux2)
do l=m,lmax
aux2 = LEGEN(l,m,cost)
aux1 = aux1 + aux2*cc(l,m)*cosa
if (.not.refy) then
aux1 = aux1 + aux2*cs(l,m)*sina
end if
end do
end do
end if
ahfgrid(i,j,k) = rp - aux1
! Check for NaN's.
if (ahfgrid(i,j,k).ge.zero) then
else if (ahfgrid(i,j,k).lt.zero) then
else
#ifdef MPI
if (myproc.eq.0) then
#endif
write(*,*) 'NaN in ahfgrid at point',i,j,k
#ifdef MPI
end if
#endif
STOP
end if
end do
end do
end do
! Synchronize.
#ifdef MPI
call synconefunc(ahfgrid)
#endif
! ***************
! *** END ***
! ***************
end subroutine AHFinder_fun
|
