Line data Source code
1 : #include "crpropa/Source.h"
2 : #include "crpropa/Random.h"
3 : #include "crpropa/Cosmology.h"
4 : #include "crpropa/Common.h"
5 : #include "crpropa/Units.h"
6 : #include "crpropa/ParticleID.h"
7 :
8 : #ifdef CRPROPA_HAVE_MUPARSER
9 : #include "muParser.h"
10 : #endif
11 :
12 : #include <sstream>
13 : #include <stdexcept>
14 :
15 : namespace crpropa {
16 :
17 : // Source ---------------------------------------------------------------------
18 19 : void Source::add(SourceFeature* property) {
19 19 : features.push_back(property);
20 19 : }
21 :
22 2103 : ref_ptr<Candidate> Source::getCandidate() const {
23 2103 : ref_ptr<Candidate> candidate = new Candidate();
24 7512 : for (int i = 0; i < features.size(); i++)
25 5409 : (*features[i]).prepareCandidate(*candidate);
26 2103 : return candidate;
27 : }
28 :
29 0 : std::string Source::getDescription() const {
30 0 : std::stringstream ss;
31 0 : ss << "Cosmic ray source\n";
32 0 : for (int i = 0; i < features.size(); i++)
33 0 : ss << " " << features[i]->getDescription();
34 0 : return ss.str();
35 0 : }
36 :
37 : // SourceList------------------------------------------------------------------
38 3 : void SourceList::add(Source* source, double weight) {
39 6 : sources.push_back(source);
40 3 : if (cdf.size() > 0)
41 1 : weight += cdf.back();
42 3 : cdf.push_back(weight);
43 3 : }
44 :
45 1002 : ref_ptr<Candidate> SourceList::getCandidate() const {
46 1002 : if (sources.size() == 0)
47 1 : throw std::runtime_error("SourceList: no sources set");
48 1001 : size_t i = Random::instance().randBin(cdf);
49 1001 : return (sources[i])->getCandidate();
50 : }
51 :
52 0 : std::string SourceList::getDescription() const {
53 0 : std::stringstream ss;
54 0 : ss << "List of cosmic ray sources\n";
55 0 : for (int i = 0; i < sources.size(); i++)
56 0 : ss << " " << sources[i]->getDescription();
57 0 : return ss.str();
58 0 : }
59 :
60 : // SourceFeature---------------------------------------------------------------
61 5407 : void SourceFeature::prepareCandidate(Candidate& candidate) const {
62 5407 : ParticleState &source = candidate.source;
63 5407 : prepareParticle(source);
64 : candidate.created = source;
65 : candidate.current = source;
66 : candidate.previous = source;
67 5407 : }
68 :
69 0 : std::string SourceFeature::getDescription() const {
70 0 : return description;
71 : }
72 :
73 : // ----------------------------------------------------------------------------
74 3 : SourceParticleType::SourceParticleType(int id) :
75 3 : id(id) {
76 3 : setDescription();
77 3 : }
78 :
79 1101 : void SourceParticleType::prepareParticle(ParticleState& particle) const {
80 1101 : particle.setId(id);
81 1101 : }
82 :
83 3 : void SourceParticleType::setDescription() {
84 3 : std::stringstream ss;
85 3 : ss << "SourceParticleType: " << id << "\n";
86 3 : description = ss.str();
87 3 : }
88 :
89 : // ----------------------------------------------------------------------------
90 0 : SourceMultipleParticleTypes::SourceMultipleParticleTypes() {
91 0 : setDescription();
92 0 : }
93 :
94 0 : void SourceMultipleParticleTypes::add(int id, double a) {
95 0 : particleTypes.push_back(id);
96 0 : if (cdf.size() > 0)
97 0 : a += cdf.back();
98 0 : cdf.push_back(a);
99 0 : setDescription();
100 0 : }
101 :
102 0 : void SourceMultipleParticleTypes::prepareParticle(ParticleState& particle) const {
103 0 : if (particleTypes.size() == 0)
104 0 : throw std::runtime_error("SourceMultipleParticleTypes: no nuclei set");
105 0 : size_t i = Random::instance().randBin(cdf);
106 0 : particle.setId(particleTypes[i]);
107 0 : }
108 :
109 0 : void SourceMultipleParticleTypes::setDescription() {
110 0 : std::stringstream ss;
111 0 : ss << "SourceMultipleParticleTypes: Random particle type\n";
112 0 : for (int i = 0; i < particleTypes.size(); i++)
113 0 : ss << " ID = " << particleTypes[i] << "\n";
114 0 : description = ss.str();
115 0 : }
116 :
117 : // ----------------------------------------------------------------------------
118 2 : SourceEnergy::SourceEnergy(double energy) :
119 2 : E(energy) {
120 2 : setDescription();
121 2 : }
122 :
123 1000 : void SourceEnergy::prepareParticle(ParticleState& p) const {
124 1000 : p.setEnergy(E);
125 1000 : }
126 :
127 2 : void SourceEnergy::setDescription() {
128 2 : std::stringstream ss;
129 2 : ss << "SourceEnergy: " << E / EeV << " EeV\n";
130 2 : description = ss.str();
131 2 : }
132 :
133 : // ----------------------------------------------------------------------------
134 4 : SourcePowerLawSpectrum::SourcePowerLawSpectrum(double Emin, double Emax,
135 4 : double index) :
136 4 : Emin(Emin), Emax(Emax), index(index) {
137 4 : setDescription();
138 4 : }
139 :
140 1102 : void SourcePowerLawSpectrum::prepareParticle(ParticleState& particle) const {
141 1102 : Random &random = Random::instance();
142 1102 : double E = random.randPowerLaw(index, Emin, Emax);
143 1102 : particle.setEnergy(E);
144 1102 : }
145 :
146 4 : void SourcePowerLawSpectrum::setDescription() {
147 4 : std::stringstream ss;
148 4 : ss << "SourcePowerLawSpectrum: Random energy ";
149 8 : ss << "E = " << Emin / EeV << " - " << Emax / EeV << " EeV, ";
150 4 : ss << "dN/dE ~ E^" << index << "\n";
151 4 : description = ss.str();
152 4 : }
153 :
154 : // ----------------------------------------------------------------------------
155 3 : SourceComposition::SourceComposition(double Emin, double Rmax, double index) :
156 3 : Emin(Emin), Rmax(Rmax), index(index) {
157 3 : setDescription();
158 3 : }
159 :
160 5 : void SourceComposition::add(int id, double weight) {
161 5 : nuclei.push_back(id);
162 5 : int A = massNumber(id);
163 5 : int Z = chargeNumber(id);
164 :
165 5 : double a = 1 + index;
166 5 : if (std::abs(a) < std::numeric_limits<double>::min())
167 5 : weight *= log(Z * Rmax / Emin);
168 : else
169 0 : weight *= (pow(Z * Rmax, a) - pow(Emin, a)) / a;
170 :
171 5 : weight *= pow(A, -a);
172 :
173 5 : if (cdf.size() > 0)
174 3 : weight += cdf.back();
175 5 : cdf.push_back(weight);
176 5 : setDescription();
177 5 : }
178 :
179 4 : void SourceComposition::add(int A, int Z, double a) {
180 4 : add(nucleusId(A, Z), a);
181 4 : }
182 :
183 3 : void SourceComposition::prepareParticle(ParticleState& particle) const {
184 3 : if (nuclei.size() == 0)
185 1 : throw std::runtime_error("SourceComposition: No source isotope set");
186 :
187 2 : Random &random = Random::instance();
188 :
189 : // draw random particle type
190 2 : size_t i = random.randBin(cdf);
191 2 : int id = nuclei[i];
192 2 : particle.setId(id);
193 :
194 : // random energy from power law
195 2 : int Z = chargeNumber(id);
196 2 : particle.setEnergy(random.randPowerLaw(index, Emin, Z * Rmax));
197 2 : }
198 :
199 8 : void SourceComposition::setDescription() {
200 8 : std::stringstream ss;
201 8 : ss << "SourceComposition: Random element and energy ";
202 16 : ss << "E = " << Emin / EeV << " - Z*" << Rmax / EeV << " EeV, ";
203 8 : ss << "dN/dE ~ E^" << index << "\n";
204 19 : for (int i = 0; i < nuclei.size(); i++)
205 11 : ss << " ID = " << nuclei[i] << "\n";
206 8 : description = ss.str();
207 8 : }
208 :
209 : // ----------------------------------------------------------------------------
210 5 : SourcePosition::SourcePosition(Vector3d position) :
211 5 : position(position) {
212 5 : setDescription();
213 5 : }
214 :
215 0 : SourcePosition::SourcePosition(double d) :
216 0 : position(Vector3d(d, 0, 0)) {
217 0 : setDescription();
218 0 : }
219 :
220 1103 : void SourcePosition::prepareParticle(ParticleState& particle) const {
221 1103 : particle.setPosition(position);
222 1103 : }
223 :
224 5 : void SourcePosition::setDescription() {
225 5 : std::stringstream ss;
226 5 : ss << "SourcePosition: " << position / Mpc << " Mpc\n";
227 5 : description = ss.str();
228 5 : }
229 :
230 : // ----------------------------------------------------------------------------
231 1 : SourceMultiplePositions::SourceMultiplePositions() {
232 1 : setDescription();
233 1 : }
234 :
235 2 : void SourceMultiplePositions::add(Vector3d pos, double weight) {
236 2 : positions.push_back(pos);
237 2 : if (cdf.size() > 0)
238 1 : weight += cdf.back();
239 2 : cdf.push_back(weight);
240 2 : }
241 :
242 10000 : void SourceMultiplePositions::prepareParticle(ParticleState& particle) const {
243 10000 : if (positions.size() == 0)
244 0 : throw std::runtime_error("SourceMultiplePositions: no position set");
245 10000 : size_t i = Random::instance().randBin(cdf);
246 10000 : particle.setPosition(positions[i]);
247 10000 : }
248 :
249 1 : void SourceMultiplePositions::setDescription() {
250 1 : std::stringstream ss;
251 1 : ss << "SourceMultiplePositions: Random position from list\n";
252 1 : for (int i = 0; i < positions.size(); i++)
253 0 : ss << " " << positions[i] / Mpc << " Mpc\n";
254 1 : description = ss.str();
255 1 : }
256 :
257 : // ----------------------------------------------------------------------------
258 1 : SourceUniformSphere::SourceUniformSphere(Vector3d center, double radius) :
259 1 : center(center), radius(radius) {
260 1 : setDescription();
261 1 : }
262 :
263 1 : void SourceUniformSphere::prepareParticle(ParticleState& particle) const {
264 1 : Random &random = Random::instance();
265 1 : double r = pow(random.rand(), 1. / 3.) * radius;
266 1 : particle.setPosition(center + random.randVector() * r);
267 1 : }
268 :
269 1 : void SourceUniformSphere::setDescription() {
270 1 : std::stringstream ss;
271 1 : ss << "SourceUniformSphere: Random position within a sphere at ";
272 1 : ss << center / Mpc << " Mpc with";
273 1 : ss << radius / Mpc << " Mpc radius\n";
274 1 : description = ss.str();
275 1 : }
276 :
277 : // ----------------------------------------------------------------------------
278 1 : SourceUniformHollowSphere::SourceUniformHollowSphere(
279 : Vector3d center,
280 : double radius_inner,
281 1 : double radius_outer) :
282 1 : center(center), radius_inner(radius_inner),
283 1 : radius_outer(radius_outer) {
284 1 : setDescription();
285 1 : }
286 :
287 100 : void SourceUniformHollowSphere::prepareParticle(ParticleState& particle) const {
288 100 : Random &random = Random::instance();
289 100 : double r = radius_inner + pow(random.rand(), 1. / 3.) * (radius_outer - radius_inner);
290 100 : particle.setPosition(center + random.randVector() * r);
291 100 : }
292 :
293 1 : void SourceUniformHollowSphere::setDescription() {
294 1 : std::stringstream ss;
295 1 : ss << "SourceUniformHollowSphere: Random position within a sphere at ";
296 1 : ss << center / Mpc << " Mpc with";
297 1 : ss << radius_inner / Mpc << " Mpc inner radius\n";
298 1 : ss << radius_outer / Mpc << " Mpc outer radius\n";
299 1 : description = ss.str();
300 1 : }
301 :
302 : // ----------------------------------------------------------------------------
303 0 : SourceUniformShell::SourceUniformShell(Vector3d center, double radius) :
304 0 : center(center), radius(radius) {
305 0 : setDescription();
306 0 : }
307 :
308 0 : void SourceUniformShell::prepareParticle(ParticleState& particle) const {
309 0 : Random &random = Random::instance();
310 0 : particle.setPosition(center + random.randVector() * radius);
311 0 : }
312 :
313 0 : void SourceUniformShell::setDescription() {
314 0 : std::stringstream ss;
315 0 : ss << "SourceUniformShell: Random position on a spherical shell at ";
316 0 : ss << center / Mpc << " Mpc with ";
317 0 : ss << radius / Mpc << " Mpc radius\n";
318 0 : description = ss.str();
319 0 : }
320 :
321 : // ----------------------------------------------------------------------------
322 1 : SourceUniformBox::SourceUniformBox(Vector3d origin, Vector3d size) :
323 1 : origin(origin), size(size) {
324 1 : setDescription();
325 1 : }
326 :
327 1 : void SourceUniformBox::prepareParticle(ParticleState& particle) const {
328 1 : Random &random = Random::instance();
329 1 : Vector3d pos(random.rand(), random.rand(), random.rand());
330 1 : particle.setPosition(pos * size + origin);
331 1 : }
332 :
333 1 : void SourceUniformBox::setDescription() {
334 1 : std::stringstream ss;
335 1 : ss << "SourceUniformBox: Random uniform position in box with ";
336 1 : ss << "origin = " << origin / Mpc << " Mpc and ";
337 1 : ss << "size = " << size / Mpc << " Mpc\n";
338 1 : description = ss.str();
339 1 : }
340 :
341 : // ---------------------------------------------------------------------------
342 1 : SourceUniformCylinder::SourceUniformCylinder(Vector3d origin, double height, double radius) :
343 1 : origin(origin), height(height), radius(radius) {
344 1 : }
345 :
346 1 : void SourceUniformCylinder::prepareParticle(ParticleState& particle) const {
347 1 : Random &random = Random::instance();
348 1 : double phi = 2*M_PI*random.rand();
349 1 : double RandRadius = radius*pow(random.rand(), 1. / 2.);
350 1 : Vector3d pos(cos(phi)*RandRadius, sin(phi)*RandRadius, (-0.5+random.rand())*height);
351 1 : particle.setPosition(pos + origin);
352 1 : }
353 :
354 0 : void SourceUniformCylinder::setDescription() {
355 0 : std::stringstream ss;
356 0 : ss << "SourceUniformCylinder: Random uniform position in cylinder with ";
357 0 : ss << "origin = " << origin / Mpc << " Mpc and ";
358 0 : ss << "radius = " << radius / Mpc << " Mpc and";
359 0 : ss << "height = " << height / Mpc << " Mpc\n";
360 0 : description = ss.str();
361 0 : }
362 :
363 : // ---------------------------------------------------------------------------
364 0 : SourceSNRDistribution::SourceSNRDistribution() :
365 0 : rEarth(8.5 * kpc), beta(3.53), zg(0.3 * kpc) {
366 0 : setAlpha(2.);
367 0 : setFrMax();
368 0 : setFzMax(0.3 * kpc);
369 0 : setRMax(20 * kpc);
370 0 : setZMax(5 * kpc);
371 0 : }
372 :
373 1 : SourceSNRDistribution::SourceSNRDistribution(double rEarth, double alpha, double beta, double zg) :
374 1 : rEarth(rEarth), beta(beta), zg(zg) {
375 1 : setAlpha(alpha);
376 1 : setFrMax();
377 1 : setFzMax(zg);
378 1 : setRMax(20 * kpc);
379 1 : setZMax(5 * kpc);
380 1 : }
381 :
382 100001 : void SourceSNRDistribution::prepareParticle(ParticleState& particle) const {
383 100001 : Random &random = Random::instance();
384 : double RPos;
385 : while (true) {
386 193062 : RPos = random.rand() * rMax;
387 193062 : double fTest = random.rand() * frMax;
388 193062 : double fR = fr(RPos);
389 193062 : if (fTest <= fR) {
390 : break;
391 : }
392 : }
393 : double ZPos;
394 : while (true) {
395 1668515 : ZPos = (random.rand() - 0.5) * 2 * zMax;
396 1668515 : double fTest = random.rand() * fzMax;
397 1668515 : double fZ=fz(ZPos);
398 1668515 : if (fTest<=fZ) {
399 : break;
400 : }
401 : }
402 100001 : double phi = random.rand() * 2 * M_PI;
403 100001 : Vector3d pos(cos(phi) * RPos, sin(phi) * RPos, ZPos);
404 100001 : particle.setPosition(pos);
405 100001 : }
406 :
407 193062 : double SourceSNRDistribution::fr(double r) const {
408 193062 : return pow(r / rEarth, alpha) * exp(- beta * (r - rEarth) / rEarth);
409 : }
410 :
411 1668515 : double SourceSNRDistribution::fz(double z) const{
412 : double Az = 1.;
413 1668515 : double f = 1. / zg * exp(- fabs(z) / zg);
414 : double fz = Az * f;
415 1668515 : return fz;
416 : }
417 :
418 0 : double SourceSNRDistribution::getAlpha() const {
419 0 : return alpha - 1; // -1 to account for the R-term in the volume element dV = R * dR * dphi * dz
420 : }
421 :
422 0 : double SourceSNRDistribution::getBeta() const {
423 0 : return beta;
424 : }
425 :
426 2 : void SourceSNRDistribution::setFrMax() {
427 2 : frMax = pow(alpha / beta, alpha) * exp(beta - alpha);
428 2 : return;
429 : }
430 :
431 1 : void SourceSNRDistribution::setFzMax(double zg) {
432 1 : fzMax = 1. / zg;
433 1 : return;
434 : }
435 :
436 2 : void SourceSNRDistribution::setRMax(double r) {
437 2 : rMax = r;
438 2 : return;
439 : }
440 :
441 1 : void SourceSNRDistribution::setZMax(double z) {
442 1 : zMax = z;
443 1 : return;
444 : }
445 :
446 0 : double SourceSNRDistribution::getFrMax() const {
447 0 : return frMax;
448 : }
449 :
450 0 : double SourceSNRDistribution::getFzMax() const {
451 0 : return fzMax;
452 : }
453 :
454 0 : double SourceSNRDistribution::getRMax() const {
455 0 : return rMax;
456 : }
457 :
458 0 : double SourceSNRDistribution::getZMax() const {
459 0 : return zMax;
460 : }
461 :
462 1 : void SourceSNRDistribution::setAlpha(double a) {
463 1 : alpha = a + 1.; // add 1 for dV = r * dR * dphi * dz
464 1 : setRMax(rMax);
465 1 : setFrMax();
466 1 : }
467 :
468 0 : void SourceSNRDistribution::setBeta(double b) {
469 0 : beta = b;
470 0 : setRMax(rMax);
471 0 : setFrMax();
472 0 : }
473 :
474 0 : void SourceSNRDistribution::setDescription() {
475 0 : std::stringstream ss;
476 0 : ss << "SourceSNRDistribution: Random position according to SNR distribution";
477 0 : ss << "rEarth = " << rEarth / kpc << " kpc and ";
478 0 : ss << "zg = " << zg / kpc << " kpc and";
479 0 : ss << "beta = " << beta << " \n";
480 0 : description = ss.str();
481 0 : }
482 :
483 : // ---------------------------------------------------------------------------
484 0 : SourcePulsarDistribution::SourcePulsarDistribution() :
485 0 : rEarth(8.5*kpc), beta(3.53), zg(0.3*kpc) {
486 0 : setFrMax(8.5*kpc, 3.53);
487 0 : setFzMax(0.3*kpc);
488 0 : setRMax(22*kpc);
489 0 : setZMax(5*kpc);
490 0 : setRBlur(0.07);
491 0 : setThetaBlur(0.35/kpc);
492 0 : }
493 :
494 0 : SourcePulsarDistribution::SourcePulsarDistribution(double rEarth, double beta, double zg, double rB, double tB) :
495 0 : rEarth(rEarth), beta(beta), zg(zg) {
496 0 : setFrMax(rEarth, beta);
497 0 : setFzMax(zg);
498 0 : setRBlur(rB);
499 0 : setThetaBlur(tB);
500 0 : setRMax(22 * kpc);
501 0 : setZMax(5 * kpc);
502 0 : }
503 :
504 0 : void SourcePulsarDistribution::prepareParticle(ParticleState& particle) const {
505 0 : Random &random = Random::instance();
506 : double Rtilde;
507 : while (true) {
508 0 : Rtilde = random.rand() * rMax;
509 0 : double fTest = random.rand() * frMax * 1.1;
510 0 : double fR = fr(Rtilde);
511 0 : if (fTest <= fR) {
512 : break;
513 : }
514 : }
515 : double ZPos;
516 : while (true) {
517 0 : ZPos = (random.rand() - 0.5) * 2 * zMax;
518 0 : double fTest = random.rand() * fzMax;
519 0 : double fZ = fz(ZPos);
520 0 : if (fTest <= fZ) {
521 : break;
522 : }
523 : }
524 :
525 0 : int i = random.randInt(3);
526 0 : double thetaTilde = ftheta(i, Rtilde);
527 0 : double RPos = blurR(Rtilde);
528 0 : double phi = blurTheta(thetaTilde, Rtilde);
529 0 : Vector3d pos(cos(phi) * RPos, sin(phi) * RPos, ZPos);
530 :
531 0 : particle.setPosition(pos);
532 0 : }
533 :
534 0 : double SourcePulsarDistribution::fr(double r) const {
535 0 : double f = r * pow(r / rEarth, 2.) * exp(-beta * (r - rEarth) / rEarth);
536 0 : return f;
537 : }
538 :
539 0 : double SourcePulsarDistribution::fz(double z) const{
540 : double Az = 1.;
541 0 : double f = 1. / zg * exp(- fabs(z) / zg);
542 : double fz = Az * f;
543 0 : return fz;
544 : }
545 :
546 0 : double SourcePulsarDistribution::ftheta(int i, double r) const {
547 0 : const double k_0[] = {4.25, 4.25, 4.89, 4.89};
548 0 : const double r_0[] = {3.48 * kpc, 3.48 * kpc, 4.9 * kpc, 4.9 * kpc};
549 0 : const double theta_0[] = {0., 3.14, 2.52, -0.62};
550 0 : double K = k_0[i];
551 0 : double R = r_0[i];
552 0 : double Theta = theta_0[i];
553 :
554 0 : double theta = K * log(r / R) + Theta;
555 :
556 0 : return theta;
557 : }
558 :
559 0 : double SourcePulsarDistribution::blurR(double rTilde) const {
560 0 : Random &random = Random::instance();
561 0 : return random.randNorm(rTilde, rBlur * rTilde);
562 : }
563 :
564 0 : double SourcePulsarDistribution::blurTheta(double thetaTilde, double rTilde) const {
565 0 : Random &random = Random::instance();
566 0 : double thetaCorr = (random.rand() - 0.5) * 2 * M_PI;
567 0 : double tau = thetaCorr * exp(- thetaBlur * rTilde);
568 0 : return thetaTilde + tau;
569 : }
570 :
571 0 : void SourcePulsarDistribution::setFrMax(double R, double b) {
572 0 : double r = 3 * R / b;
573 0 : frMax = fr(r);
574 0 : }
575 :
576 0 : void SourcePulsarDistribution::setFzMax(double zg) {
577 0 : fzMax = 1. / zg;
578 0 : return;
579 : }
580 :
581 0 : void SourcePulsarDistribution::setRMax(double r) {
582 0 : rMax = r;
583 0 : return;
584 : }
585 :
586 0 : void SourcePulsarDistribution::setZMax(double z) {
587 0 : zMax = z;
588 0 : return;
589 : }
590 :
591 0 : void SourcePulsarDistribution::setRBlur(double r) {
592 0 : rBlur = r;
593 0 : return;
594 : }
595 :
596 0 : void SourcePulsarDistribution::setThetaBlur(double theta) {
597 0 : thetaBlur = theta;
598 0 : return;
599 : }
600 :
601 0 : double SourcePulsarDistribution::getFrMax() {
602 0 : return frMax;
603 : }
604 :
605 0 : double SourcePulsarDistribution::getFzMax() {
606 0 : return fzMax;
607 : }
608 :
609 0 : double SourcePulsarDistribution::getRMax() {
610 0 : return rMax;
611 : }
612 :
613 0 : double SourcePulsarDistribution::getZMax() {
614 0 : return zMax;
615 : }
616 :
617 0 : double SourcePulsarDistribution::getRBlur() {
618 0 : return rBlur;
619 : }
620 :
621 0 : double SourcePulsarDistribution::getThetaBlur() {
622 0 : return thetaBlur;
623 : }
624 :
625 :
626 0 : void SourcePulsarDistribution::setDescription() {
627 0 : std::stringstream ss;
628 0 : ss << "SourcePulsarDistribution: Random position according to pulsar distribution";
629 0 : ss << "rEarth = " << rEarth / kpc << " kpc and ";
630 0 : ss << "zg = " << zg / kpc << " kpc and ";
631 0 : ss << "beta = " << beta << " and ";
632 0 : ss << "r_blur = " << rBlur << " and ";
633 0 : ss << "theta blur = " << thetaBlur << "\n";
634 0 : description = ss.str();
635 0 : }
636 :
637 : // ----------------------------------------------------------------------------
638 1 : SourceUniform1D::SourceUniform1D(double minD, double maxD, bool withCosmology) {
639 1 : this->withCosmology = withCosmology;
640 1 : if (withCosmology) {
641 1 : this->minD = comoving2LightTravelDistance(minD);
642 1 : this->maxD = comoving2LightTravelDistance(maxD);
643 : } else {
644 0 : this->minD = minD;
645 0 : this->maxD = maxD;
646 : }
647 1 : setDescription();
648 1 : }
649 :
650 1 : void SourceUniform1D::prepareParticle(ParticleState& particle) const {
651 1 : Random& random = Random::instance();
652 1 : double d = random.rand() * (maxD - minD) + minD;
653 1 : if (withCosmology)
654 1 : d = lightTravel2ComovingDistance(d);
655 1 : particle.setPosition(Vector3d(d, 0, 0));
656 1 : }
657 :
658 1 : void SourceUniform1D::setDescription() {
659 1 : std::stringstream ss;
660 1 : ss << "SourceUniform1D: Random uniform position in D = ";
661 2 : ss << minD / Mpc << " - " << maxD / Mpc << " Mpc";
662 1 : if (withCosmology)
663 1 : ss << " (including cosmology)";
664 1 : ss << "\n";
665 1 : description = ss.str();
666 1 : }
667 :
668 : // ----------------------------------------------------------------------------
669 2 : SourceDensityGrid::SourceDensityGrid(ref_ptr<Grid1f> grid) :
670 2 : grid(grid) {
671 : float sum = 0;
672 14 : for (int ix = 0; ix < grid->getNx(); ix++) {
673 116 : for (int iy = 0; iy < grid->getNy(); iy++) {
674 1112 : for (int iz = 0; iz < grid->getNz(); iz++) {
675 1008 : sum += grid->get(ix, iy, iz);
676 1008 : grid->get(ix, iy, iz) = sum;
677 : }
678 : }
679 : }
680 2 : setDescription();
681 2 : }
682 :
683 10001 : void SourceDensityGrid::prepareParticle(ParticleState& particle) const {
684 10001 : Random &random = Random::instance();
685 :
686 : // draw random bin
687 10001 : size_t i = random.randBin(grid->getGrid());
688 10001 : Vector3d pos = grid->positionFromIndex(i);
689 :
690 : // draw uniform position within bin
691 10001 : double dx = random.rand() - 0.5;
692 10001 : double dy = random.rand() - 0.5;
693 10001 : double dz = random.rand() - 0.5;
694 : pos += Vector3d(dx, dy, dz) * grid->getSpacing();
695 :
696 10001 : particle.setPosition(pos);
697 10001 : }
698 :
699 2 : void SourceDensityGrid::setDescription() {
700 2 : description = "SourceDensityGrid: 3D source distribution according to density grid\n";
701 2 : }
702 :
703 : // ----------------------------------------------------------------------------
704 2 : SourceDensityGrid1D::SourceDensityGrid1D(ref_ptr<Grid1f> grid) :
705 2 : grid(grid) {
706 2 : if (grid->getNy() != 1)
707 0 : throw std::runtime_error("SourceDensityGrid1D: Ny != 1");
708 2 : if (grid->getNz() != 1)
709 0 : throw std::runtime_error("SourceDensityGrid1D: Nz != 1");
710 :
711 : float sum = 0;
712 22 : for (int ix = 0; ix < grid->getNx(); ix++) {
713 20 : sum += grid->get(ix, 0, 0);
714 20 : grid->get(ix, 0, 0) = sum;
715 : }
716 2 : setDescription();
717 2 : }
718 :
719 101 : void SourceDensityGrid1D::prepareParticle(ParticleState& particle) const {
720 101 : Random &random = Random::instance();
721 :
722 : // draw random bin
723 101 : size_t i = random.randBin(grid->getGrid());
724 101 : Vector3d pos = grid->positionFromIndex(i);
725 :
726 : // draw uniform position within bin
727 101 : double dx = random.rand() - 0.5;
728 101 : pos.x += dx * grid->getSpacing().x;
729 :
730 101 : particle.setPosition(pos);
731 101 : }
732 :
733 2 : void SourceDensityGrid1D::setDescription() {
734 2 : description = "SourceDensityGrid1D: 1D source distribution according to density grid\n";
735 2 : }
736 :
737 : // ----------------------------------------------------------------------------
738 3 : SourceIsotropicEmission::SourceIsotropicEmission() {
739 3 : setDescription();
740 3 : }
741 :
742 1101 : void SourceIsotropicEmission::prepareParticle(ParticleState& particle) const {
743 1101 : Random &random = Random::instance();
744 1101 : particle.setDirection(random.randVector());
745 1101 : }
746 :
747 3 : void SourceIsotropicEmission::setDescription() {
748 3 : description = "SourceIsotropicEmission: Random isotropic direction\n";
749 3 : }
750 :
751 : // ----------------------------------------------------------------------------
752 1 : SourceDirectedEmission::SourceDirectedEmission(Vector3d mu, double kappa): mu(mu), kappa(kappa) {
753 1 : if (kappa <= 0)
754 0 : throw std::runtime_error("The concentration parameter kappa should be larger than 0.");
755 1 : setDescription();
756 1 : }
757 :
758 1000 : void SourceDirectedEmission::prepareCandidate(Candidate &candidate) const {
759 1000 : Random &random = Random::instance();
760 :
761 : Vector3d muvec = mu / mu.getR();
762 1000 : Vector3d v = random.randFisherVector(muvec, kappa);
763 :
764 1000 : v = v.getUnitVector();
765 1000 : candidate.source.setDirection(v);
766 1000 : candidate.created.setDirection(v);
767 1000 : candidate.previous.setDirection(v);
768 1000 : candidate.current.setDirection(v);
769 :
770 : //set the weight of the particle, see eq. 3.1 of PoS(ICRC2019)447
771 1000 : double pdfVonMises = kappa / (2. * M_PI * (1. - exp(-2. * kappa))) * exp(-kappa * (1. - v.dot(mu)));
772 1000 : double weight = 1. / (4. * M_PI * pdfVonMises);
773 1000 : candidate.setWeight(weight);
774 1000 : }
775 :
776 1 : void SourceDirectedEmission::setDescription() {
777 1 : std::stringstream ss;
778 1 : ss << "SourceDirectedEmission: Random directed emission following the von-Mises-Fisher distribution with mean direction ";
779 1 : ss << mu << " and concentration parameter kappa = ";
780 1 : ss << kappa << "\n";
781 1 : description = ss.str();
782 1 : }
783 :
784 : // ----------------------------------------------------------------------------
785 0 : SourceLambertDistributionOnSphere::SourceLambertDistributionOnSphere(const Vector3d ¢er, double radius, bool inward) :
786 0 : center(center), radius(radius) {
787 0 : this->inward = inward;
788 0 : setDescription();
789 0 : }
790 :
791 0 : void SourceLambertDistributionOnSphere::prepareParticle(ParticleState& particle) const {
792 0 : Random &random = Random::instance();
793 0 : Vector3d normalVector = random.randVector();
794 0 : particle.setPosition(center + normalVector * radius);
795 0 : double sign = inward ? -1 : 1; // negative (positive) Lamberts vector for inward (outward) directed emission
796 0 : particle.setDirection(Vector3d(0, 0, 0) + sign * random.randVectorLamberts(normalVector));
797 0 : }
798 :
799 0 : void SourceLambertDistributionOnSphere::setDescription() {
800 0 : std::stringstream ss;
801 0 : ss << "SourceLambertDistributionOnSphere: Random position and direction on a Sphere with center ";
802 0 : ss << center / kpc << " kpc and ";
803 0 : ss << radius / kpc << " kpc radius\n";
804 0 : description = ss.str();
805 0 : }
806 :
807 : // ----------------------------------------------------------------------------
808 0 : SourceDirection::SourceDirection(Vector3d direction) :
809 0 : direction(direction) {
810 0 : setDescription();
811 0 : }
812 :
813 0 : void SourceDirection::prepareParticle(ParticleState& particle) const {
814 0 : particle.setDirection(direction);
815 0 : }
816 :
817 0 : void SourceDirection::setDescription() {
818 0 : std::stringstream ss;
819 0 : ss << "SourceDirection: Emission direction = " << direction << "\n";
820 0 : description = ss.str();
821 0 : }
822 :
823 : // ----------------------------------------------------------------------------
824 0 : SourceEmissionMap::SourceEmissionMap(EmissionMap *emissionMap) : emissionMap(emissionMap) {
825 0 : setDescription();
826 0 : }
827 :
828 0 : void SourceEmissionMap::prepareCandidate(Candidate &candidate) const {
829 0 : if (emissionMap) {
830 0 : bool accept = emissionMap->checkDirection(candidate.source);
831 0 : candidate.setActive(accept);
832 : }
833 0 : }
834 :
835 0 : void SourceEmissionMap::setDescription() {
836 0 : description = "SourceEmissionMap: accept only directions from emission map\n";
837 0 : }
838 :
839 0 : void SourceEmissionMap::setEmissionMap(EmissionMap *emissionMap) {
840 0 : this->emissionMap = emissionMap;
841 0 : }
842 :
843 : // ----------------------------------------------------------------------------
844 1 : SourceEmissionCone::SourceEmissionCone(Vector3d direction, double aperture) :
845 1 : aperture(aperture) {
846 1 : setDirection(direction);
847 1 : setDescription();
848 :
849 1 : }
850 :
851 1 : void SourceEmissionCone::prepareParticle(ParticleState& particle) const {
852 1 : Random &random = Random::instance();
853 1 : particle.setDirection(random.randConeVector(direction, aperture));
854 1 : }
855 :
856 1 : void SourceEmissionCone::setDirection(Vector3d dir) {
857 1 : if (dir.getR() == 0) {
858 0 : throw std::runtime_error("SourceEmissionCone: The direction vector was a null vector.");
859 : } else {
860 1 : direction = dir.getUnitVector();
861 : }
862 1 : }
863 :
864 1 : void SourceEmissionCone::setDescription() {
865 1 : std::stringstream ss;
866 1 : ss << "SourceEmissionCone: Jetted emission in ";
867 1 : ss << "direction = " << direction << " with ";
868 1 : ss << "half-opening angle = " << aperture << " rad\n";
869 1 : description = ss.str();
870 1 : }
871 :
872 : // ----------------------------------------------------------------------------
873 1 : SourceRedshift::SourceRedshift(double z) :
874 1 : z(z) {
875 1 : setDescription();
876 1 : }
877 :
878 1 : void SourceRedshift::prepareCandidate(Candidate& candidate) const {
879 1 : candidate.setRedshift(z);
880 1 : }
881 :
882 1 : void SourceRedshift::setDescription() {
883 1 : std::stringstream ss;
884 1 : ss << "SourceRedshift: Redshift z = " << z << "\n";
885 1 : description = ss.str();
886 1 : }
887 :
888 : // ----------------------------------------------------------------------------
889 0 : SourceUniformRedshift::SourceUniformRedshift(double zmin, double zmax) :
890 0 : zmin(zmin), zmax(zmax) {
891 0 : setDescription();
892 0 : }
893 :
894 0 : void SourceUniformRedshift::prepareCandidate(Candidate& candidate) const {
895 0 : double z = Random::instance().randUniform(zmin, zmax);
896 0 : candidate.setRedshift(z);
897 0 : }
898 :
899 0 : void SourceUniformRedshift::setDescription() {
900 0 : std::stringstream ss;
901 0 : ss << "SourceUniformRedshift: Uniform redshift in z = ";
902 0 : ss << zmin << " - " << zmax << "\n";
903 0 : description = ss.str();
904 0 : }
905 :
906 : // ----------------------------------------------------------------------------
907 2 : SourceRedshiftEvolution::SourceRedshiftEvolution(double m, double zmin, double zmax) : m(m), zmin(zmin), zmax(zmax) {
908 2 : std::stringstream ss;
909 : ss << "SourceRedshiftEvolution: (1+z)^m, m = " << m;
910 2 : ss << ", z = " << zmin << " - " << zmax << "\n";
911 2 : description = ss.str();
912 2 : }
913 :
914 200 : void SourceRedshiftEvolution::prepareCandidate(Candidate& candidate) const {
915 200 : double x = Random::instance().randUniform(0, 1);
916 : double norm, z;
917 :
918 : // special case: m=-1
919 200 : if ((std::abs(m+1)) < std::numeric_limits<double>::epsilon()) {
920 100 : norm = log1p(zmax) - log1p(zmin);
921 100 : z = exp(norm*x) * (1+zmin) - 1;
922 : } else {
923 100 : norm = ( pow(1+zmax, m+1) - pow(1+zmin, m+1) ) / (m+1);
924 100 : z = pow( norm*(m+1)*x + pow(1+zmin, m+1), 1./(m+1)) - 1;
925 : }
926 200 : candidate.setRedshift(z);
927 200 : }
928 :
929 : // ----------------------------------------------------------------------------
930 1 : SourceRedshift1D::SourceRedshift1D() {
931 1 : setDescription();
932 1 : }
933 :
934 1 : void SourceRedshift1D::prepareCandidate(Candidate& candidate) const {
935 1 : double d = candidate.source.getPosition().getR();
936 1 : double z = comovingDistance2Redshift(d);
937 1 : candidate.setRedshift(z);
938 1 : }
939 :
940 1 : void SourceRedshift1D::setDescription() {
941 1 : description = "SourceRedshift1D: Redshift according to source distance\n";
942 1 : }
943 :
944 : // ----------------------------------------------------------------------------
945 : #ifdef CRPROPA_HAVE_MUPARSER
946 1 : SourceGenericComposition::SourceGenericComposition(double Emin, double Emax, std::string expression, size_t bins) :
947 1 : Emin(Emin), Emax(Emax), expression(expression), bins(bins) {
948 :
949 : // precalculate energy bins
950 1 : double logEmin = ::log10(Emin);
951 1 : double logEmax = ::log10(Emax);
952 1 : double logStep = (logEmax - logEmin) / bins;
953 1 : energy.resize(bins + 1);
954 1026 : for (size_t i = 0; i <= bins; i++) {
955 1025 : energy[i] = ::pow(10, logEmin + i * logStep);
956 : }
957 1 : setDescription();
958 1 : }
959 :
960 2 : void SourceGenericComposition::add(int id, double weight) {
961 2 : int A = massNumber(id);
962 2 : int Z = chargeNumber(id);
963 :
964 : Nucleus n;
965 2 : n.id = id;
966 :
967 : // calculate nuclei cdf
968 2 : mu::Parser p;
969 : double E;
970 2 : p.DefineVar("E", &E);
971 2 : p.DefineConst("Emin", Emin);
972 2 : p.DefineConst("Emax", Emax);
973 2 : p.DefineConst("bins", bins);
974 2 : p.DefineConst("A", (double)A);
975 2 : p.DefineConst("Z", (double)Z);
976 :
977 2 : p.DefineConst("MeV", MeV);
978 2 : p.DefineConst("GeV", GeV);
979 2 : p.DefineConst("TeV", TeV);
980 2 : p.DefineConst("PeV", PeV);
981 2 : p.DefineConst("EeV", EeV);
982 :
983 2 : p.SetExpr(expression);
984 :
985 : // calculate pdf
986 2 : n.cdf.resize(bins + 1);
987 :
988 2052 : for (std::size_t i=0; i<=bins; ++i) {
989 2050 : E = energy[i];
990 2050 : n.cdf[i] = p.Eval();
991 : }
992 :
993 : // integrate
994 2050 : for (std::size_t i=bins; i>0; --i) {
995 2048 : n.cdf[i] = (n.cdf[i-1] + n.cdf[i]) * (energy[i] - energy[i-1]) / 2;
996 : }
997 2 : n.cdf[0] = 0;
998 :
999 : // cumulate
1000 2050 : for (std::size_t i=1; i<=bins; ++i) {
1001 2048 : n.cdf[i] += n.cdf[i-1];
1002 : }
1003 :
1004 2 : nuclei.push_back(n);
1005 :
1006 : // update composition cdf
1007 2 : if (cdf.size() == 0)
1008 1 : cdf.push_back(weight * n.cdf.back());
1009 : else
1010 1 : cdf.push_back(cdf.back() + weight * n.cdf.back());
1011 2 : }
1012 :
1013 0 : void SourceGenericComposition::add(int A, int Z, double a) {
1014 0 : add(nucleusId(A, Z), a);
1015 0 : }
1016 :
1017 100000 : void SourceGenericComposition::prepareParticle(ParticleState& particle) const {
1018 100000 : if (nuclei.size() == 0)
1019 0 : throw std::runtime_error("SourceComposition: No source isotope set");
1020 :
1021 100000 : Random &random = Random::instance();
1022 :
1023 :
1024 : // draw random particle type
1025 100000 : size_t iN = random.randBin(cdf);
1026 : const Nucleus &n = nuclei.at(iN);
1027 100000 : particle.setId(n.id);
1028 :
1029 : // random energy
1030 100000 : double E = interpolate(random.rand() * n.cdf.back(), n.cdf, energy);
1031 100000 : particle.setEnergy(E);
1032 100000 : }
1033 :
1034 1 : void SourceGenericComposition::setDescription() {
1035 1 : description = "Generice source composition" + expression;
1036 1 : }
1037 :
1038 : #endif
1039 :
1040 : // ----------------------------------------------------------------------------
1041 :
1042 1 : SourceTag::SourceTag(std::string tag) {
1043 1 : setTag(tag);
1044 1 : }
1045 :
1046 1 : void SourceTag::prepareCandidate(Candidate &cand) const {
1047 1 : cand.setTagOrigin(sourceTag);
1048 1 : }
1049 :
1050 1 : void SourceTag::setDescription() {
1051 1 : description = "SourceTag: " + sourceTag;
1052 1 : }
1053 :
1054 1 : void SourceTag::setTag(std::string tag) {
1055 1 : sourceTag = tag;
1056 1 : setDescription();
1057 1 : }
1058 :
1059 : // ----------------------------------------------------------------------------
1060 :
1061 0 : SourceMassDistribution::SourceMassDistribution(ref_ptr<Density> density, double max, double x, double y, double z) :
1062 0 : density(density), maxDensity(max), xMin(-x), xMax(x), yMin(-y), yMax(y), zMin(-z), zMax(z) {}
1063 :
1064 0 : void SourceMassDistribution::setMaximalDensity(double maxDensity) {
1065 0 : if (maxDensity <= 0) {
1066 0 : KISS_LOG_WARNING << "SourceMassDistribution: maximal density must be larger than 0. Nothing changed.\n";
1067 0 : return;
1068 : }
1069 0 : this->maxDensity = maxDensity;
1070 : }
1071 :
1072 0 : void SourceMassDistribution::setXrange(double xMin, double xMax) {
1073 0 : if (xMin > xMax) {
1074 0 : KISS_LOG_WARNING << "SourceMassDistribution: minimal x-value must not exceed the maximal one\n";
1075 0 : return;
1076 : }
1077 0 : this -> xMin = xMin;
1078 0 : this -> xMax = xMax;
1079 : }
1080 :
1081 0 : void SourceMassDistribution::setYrange(double yMin, double yMax) {
1082 0 : if (yMin > yMax) {
1083 0 : KISS_LOG_WARNING << "SourceMassDistribution: minimal y-value must not exceed the maximal one\n";
1084 0 : return;
1085 : }
1086 0 : this -> yMin = yMin;
1087 0 : this -> yMax = yMax;
1088 : }
1089 :
1090 0 : void SourceMassDistribution::setZrange(double zMin, double zMax) {
1091 0 : if (zMin > zMax) {
1092 0 : KISS_LOG_WARNING << "SourceMassDistribution: minimal z-value must not exceed the maximal one\n";
1093 0 : return;
1094 : }
1095 0 : this -> zMin = zMin;
1096 0 : this -> zMax = zMax;
1097 : }
1098 :
1099 0 : Vector3d SourceMassDistribution::samplePosition() const {
1100 : Vector3d pos;
1101 0 : Random &rand = Random::instance();
1102 :
1103 0 : for (int i = 0; i < maxTries; i++) {
1104 0 : pos.x = rand.randUniform(xMin, xMax);
1105 0 : pos.y = rand.randUniform(yMin, yMax);
1106 0 : pos.z = rand.randUniform(zMin, zMax);
1107 :
1108 0 : double n_density = density->getDensity(pos) / maxDensity;
1109 0 : double n_test = rand.rand();
1110 0 : if (n_test < n_density) {
1111 : return pos;
1112 : }
1113 : }
1114 0 : KISS_LOG_WARNING << "SourceMassDistribution: sampling a position was not possible within "
1115 0 : << maxTries << " tries. Please check the maximum density or increse the number of maximal tries. \n";
1116 : return Vector3d(0.);
1117 : }
1118 :
1119 0 : void SourceMassDistribution::prepareParticle(ParticleState &state) const {
1120 0 : Vector3d pos = samplePosition();
1121 0 : state.setPosition(pos);
1122 0 : }
1123 :
1124 0 : void SourceMassDistribution::setMaximalTries(int tries) {
1125 0 : this -> maxTries = tries;
1126 0 : }
1127 :
1128 0 : std::string SourceMassDistribution::getDescription() {
1129 0 : std::stringstream ss;
1130 0 : ss << "SourceMassDistribuion: following the density distribution :\n";
1131 0 : ss << "\t" << density -> getDescription();
1132 0 : ss << "with a maximal density of " << maxDensity << " / m^3 \n";
1133 0 : ss << "using the sampling range: \n";
1134 0 : ss << "\t x in [" << xMin / kpc << " ; " << xMax / kpc << "] kpc \n";
1135 0 : ss << "\t y in [" << yMin / kpc << " ; " << yMax / kpc << "] kpc \n";
1136 0 : ss << "\t z in [" << zMin / kpc << " ; " << zMax / kpc << "] kpc \n";
1137 0 : ss << "with maximal number of tries for sampling of " << maxTries << "\n";
1138 :
1139 0 : return ss.str();
1140 0 : }
1141 :
1142 : } // namespace crpropa
|
