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
|
// coords.hh -- classes for mpe coordinates
// $Id$
//
// coords - misc coordinates-stuff namespace
//
//
// prerequisites:
// <jt/stdc.h> -- for DEGREES_TO_RADIANS and RADIANS_TO_DEGREES
// fp.hh -- basic floating point stuff
//
//*****************************************************************************
//
// This namespace contains coordinate conversion functions, focusing
// on the angular coordinates. It also contains coordinate bit masks
// and generic coordinate mnemonics.
//
// The angular coordinates are
// mu = rotation about x axis = arctan(y/z)
// nu = rotation about y axis = arctan(x/z)
// phi = rotation about z axis = arctan(y/x)
//
// Terminology:
// (a,b,c) == any one of (a,b,c) == a | b | c
// ((a,b,c)) == any two of (a,b,c) == (a,b) | (a,c) | (b,c)
// ((r,a,b,c)) == r and any two of (a,b,c) == (r,a,b) | (r,a,c) | (r,b,c)
// {{a,b,c}} == any two of (a,b,c) separated by an underscore
// == a_b | a_c | b_c
// {{r,a,b,c}} == r, an underscore, then any two of (a,b,c) separated
// by an underscore
// == r_a_b | r_a_c | r_b_c
//
namespace coords
{
//
// coords:: is a *namespace*, not a class, so we need
// explicit inline declarations on inline functions to
// avoid (conflicting) duplicate defs in separate compilation
// units.
//
//
// FIXME:
// should the reference arguments in the following functions
// have restrict qualifiers?
//
// Bugs:
// The functions aren't optimally efficient -- they have
// lots of could-be-optimized-away divisions and trig calls.
// But in practice this doesn't matter, as we don't call
// these functions inside inner loops. (For example, we
// cache all the partial derivatives needed for interpatch
// tensor transformations, so these functions are only
// called to initialize the cache when setting up the
// patch system.)
//
// ((mu,nu,phi)) --> the 3rd
fp phi_of_mu_nu(fp mu, fp nu);
fp nu_of_mu_phi(fp mu, fp phi);
fp mu_of_nu_phi(fp nu, fp phi);
// partial derivatives of any of (mu,nu,phi)
// wrt any other with the 3rd held constant
// ... partial phi / partial (mu,nu)
fp partial_phi_wrt_mu_at_const_nu(fp mu, fp nu );
fp partial_phi_wrt_nu_at_const_mu(fp mu, fp nu );
// ... partial mu / partial (nu,phi)
fp partial_mu_wrt_nu_at_const_phi(fp nu, fp phi);
fp partial_mu_wrt_phi_at_const_nu(fp nu, fp phi);
// ... partial nu / partial (mu,phi)
fp partial_nu_wrt_mu_at_const_phi(fp mu, fp phi);
fp partial_nu_wrt_phi_at_const_mu(fp mu, fp phi);
// 2nd partial derivatives of (mu,nu,phi) among themselves
// ... partial^2 phi / partial (mu,nu)^2
fp partial2_phi_wrt_mu2_at_const_nu (fp mu, fp nu);
fp partial2_phi_wrt_mu_nu_at_const_nu_mu (fp mu, fp nu);
fp partial2_phi_wrt_nu2_at_const_mu (fp mu, fp nu);
// ... partial^2 mu / partial (nu,phi)^2
fp partial2_mu_wrt_nu2_at_const_phi (fp nu, fp phi);
fp partial2_mu_wrt_nu_phi_at_const_phi_nu(fp nu, fp phi);
fp partial2_mu_wrt_phi2_at_const_nu (fp nu, fp phi);
// ... partial^2 nu / partial (mu,phi)^2
fp partial2_nu_wrt_mu2_at_const_phi (fp mu, fp phi);
fp partial2_nu_wrt_mu_phi_at_const_phi_mu(fp mu, fp phi);
fp partial2_nu_wrt_phi2_at_const_mu (fp mu, fp phi);
// (r,(mu,nu,phi)) <--> (x,y,z)
// ... on unit sphere
void xyz_of_mu_nu(fp mu, fp nu, fp &x, fp &y, fp &z);
void xyz_of_mu_phi(fp mu, fp phi, fp &x, fp &y, fp &z);
void xyz_of_nu_phi(fp nu, fp phi, fp &x, fp &y, fp &z);
void mu_nu_of_xyz(fp x, fp y, fp z, fp &mu, fp &nu);
void mu_phi_of_xyz(fp x, fp y, fp z, fp &mu, fp &phi);
void nu_phi_of_xyz(fp x, fp y, fp z, fp &nu, fp &phi);
// ... generic r
void xyz_of_r_mu_nu(fp r, fp mu, fp nu, fp &x, fp &y, fp &z);
void xyz_of_r_mu_phi(fp r, fp mu, fp phi, fp &x, fp &y, fp &z);
void xyz_of_r_nu_phi(fp r, fp nu, fp phi, fp &x, fp &y, fp &z);
void r_mu_nu_of_xyz(fp x, fp y, fp z, fp &r, fp &mu, fp &nu);
void r_mu_phi_of_xyz(fp x, fp y, fp z, fp &r, fp &mu, fp &phi);
void r_nu_phi_of_xyz(fp x, fp y, fp z, fp &r, fp &nu, fp &phi);
// usual polar spherical (r,theta,phi) <--> (x,y,z)
// ... on unit sphere
void xyz_of_theta_phi(fp theta, fp phi, fp &x, fp &y, fp &z);
void theta_phi_of_xyz(fp x, fp y, fp z, fp &theta, fp &phi);
// ... generic r
void xyz_of_r_theta_phi(fp r, fp theta, fp phi, fp &x, fp &y, fp &z);
void r_theta_phi_of_xyz(fp x, fp y, fp z, fp &r, fp &theta, fp &phi);
// ((mu,nu,phi)) <--> usual polar spherical (theta,phi)
// ... note phi is the same coordinate in both systems
void theta_phi_of_mu_nu(fp mu, fp nu, fp &ps_theta, fp &ps_phi);
void theta_phi_of_mu_phi(fp mu, fp phi, fp &ps_theta, fp &ps_phi);
void theta_phi_of_nu_phi(fp nu, fp phi, fp &ps_theta, fp &ps_phi);
void mu_nu_of_theta_phi(fp ps_theta, fp ps_phi, fp &mu, fp &nu);
void mu_phi_of_theta_phi(fp ps_theta, fp ps_phi, fp &mu, fp &phi);
void nu_phi_of_theta_phi(fp ps_theta, fp ps_phi, fp &nu, fp &phi);
//*************************************
//
// We want to manipulate tensor indices, and also do calculations
// like "which coordinate do these two patches have in common".
// We implement the former in the obvious way (a tensor index is
// just an integer), and the latter by Boolean operations on
// integers using the coords_set bit vectors below:
//
//
// (mu,nu,phi) coordinates
// ... in a non-namespace context we use the prefix mnp_
// to identify things relating to these coordinates
//
namespace mu_nu_phi
{
static const int N_coords = 3;
// tensor indices
static const int index_mu = 0;
static const int index_nu = 1;
static const int index_phi = 2;
// sets of coordinates
typedef int coords_set;
// singleton coordinate set for a given index
inline coords_set coord_set_of_index(int i) { return 0x1 << i; }
static const coords_set set_mu = coord_set_of_index(index_mu );
static const coords_set set_nu = coord_set_of_index(index_nu );
static const coords_set set_phi = coord_set_of_index(index_phi);
static const coords_set set_empty = 0;
static const coords_set set_all = set_mu | set_nu | set_phi;
// human-readable coordinate names for debugging etc
const char *name_of_index(int c);
const char *name_of_coords_set(coords_set S);
// set complement of coordinates
inline coords_set coords_set_not(coords_set S)
{ return set_all - S; }
};
//
// (r,rho,sigma) coordinates
// ... in a non-namespace context we use the prefix rrs_
// to identify things relating to these coordinates
//
namespace r_rho_sigma
{
static const int N_coords = 3;
static const int N_angular_coords = 2;
// tensor indices
static const int index_r = 0;
static const int index_rho = 1;
static const int index_sigma = 2;
// human-readable coordinate names for debugging etc
const char *name_of_index(int c);
};
//*************************************
//
// degrees <--> radians conversions
//
// ... note coords:: is a *namespace*, not a class,
// so we need an explicit inline declaration to avoid
// conflicting duplicate defs in separate compilation units
inline fp dang_of_ang(fp ang) { return RADIANS_TO_DEGREES * ang; }
inline fp ang_of_dang(fp dang) { return DEGREES_TO_RADIANS * dang; }
};
|
