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
|
// miscfp.cc -- misc floating-point functions
//
// jtutil::signum - sign of a floating point number
// jtutil::hypot3 - 3D Euclidean distance
// jtutil::arctan_xy - 4-quadrant arc tangent
// jtutil::modulo_reduce - reduce an angle modulo 2*pi radians (360 degrees)
// jtutil::zero_C_array - set a C-style array to all zeros
//
// ***** template instantiations *****
//
#include <math.h>
#include <stdlib.h>
// we want to instantiate templates with CCTK_* types
#ifdef STANDALONE_TEST
#include "fake_cctk.h"
#else
#include "cctk.h"
#endif
#include "stdc.h"
#include "util.hh"
//******************************************************************************
//******************************************************************************
//******************************************************************************
//
// This function computes the floating point "signum" function (as in APL),
// signum(x) = -1.0, 0.0, or +1.0, according to the sign of x.
//
namespace jtutil
{
double signum(double x)
{
if (x == 0.0)
then return 0.0;
else return (x > 0.0) ? 1.0 : -1.0;
}
} // namespace jtutil::
//******************************************************************************
//
// This function computes the 3-dimensional Euclidean distance,
// analogously to the standard math library function hypot(2) .
//
// Arguments:
// (x,y,z) = (in) The rectangular coordinates of a point in $\Re^3$.
//
// Result:
// The function returns the Euclidean distance of (x,y,z) from the origin.
//
// Bugs:
// Unlike hypot(2), this function takes no special care to avoid
// unwarranted IEEE exceptions if any of |x|, |y|, or |z| is close to
// the overflow and/or underflow threshold.
//
namespace jtutil
{
double hypot3(double x, double y, double z)
{
return sqrt(x*x + y*y + z*z);
}
} // namespace jtutil::
//******************************************************************************
//
// This function computes a four-quadrant arctangent, using what I think
// is the "right" ordering of the arguments, and returning 0.0 if both
// arguments are 0.0.
//
// Arguments:
// (x,y) = (in) The rectangular coordinates of a point in $\Re^2$.
//
// Result:
// The function returns a value $\theta$ such that $x + iy = R e^{i\theta}$ for
// some real $R$, i.e. it returns the angle between the positive $x$ axis and
// the line joining the origin and the point $(x,y)$.
//
namespace jtutil
{
double arctan_xy(double x, double y)
{
// note reversed argument order (y,x) in std::atan2() function
return ((x == 0.0) && (y == 0.0)) ? 0.0 : atan2(y,x);
}
} // namespace jtutil::
//******************************************************************************
//
// This function reduces x modulo xmod to be (fuzzily) in the range
// [xmin, xmax] , or does an error_exit() if no such value exists.
//
namespace jtutil
{
double modulo_reduce(double x, double xmod, double xmin, double xmax)
{
double xx = x;
while (fuzzy<double>::LT(xx, xmin))
{
xx += xmod;
}
while (fuzzy<double>::GT(xx, xmax))
{
xx -= xmod;
}
if (! (fuzzy<double>::GE(xx, xmin) && fuzzy<double>::LE(xx, xmax)) )
then jtutil::error_exit(ERROR_EXIT,
"***** modulo_reduce(): no modulo value is fuzzily within specified range!\n"
" x = %g xmod = %g\n"
" [xmin,xmax] = [%g,%g]\n"
" ==> xx = %g\n"
,
x, xmod,
xmin, xmax,
xx); /*NOTREACHED*/
return xx;
}
} // namespace jtutil::
//******************************************************************************
//
// This function sets a C-style array to all zeros.
//
namespace jtutil
{
template <typename fp_t>
void zero_C_array(int N, fp_t array[])
{
for (int i = 0 ; i < N ; ++i)
{
array[i] = 0;
}
}
} // namespace jtutil::
//******************************************************************************
//******************************************************************************
//******************************************************************************
//
// ***** template instantiations *****
//
namespace jtutil
{
template void zero_C_array<CCTK_REAL>(int, CCTK_REAL[]);
} // namespace jtutil::
|
