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
|
/*@@
@file AHFinder_fun.F
@date April 1998
@author Miguel Alcubierre
@desc
Find horizon function.
@enddesc
@@*/
#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
DECLARE_CCTK_FUNCTIONS
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
CCTK_REAL red_tmp
! **************************
! *** DEFINE NUMBERS ***
! **************************
zero = 0.0D0
half = 0.5D0
one = 1.0D0
two = 2.0D0
pi = acos(-one)
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 even if I do not 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
end do
end do
end do
! ***************
! *** END ***
! ***************
end subroutine AHFinder_fun
|