Line data Source code
1 : #include "crpropa/module/Boundary.h"
2 : #include "crpropa/Units.h"
3 :
4 : #include <sstream>
5 :
6 : namespace crpropa {
7 :
8 0 : PeriodicBox::PeriodicBox() :
9 0 : origin(Vector3d(0, 0, 0)), size(Vector3d(0, 0, 0)) {
10 0 : }
11 :
12 2 : PeriodicBox::PeriodicBox(Vector3d o, Vector3d s) :
13 2 : origin(o), size(s) {
14 2 : }
15 :
16 2 : void PeriodicBox::process(Candidate *c) const {
17 2 : Vector3d pos = c->current.getPosition();
18 2 : Vector3d n = ((pos - origin) / size).floor();
19 :
20 2 : if ((n.x == 0) and (n.y == 0) and (n.z == 0))
21 : return; // do nothing if candidate is inside the box
22 :
23 2 : c->current.setPosition(pos - n * size);
24 2 : c->previous.setPosition(c->previous.getPosition() - n * size);
25 2 : c->source.setPosition(c->source.getPosition() - n * size);
26 2 : c->created.setPosition(c->created.getPosition() - n * size);
27 : }
28 :
29 0 : void PeriodicBox::setOrigin(Vector3d o) {
30 : origin = o;
31 0 : }
32 0 : void PeriodicBox::setSize(Vector3d s) {
33 : size = s;
34 0 : }
35 :
36 0 : std::string PeriodicBox::getDescription() const {
37 0 : std::stringstream s;
38 0 : s << "Periodic box: origin " << origin / Mpc << " Mpc, ";
39 0 : s << "size " << size / Mpc << " Mpc";
40 0 : return s.str();
41 0 : }
42 :
43 :
44 1 : ReflectiveShell::ReflectiveShell(Vector3d center, double r) :
45 1 : center(center), radius(r) {
46 1 : }
47 :
48 2 : double ReflectiveShell::distance(const Vector3d &point) const {
49 : Vector3d dR = point - center;
50 2 : return dR.getR() - radius;
51 : }
52 :
53 1 : Vector3d ReflectiveShell::normal(const Vector3d& point) const {
54 : Vector3d d = point - center;
55 1 : return d.getUnitVector();
56 : }
57 :
58 1 : void ReflectiveShell::process(Candidate *c) const {
59 1 : double currentDistance = distance(c->current.getPosition());
60 1 : double previousDistance = distance(c->previous.getPosition());
61 : // check if cosmic ray crossed boundary in last step
62 1 : if (currentDistance * previousDistance < 0){
63 1 : Vector3d currentDirection = c->current.getDirection();
64 1 : Vector3d previousPosition = c->previous.getPosition();
65 : // get point where trajectory intersects the shell boundary
66 1 : double p_half = previousPosition.dot(currentDirection) / currentDirection.getR2();
67 1 : double q = (previousPosition.getR2() - radius * radius) / currentDirection.getR2();
68 1 : double k1 = - p_half + sqrt(p_half * p_half - q);
69 : Vector3d intersectPoint = previousPosition + k1 * currentDirection;
70 : // flip component of velocity normal to surface
71 1 : Vector3d surfaceNormal = normal(intersectPoint);
72 1 : Vector3d newDirection = currentDirection - 2 * currentDirection.dot(surfaceNormal) * surfaceNormal;
73 1 : c->current.setDirection(newDirection.getUnitVector());
74 : // also update position
75 1 : c->current.setPosition(intersectPoint + newDirection * (c->getCurrentStep() - (k1 * currentDirection).getR()) / newDirection.getR());
76 : }
77 1 : }
78 :
79 0 : void ReflectiveShell::setCenter(Vector3d c) {
80 : center = c;
81 0 : }
82 0 : void ReflectiveShell::setRadius(double r) {
83 0 : radius = r;
84 0 : }
85 :
86 0 : std::string ReflectiveShell::getDescription() const {
87 0 : std::stringstream s;
88 0 : s << "Reflective sphere: center: " << center / Mpc << " Mpc, ";
89 0 : s << "radius: " << radius / Mpc << " Mpc";
90 0 : return s.str();
91 0 : };
92 :
93 :
94 0 : ReflectiveBox::ReflectiveBox() :
95 0 : origin(Vector3d(0, 0, 0)), size(Vector3d(0, 0, 0)) {
96 0 : }
97 :
98 1 : ReflectiveBox::ReflectiveBox(Vector3d o, Vector3d s) :
99 1 : origin(o), size(s) {
100 1 : }
101 :
102 1 : void ReflectiveBox::process(Candidate *c) const {
103 1 : Vector3d cur = (c->current.getPosition() - origin) / size; // current position in cell units
104 1 : Vector3d n = cur.floor();
105 :
106 1 : if ((n.x == 0) and (n.y == 0) and (n.z == 0))
107 : return; // do nothing if candidate is inside the box
108 :
109 : // flip direction
110 1 : Vector3d nReflect(pow(-1, n.x), pow(-1, n.y), pow(-1, n.z));
111 1 : c->current.setDirection(c->current.getDirection() * nReflect);
112 1 : c->previous.setDirection(c->previous.getDirection() * nReflect);
113 1 : c->created.setDirection(c->created.getDirection() * nReflect);
114 1 : c->source.setDirection(c->source.getDirection() * nReflect);
115 :
116 1 : Vector3d src = (c->source.getPosition() - origin) / size; // initial position in cell units
117 1 : Vector3d cre = (c->created.getPosition() - origin) / size; // initial position in cell units
118 1 : Vector3d prv = (c->previous.getPosition() - origin) / size; // previous position in cell units
119 :
120 : // repeatedly translate until the current position is inside the cell
121 1 : while ((cur.x < 0) or (cur.x > 1)) {
122 0 : double t = 2 * (cur.x > 1);
123 0 : src.x = t - src.x;
124 0 : cre.x = t - cre.x;
125 0 : prv.x = t - prv.x;
126 0 : cur.x = t - cur.x;
127 : }
128 1 : while ((cur.y < 0) or (cur.y > 1)) {
129 0 : double t = 2 * (cur.y > 1);
130 0 : src.y = t - src.y;
131 0 : cre.y = t - cre.y;
132 0 : prv.y = t - prv.y;
133 0 : cur.y = t - cur.y;
134 : }
135 2 : while ((cur.z < 0) or (cur.z > 1)) {
136 1 : double t = 2 * (cur.z > 1);
137 1 : src.z = t - src.z;
138 1 : cre.z = t - cre.z;
139 1 : prv.z = t - prv.z;
140 1 : cur.z = t - cur.z;
141 : }
142 :
143 1 : c->current.setPosition(cur * size + origin);
144 1 : c->source.setPosition(src * size + origin);
145 1 : c->created.setPosition(cre * size + origin);
146 1 : c->previous.setPosition(prv * size + origin);
147 : }
148 :
149 0 : void ReflectiveBox::setOrigin(Vector3d o) {
150 : origin = o;
151 0 : }
152 0 : void ReflectiveBox::setSize(Vector3d s) {
153 : size = s;
154 0 : }
155 :
156 0 : std::string ReflectiveBox::getDescription() const {
157 0 : std::stringstream s;
158 0 : s << "Reflective box: origin " << origin / Mpc << " Mpc, ";
159 0 : s << "size " << size / Mpc << " Mpc";
160 0 : return s.str();
161 0 : }
162 :
163 0 : CubicBoundary::CubicBoundary() :
164 0 : origin(Vector3d(0, 0, 0)), size(0), limitStep(true), margin(0.1 * kpc) {
165 0 : }
166 :
167 4 : CubicBoundary::CubicBoundary(Vector3d o, double s) :
168 4 : origin(o), size(s), limitStep(true), margin(0.1 * kpc) {
169 4 : }
170 :
171 4 : void CubicBoundary::process(Candidate *c) const {
172 4 : Vector3d r = c->current.getPosition() - origin;
173 : double lo = r.min();
174 : double hi = r.max();
175 4 : if ((lo <= 0) or (hi >= size)) {
176 1 : reject(c);
177 : }
178 4 : if (limitStep) {
179 4 : c->limitNextStep(lo + margin);
180 4 : c->limitNextStep(size - hi + margin);
181 : }
182 4 : }
183 :
184 0 : void CubicBoundary::setOrigin(Vector3d o) {
185 : origin = o;
186 0 : }
187 0 : void CubicBoundary::setSize(double s) {
188 0 : size = s;
189 0 : }
190 2 : void CubicBoundary::setMargin(double m) {
191 2 : margin = m;
192 2 : }
193 2 : void CubicBoundary::setLimitStep(bool b) {
194 2 : limitStep = b;
195 2 : }
196 :
197 0 : std::string CubicBoundary::getDescription() const {
198 0 : std::stringstream s;
199 0 : s << "Cubic Boundary: origin " << origin / Mpc << " Mpc, ";
200 0 : s << "size " << size / Mpc << " Mpc, ";
201 0 : s << "Flag: '" << rejectFlagKey << "' -> '" << rejectFlagValue << "', ";
202 0 : s << "MakeInactive: " << (makeRejectedInactive ? "yes" : "no");
203 0 : if (rejectAction.valid())
204 0 : s << ", Action: " << rejectAction->getDescription();
205 0 : return s.str();
206 0 : }
207 :
208 0 : SphericalBoundary::SphericalBoundary() :
209 0 : center(Vector3d(0, 0, 0)), radius(0), limitStep(true), margin(0.1 * kpc) {
210 0 : }
211 :
212 3 : SphericalBoundary::SphericalBoundary(Vector3d c, double r) :
213 3 : center(c), radius(r), limitStep(true), margin(0.1 * kpc) {
214 3 : }
215 :
216 3 : void SphericalBoundary::process(Candidate *c) const {
217 3 : double d = (c->current.getPosition() - center).getR();
218 3 : if (d >= radius) {
219 1 : reject(c);
220 : }
221 3 : if (limitStep)
222 3 : c->limitNextStep(radius - d + margin);
223 3 : }
224 :
225 0 : void SphericalBoundary::setCenter(Vector3d c) {
226 : center = c;
227 0 : }
228 0 : void SphericalBoundary::setRadius(double r) {
229 0 : radius = r;
230 0 : }
231 1 : void SphericalBoundary::setMargin(double m) {
232 1 : margin = m;
233 1 : }
234 1 : void SphericalBoundary::setLimitStep(bool b) {
235 1 : limitStep = b;
236 1 : }
237 :
238 0 : std::string SphericalBoundary::getDescription() const {
239 0 : std::stringstream s;
240 0 : s << "Spherical Boundary: radius " << radius / Mpc << " Mpc, ";
241 0 : s << "around " << center / Mpc << " Mpc, ";
242 0 : s << "Flag: '" << rejectFlagKey << "' -> '" << rejectFlagValue << "', ";
243 0 : s << "MakeInactive: " << (makeRejectedInactive ? "yes" : "no");
244 0 : if (rejectAction.valid())
245 0 : s << ", Action: " << rejectAction->getDescription();
246 0 : return s.str();
247 0 : }
248 :
249 0 : EllipsoidalBoundary::EllipsoidalBoundary() :
250 : focalPoint1(Vector3d(0, 0, 0)), focalPoint2(Vector3d(0, 0, 0)),
251 0 : majorAxis(0), limitStep(true), margin(0.1 * kpc) {
252 0 : }
253 :
254 3 : EllipsoidalBoundary::EllipsoidalBoundary(Vector3d f1, Vector3d f2, double a) :
255 3 : focalPoint1(f1), focalPoint2(f2), majorAxis(a), limitStep(true),
256 3 : margin(0.1 * kpc) {
257 3 : }
258 :
259 3 : void EllipsoidalBoundary::process(Candidate *c) const {
260 3 : Vector3d pos = c->current.getPosition();
261 3 : double d = pos.getDistanceTo(focalPoint1) + pos.getDistanceTo(focalPoint2);
262 3 : if (d >= majorAxis) {
263 1 : reject(c);
264 : }
265 3 : if (limitStep)
266 3 : c->limitNextStep(majorAxis - d + margin);
267 3 : }
268 :
269 0 : void EllipsoidalBoundary::setFocalPoints(Vector3d f1, Vector3d f2) {
270 : focalPoint1 = f1;
271 : focalPoint2 = f2;
272 0 : }
273 0 : void EllipsoidalBoundary::setMajorAxis(double a) {
274 0 : majorAxis = a;
275 0 : }
276 1 : void EllipsoidalBoundary::setMargin(double m) {
277 1 : margin = m;
278 1 : }
279 1 : void EllipsoidalBoundary::setLimitStep(bool b) {
280 1 : limitStep = b;
281 1 : }
282 :
283 0 : std::string EllipsoidalBoundary::getDescription() const {
284 0 : std::stringstream s;
285 0 : s << "Ellipsoidal Boundary: F1 = " << focalPoint1 / Mpc << " Mpc, ";
286 0 : s << "F2 = " << focalPoint2 / Mpc << " Mpc, ";
287 0 : s << "major axis " << majorAxis / Mpc << " Mpc; ";
288 0 : s << "Flag: '" << rejectFlagKey << "' -> '" << rejectFlagValue << "', ";
289 0 : s << "MakeInactive: " << (makeRejectedInactive ? "yes" : "no");
290 0 : if (rejectAction.valid())
291 0 : s << ", Action: " << rejectAction->getDescription();
292 0 : return s.str();
293 0 : }
294 :
295 0 : CylindricalBoundary::CylindricalBoundary() :
296 0 : origin(Vector3d(0,0,0)), height(0), radius(0), limitStep(false), margin(0) {
297 0 : }
298 :
299 3 : CylindricalBoundary::CylindricalBoundary(Vector3d o, double h, double r) :
300 3 : origin(o), height(h), radius(r), limitStep(false) , margin(0){
301 3 : }
302 :
303 3 : void CylindricalBoundary::process(Candidate *c) const {
304 3 : Vector3d d = c->current.getPosition() - origin;
305 3 : double R2 = pow(d.x, 2.)+pow(d.y, 2.);
306 3 : double Z = fabs(d.z);
307 3 : if ( R2 < pow(radius, 2.) and Z < height/2.) {
308 2 : if(limitStep) {
309 2 : c->limitNextStep(std::min(radius - pow(R2, 0.5), height/2. - Z) + margin);
310 : }
311 : return;
312 : }
313 1 : reject(c);
314 : }
315 :
316 0 : void CylindricalBoundary::setOrigin(Vector3d o) {
317 : origin = o;
318 0 : }
319 0 : void CylindricalBoundary::setHeight(double h) {
320 0 : height = h;
321 0 : }
322 0 : void CylindricalBoundary::setRadius(double r) {
323 0 : radius = r;
324 0 : }
325 1 : void CylindricalBoundary::setLimitStep(bool b) {
326 1 : limitStep = b;
327 1 : }
328 1 : void CylindricalBoundary::setMargin(double m) {
329 1 : margin = m;
330 1 : }
331 :
332 :
333 0 : std::string CylindricalBoundary::getDescription() const {
334 0 : std::stringstream s;
335 0 : s << "Cylindrical Boundary: origin = " << origin / kpc << " kpc, ";
336 0 : s << "radius = " << radius / kpc << " kpc, ";
337 0 : s << "height" << height / kpc << " kpc; ";
338 0 : s << "Flag: '" << rejectFlagKey << "' -> '" << rejectFlagValue << "', ";
339 0 : s << "MakeInactive: " << (makeRejectedInactive ? "yes" : "no");
340 0 : if (rejectAction.valid())
341 0 : s << ", Action: " << rejectAction->getDescription();
342 0 : return s.str();
343 0 : }
344 :
345 : } // namespace crpropa
|
