Line data Source code
1 : #include "crpropa/advectionField/AdvectionField.h"
2 :
3 :
4 : namespace crpropa {
5 :
6 :
7 :
8 3 : void AdvectionFieldList::addField(ref_ptr<AdvectionField> field) {
9 3 : fields.push_back(field);
10 3 : }
11 :
12 1 : Vector3d AdvectionFieldList::getField(const Vector3d &position) const {
13 : Vector3d b(0.);
14 4 : for (int i = 0; i < fields.size(); i++)
15 3 : b += fields[i]->getField(position);
16 1 : return b;
17 : }
18 :
19 1 : double AdvectionFieldList::getDivergence(const Vector3d &position) const {
20 : double D=0.;
21 : // Work on default values for divergence or an error handling
22 4 : for (int i = 0; i < fields.size(); i++)
23 3 : D += fields[i]->getDivergence(position);
24 1 : return D;
25 : }
26 :
27 :
28 : //----------------------------------------------------------------
29 8 : UniformAdvectionField::UniformAdvectionField(const Vector3d &value) :
30 8 : value(value) {
31 8 : }
32 :
33 120005 : Vector3d UniformAdvectionField::getField(const Vector3d &position) const {
34 120005 : return value;
35 : }
36 :
37 5 : double UniformAdvectionField::getDivergence(const Vector3d &position) const {
38 5 : return 0.;
39 : }
40 :
41 0 : std::string UniformAdvectionField::getDescription() const {
42 0 : std::stringstream s;
43 0 : s << "v0: " << value / km * sec << " km/s, ";
44 0 : return s.str();
45 0 : }
46 :
47 : //----------------------------------------------------------------
48 :
49 2 : ConstantSphericalAdvectionField::ConstantSphericalAdvectionField(const Vector3d origin, double vWind) {
50 2 : setOrigin(origin);
51 2 : setVWind(vWind);
52 2 : }
53 :
54 1 : Vector3d ConstantSphericalAdvectionField::getField(const Vector3d &position) const {
55 : Vector3d Pos = position-origin;
56 1 : return vWind * Pos.getUnitVector();
57 : }
58 :
59 2 : double ConstantSphericalAdvectionField::getDivergence(const Vector3d &position) const {
60 : double R = (position-origin).getR();
61 2 : return 2*vWind/R;
62 : }
63 :
64 2 : void ConstantSphericalAdvectionField::setOrigin(const Vector3d o) {
65 : origin=o;
66 2 : return;
67 : }
68 :
69 2 : void ConstantSphericalAdvectionField::setVWind(double v) {
70 2 : vWind = v;
71 2 : return;
72 : }
73 :
74 3 : Vector3d ConstantSphericalAdvectionField::getOrigin() const {
75 3 : return origin;
76 : }
77 :
78 1 : double ConstantSphericalAdvectionField::getVWind() const {
79 1 : return vWind;
80 : }
81 :
82 0 : std::string ConstantSphericalAdvectionField::getDescription() const {
83 0 : std::stringstream s;
84 0 : s << "Origin: " << origin / kpc << " kpc, ";
85 0 : s << "v0: " << vWind / km * sec << " km/s, ";
86 0 : return s.str();
87 0 : }
88 :
89 :
90 :
91 : //----------------------------------------------------------------
92 :
93 1 : SphericalAdvectionField::SphericalAdvectionField(const Vector3d origin, double radius, double vMax, double tau, double alpha) {
94 1 : setOrigin(origin);
95 1 : setRadius(radius);
96 1 : setVMax(vMax);
97 1 : setTau(tau);
98 1 : setAlpha(alpha);
99 1 : }
100 :
101 3 : Vector3d SphericalAdvectionField::getField(const Vector3d &position) const {
102 : Vector3d Pos = position-origin;
103 3 : double R = Pos.getR();
104 3 : if (R>radius) {
105 : return Vector3d(0.);
106 : }
107 2 : double v_R = getV(R);
108 2 : return v_R * Pos.getUnitVector();
109 : }
110 :
111 3 : double SphericalAdvectionField::getDivergence(const Vector3d &position) const {
112 : double R = (position-origin).getR();
113 3 : if (R>radius) {
114 : return 0.;
115 : }
116 2 : double D = 2*vMax/R * ( 1-( 1-alpha*(pow(R, alpha)/(2*tau)) )*exp(-( pow(R, alpha)/tau )) );
117 2 : return D;
118 : }
119 :
120 2 : double SphericalAdvectionField::getV(const double &r) const {
121 2 : double f = vMax * (1-exp(-(pow(r, alpha)/tau)));
122 2 : return f;
123 : }
124 :
125 1 : void SphericalAdvectionField::setOrigin(const Vector3d o) {
126 : origin = o;
127 1 : return;
128 : }
129 :
130 1 : void SphericalAdvectionField::setRadius(double r) {
131 1 : radius = r;
132 1 : return;
133 : }
134 :
135 1 : void SphericalAdvectionField::setVMax(double v) {
136 1 : vMax = v;
137 1 : return;
138 : }
139 :
140 1 : void SphericalAdvectionField::setTau(double t) {
141 1 : tau = t;
142 1 : return;
143 : }
144 :
145 1 : void SphericalAdvectionField::setAlpha(double a) {
146 1 : alpha = a;
147 1 : return;
148 : }
149 :
150 3 : Vector3d SphericalAdvectionField::getOrigin() const {
151 3 : return origin;
152 : }
153 :
154 1 : double SphericalAdvectionField::getRadius() const {
155 1 : return radius;
156 : }
157 :
158 1 : double SphericalAdvectionField::getVMax() const {
159 1 : return vMax;
160 : }
161 :
162 1 : double SphericalAdvectionField::getTau() const {
163 1 : return tau;
164 : }
165 :
166 1 : double SphericalAdvectionField::getAlpha() const {
167 1 : return alpha;
168 : }
169 :
170 0 : std::string SphericalAdvectionField::getDescription() const {
171 0 : std::stringstream s;
172 0 : s << "Origin: " << origin / kpc << " kpc, ";
173 0 : s << "Radius: " << radius / kpc << " kpc, ";
174 0 : s << "vMax: " << vMax / km * sec << " km/s, ";
175 0 : s << "tau: " << tau << ", ";
176 0 : s << "alpha: " << alpha << "\n";
177 0 : return s.str();
178 0 : }
179 :
180 : //----------------------------------------------------------------
181 :
182 0 : OneDimensionalCartesianShock::OneDimensionalCartesianShock(double compressionRatio, double vUp, double lShock){
183 0 : setComp(compressionRatio);
184 0 : setVup(vUp);
185 0 : setShockwidth(lShock);
186 0 : }
187 :
188 0 : Vector3d OneDimensionalCartesianShock::getField(const Vector3d &position) const {
189 0 : double x = position.x;
190 0 : double vDown = vUp / compressionRatio;
191 :
192 0 : double a = (vUp + vDown) * 0.5;
193 0 : double b = (vUp - vDown) * 0.5;
194 :
195 : Vector3d v(0.);
196 0 : v.x = a - b * tanh(x / lShock);
197 0 : return v;
198 :
199 : }
200 :
201 0 : double OneDimensionalCartesianShock::getDivergence(const Vector3d &position) const {
202 0 : double x = position.x;
203 0 : double vDown = vUp / compressionRatio;
204 :
205 : double a = (vUp + vDown) * 0.5;
206 0 : double b = (vUp - vDown) * 0.5;
207 0 : return -b / lShock * (1 - tanh(x / lShock) * tanh(x / lShock));
208 : }
209 :
210 0 : void OneDimensionalCartesianShock::setComp(double r) {
211 0 : compressionRatio = r;
212 0 : return;
213 : }
214 :
215 0 : void OneDimensionalCartesianShock::setVup(double v) {
216 0 : vUp = v;
217 0 : return;
218 : }
219 0 : void OneDimensionalCartesianShock::setShockwidth(double w) {
220 0 : lShock = w;
221 0 : return;
222 : }
223 :
224 0 : double OneDimensionalCartesianShock::getComp() const {
225 0 : return compressionRatio;
226 : }
227 :
228 0 : double OneDimensionalCartesianShock::getVup() const {
229 0 : return vUp;
230 : }
231 :
232 0 : double OneDimensionalCartesianShock::getShockwidth() const {
233 0 : return lShock;
234 : }
235 :
236 0 : std::string OneDimensionalCartesianShock::getDescription() const {
237 0 : std::stringstream s;
238 0 : s << "Shock width: " << lShock / km << " km, ";
239 0 : s << "Vup: " << vUp / km * sec << " km/s, ";
240 0 : s << "Compression: " << compressionRatio;
241 0 : return s.str();
242 0 : }
243 :
244 : //----------------------------------------------------------------
245 :
246 0 : OneDimensionalSphericalShock::OneDimensionalSphericalShock(double rShock, double vUp, double compressionRatio, double lShock, bool coolUpstream ){
247 0 : setComp(compressionRatio);
248 0 : setVup(vUp);
249 0 : setShockwidth(lShock);
250 0 : setShockRadius(rShock);
251 0 : setCooling(coolUpstream);
252 0 : }
253 :
254 0 : Vector3d OneDimensionalSphericalShock::getField(const Vector3d &position) const {
255 : double r = position.getR();
256 0 : Vector3d e_r = position.getUnitVector();
257 :
258 0 : double vDown = vUp / compressionRatio;
259 0 : double a = (vUp + vDown) * 0.5;
260 0 : double b = (vUp - vDown) * 0.5;
261 :
262 : double v;
263 0 : if (coolUpstream == true){
264 :
265 0 : if (r <= rShock)
266 0 : v = a - b * tanh((r-rShock) / lShock);
267 : else
268 0 : v = (a - b * tanh((r-rShock) / lShock)) * (rShock / r) * (rShock / r);
269 :
270 : }
271 : else
272 0 : v = (a - b * tanh((r-rShock) / lShock)) * (rShock / r) * (rShock / r);
273 :
274 0 : return v * e_r;
275 : }
276 :
277 0 : double OneDimensionalSphericalShock::getDivergence(const Vector3d &position) const {
278 : double r = position.getR();
279 :
280 0 : double vDown = vUp / compressionRatio;
281 0 : double a = (vUp + vDown) * 0.5;
282 0 : double b = (vUp - vDown) * 0.5;
283 :
284 0 : double c = tanh((r-rShock) / lShock);
285 :
286 0 : if (coolUpstream == true){
287 0 : if (r <= rShock)
288 0 : return 2 * a / r - 2 * b / r * c - b / lShock * (1 - c * c);
289 : else
290 0 : return -(rShock / r) * (rShock / r) * b / lShock * (1 - c * c);
291 : }
292 : else
293 0 : return -(rShock / r) * (rShock / r) * b / lShock * (1 - c * c);
294 :
295 : }
296 :
297 0 : void OneDimensionalSphericalShock::setComp(double r) {
298 0 : compressionRatio = r;
299 0 : return;
300 : }
301 :
302 0 : void OneDimensionalSphericalShock::setVup(double v) {
303 0 : vUp = v;
304 0 : return;
305 : }
306 0 : void OneDimensionalSphericalShock::setShockwidth(double w) {
307 0 : lShock = w;
308 0 : return;
309 : }
310 :
311 0 : void OneDimensionalSphericalShock::setShockRadius(double r) {
312 0 : rShock = r;
313 0 : return;
314 : }
315 :
316 0 : void OneDimensionalSphericalShock::setCooling(bool c) {
317 0 : coolUpstream = c;
318 0 : return;
319 : }
320 :
321 0 : double OneDimensionalSphericalShock::getComp() const {
322 0 : return compressionRatio;
323 : }
324 :
325 0 : double OneDimensionalSphericalShock::getVup() const {
326 0 : return vUp;
327 : }
328 :
329 0 : double OneDimensionalSphericalShock::getShockwidth() const {
330 0 : return lShock;
331 : }
332 :
333 0 : double OneDimensionalSphericalShock::getShockRadius() const {
334 0 : return rShock;
335 : }
336 :
337 0 : bool OneDimensionalSphericalShock::getCooling() const {
338 0 : return coolUpstream;
339 : }
340 :
341 0 : std::string OneDimensionalSphericalShock::getDescription() const {
342 0 : std::stringstream s;
343 0 : s << "Shock width: " << lShock / km << " km, ";
344 0 : s << "Shock radius: " << rShock / km << " km, ";
345 0 : s << "Vup: " << vUp / km * sec << " km/s, ";
346 0 : s << "Comp: " << compressionRatio;
347 :
348 0 : return s.str();
349 0 : }
350 :
351 : //----------------------------------------------------------------
352 :
353 0 : ObliqueAdvectionShock::ObliqueAdvectionShock(double compressionRatio, double vXUp, double vY, double lShock) {
354 0 : setComp(compressionRatio);
355 0 : setVup(vXUp);
356 0 : setVy(vY);
357 0 : setShockwidth(lShock);
358 0 : }
359 :
360 0 : Vector3d ObliqueAdvectionShock::getField(const Vector3d &position) const {
361 0 : double x = position.x;
362 0 : double vXDown = vXUp / compressionRatio;
363 :
364 0 : double a = (vXUp + vXDown) * 0.5;
365 0 : double b = (vXUp - vXDown) * 0.5;
366 :
367 : Vector3d v(0.);
368 0 : v.x = a - b * tanh(x / lShock);
369 0 : v.y = vY;
370 :
371 0 : return v;
372 : }
373 :
374 0 : double ObliqueAdvectionShock::getDivergence(const Vector3d &position) const {
375 0 : double x = position.x;
376 0 : double vXDown = vXUp / compressionRatio;
377 : // vy = const
378 :
379 : double a = (vXUp + vXDown) * 0.5;
380 0 : double b = (vXUp - vXDown) * 0.5;
381 :
382 0 : return -b / lShock * (1 - tanh(x / lShock) * tanh(x / lShock));
383 : }
384 :
385 0 : void ObliqueAdvectionShock::setComp(double r) {
386 0 : compressionRatio = r;
387 0 : return;
388 : }
389 :
390 0 : void ObliqueAdvectionShock::setVup(double v) {
391 0 : vXUp = v;
392 0 : return;
393 : }
394 :
395 0 : void ObliqueAdvectionShock::setVy(double v) {
396 0 : vY = v;
397 0 : return;
398 : }
399 0 : void ObliqueAdvectionShock::setShockwidth(double w) {
400 0 : lShock = w;
401 0 : return;
402 : }
403 :
404 0 : double ObliqueAdvectionShock::getComp() const {
405 0 : return compressionRatio;
406 : }
407 :
408 0 : double ObliqueAdvectionShock::getVup() const {
409 0 : return vXUp;
410 : }
411 :
412 0 : double ObliqueAdvectionShock::getVy() const {
413 0 : return vY;
414 : }
415 :
416 0 : double ObliqueAdvectionShock::getShockwidth() const {
417 0 : return lShock;
418 : }
419 :
420 0 : std::string ObliqueAdvectionShock::getDescription() const {
421 0 : std::stringstream s;
422 0 : s << "Shock width: " << lShock / km << " km, ";
423 0 : s << "Vx_up: " << vXUp / km * sec << " km/s, ";
424 0 : s << "Vy: " << vY / km * sec << " km/s, ";
425 0 : s << "Comp: " << compressionRatio;
426 :
427 0 : return s.str();
428 0 : }
429 :
430 : //-----------------------------------------------------------------
431 :
432 1 : SphericalAdvectionShock::SphericalAdvectionShock(const Vector3d origin, double r_0, double v_0, double l) {
433 1 : setOrigin(origin);
434 1 : setR0(r_0);
435 1 : setV0(v_0);
436 1 : setLambda(l);
437 1 : setRRot(r_0);
438 1 : setAzimuthalSpeed(0.);
439 1 : }
440 :
441 :
442 3 : Vector3d SphericalAdvectionShock::getField(const Vector3d &pos) const {
443 : Vector3d R = pos-origin;
444 3 : Vector3d e_r = R.getUnitVector();
445 3 : Vector3d e_phi = R.getUnitVectorPhi();
446 : double r = R.getR();
447 :
448 3 : double v_r = v_0 * ( 1 + (pow(r_0/(2*r), 2.) -1 ) * g(r));
449 3 : double v_p = v_phi * (r_rot/r);
450 :
451 3 : return v_r * e_r + v_p * e_phi;
452 : }
453 :
454 :
455 2 : double SphericalAdvectionShock::getDivergence(const Vector3d &pos) const {
456 : double r = (pos-origin).getR();
457 :
458 2 : double d1 = 2./r*(1-g(r));
459 2 : double d2 = (pow(r_0/(2*r), 2.)-1)*g_prime(r);
460 :
461 2 : return v_0 * (d1+d2);
462 : }
463 :
464 :
465 5 : double SphericalAdvectionShock::g(double r) const {
466 5 : double a = (r-r_0)/lambda;
467 5 : return 1. / (1+exp(-a));
468 : }
469 :
470 2 : double SphericalAdvectionShock::g_prime(double r) const {
471 2 : double a = (r-r_0)/lambda;
472 2 : return 1. / (2*lambda*(1+cosh(-a)));
473 : }
474 :
475 1 : void SphericalAdvectionShock::setOrigin(const Vector3d o) {
476 : origin = o;
477 1 : }
478 :
479 1 : void SphericalAdvectionShock::setR0(double r) {
480 1 : r_0 = r;
481 1 : }
482 :
483 1 : void SphericalAdvectionShock::setV0(double v) {
484 1 : v_0 = v;
485 1 : }
486 :
487 1 : void SphericalAdvectionShock::setLambda(double l) {
488 1 : lambda = l;
489 1 : }
490 :
491 2 : void SphericalAdvectionShock::setRRot(double r) {
492 2 : r_rot = r;
493 2 : }
494 :
495 2 : void SphericalAdvectionShock::setAzimuthalSpeed(double v) {
496 2 : v_phi = v;
497 2 : }
498 :
499 3 : Vector3d SphericalAdvectionShock::getOrigin() const {
500 3 : return origin;
501 : }
502 :
503 1 : double SphericalAdvectionShock::getR0() const {
504 1 : return r_0;
505 : }
506 :
507 1 : double SphericalAdvectionShock::getV0() const {
508 1 : return v_0;
509 : }
510 :
511 1 : double SphericalAdvectionShock::getLambda() const {
512 1 : return lambda;
513 : }
514 :
515 2 : double SphericalAdvectionShock::getRRot() const {
516 2 : return r_rot;
517 : }
518 :
519 2 : double SphericalAdvectionShock::getAzimuthalSpeed() const {
520 2 : return v_phi;
521 : }
522 :
523 0 : std::string SphericalAdvectionShock::getDescription() const {
524 0 : std::stringstream s;
525 0 : s << "Origin: " << origin / kpc << " kpc, ";
526 0 : s << "r0 (shock radius): " << r_0 / kpc << " kpc, ";
527 0 : s << "r_rot (norm. azimuthal velocity): " << r_rot / kpc << " kpc, ";
528 0 : s << "v0 (maximum radial speed): " << v_0 / km * sec << " km/s, ";
529 0 : s << "v_phi (azimuthal speed @ r_rot): " << v_phi / km * sec << " km/s, ";
530 0 : s << "lambda: " << lambda / pc << " pc";
531 0 : return s.str();
532 0 : }
533 :
534 : } // namespace crpropa
|
