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
|
/* $Header$ */
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
static
CCTK_REAL min (CCTK_REAL const x, CCTK_REAL const y)
{
return x < y ? x : y;
}
static
CCTK_REAL max (CCTK_REAL const x, CCTK_REAL const y)
{
return x > y ? x : y;
}
void SphericalSurface_Set (CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
CCTK_REAL const pi = 3.1415926535897932384626433832795028841971693993751;
int n;
int i, j;
for (n=0; n<nsurfaces; ++n) {
if (set_spherical[n]) {
sf_valid[n] = +1;
sf_area[n] = 4 * pi * pow (radius[n], 2);
sf_mean_radius[n] = radius[n];
sf_centroid_x[n] = origin_x[n];
sf_centroid_y[n] = origin_y[n];
sf_centroid_z[n] = origin_z[n];
sf_quadrupole_xx[n] = 0.0;
sf_quadrupole_xy[n] = 0.0;
sf_quadrupole_xz[n] = 0.0;
sf_quadrupole_yy[n] = 0.0;
sf_quadrupole_yz[n] = 0.0;
sf_quadrupole_zz[n] = 0.0;
sf_min_radius[n] = radius[n];
sf_max_radius[n] = radius[n];
sf_min_x[n] = origin_x[n] - radius[n];
sf_min_y[n] = origin_y[n] - radius[n];
sf_min_z[n] = origin_z[n] - radius[n];
sf_max_x[n] = origin_x[n] + radius[n];
sf_max_y[n] = origin_y[n] + radius[n];
sf_max_z[n] = origin_z[n] + radius[n];
for (j=0; j<nphi[n]; ++j) {
for (i=0; i<ntheta[n]; ++i) {
sf_radius[n] = radius[n];
}
}
sf_origin_x[n] = origin_x[n];
sf_origin_y[n] = origin_y[n];
sf_origin_z[n] = origin_z[n];
} else if (set_elliptic[n]) {
/*
* Equation for an ellipsoid:
* x^2 / rx^2 + y^2 / ry^2 + z^2/ rz^2 = 1
* where rx, ry, rz are the radii in the x, y, and z directions.
*
* With
* x = r (sin theta) (cos phi)
* y = r (sin theta) (sin phi)
* z = r (cos theta)
*
* we get the solution
* 1/r^2 = x^2/rx^2 + y^2/ry^2 + z^2/rz^2
*
* With the form
* x A x = 1
* we obtain
* A = diag [1/rx^2, 1/ry^2, 1/rz^2]
*/
CCTK_REAL const rx2 = pow (radius_x[n], 2);
CCTK_REAL const ry2 = pow (radius_y[n], 2);
CCTK_REAL const rz2 = pow (radius_z[n], 2);
sf_valid[n] = +1;
sf_area[n] = 0 * 4 * pi * pow (radius[n], 2);
sf_mean_radius[n] = 0 * radius[n];
sf_centroid_x[n] = origin_x[n];
sf_centroid_y[n] = origin_y[n];
sf_centroid_z[n] = origin_z[n];
sf_quadrupole_xx[n] = 1.0 / rz2;
sf_quadrupole_xy[n] = 0.0;
sf_quadrupole_xz[n] = 0.0;
sf_quadrupole_yy[n] = 1.0 / ry2;
sf_quadrupole_yz[n] = 0.0;
sf_quadrupole_zz[n] = 1.0 / rx2;
sf_min_radius[n] = min (radius_x[n], min (radius_y[n], radius_z[n]));
sf_max_radius[n] = max (radius_x[n], max (radius_y[n], radius_z[n]));
sf_min_x[n] = origin_x[n] - radius_x[n];
sf_min_y[n] = origin_y[n] - radius_y[n];
sf_min_z[n] = origin_z[n] - radius_z[n];
sf_max_x[n] = origin_x[n] + radius_x[n];
sf_max_y[n] = origin_y[n] + radius_y[n];
sf_max_z[n] = origin_z[n] + radius_z[n];
for (j=0; j<nphi[n]; ++j) {
for (i=0; i<ntheta[n]; ++i) {
CCTK_REAL const theta = sf_origin_theta[n] + i * sf_delta_theta[n];
CCTK_REAL const phi = sf_origin_phi[n] + j * sf_delta_phi[n];
CCTK_REAL const x2 = pow (sin(theta) * cos(phi), 2);
CCTK_REAL const y2 = pow (sin(theta) * sin(phi), 2);
CCTK_REAL const z2 = pow (cos(theta) , 2);
sf_radius[n] = 1.0 / sqrt (x2 / rx2 + y2 / ry2 + z2 / rz2);
}
}
}
sf_origin_x[n] = origin_x[n];
sf_origin_y[n] = origin_y[n];
sf_origin_z[n] = origin_z[n];
} /* for n */
}
|
