Line data Source code
1 : #include "crpropa/module/DiffusionSDE.h"
2 :
3 :
4 : using namespace crpropa;
5 :
6 : // Defining Cash-Karp coefficients
7 : const double a[] = { 0., 0., 0., 0., 0., 0., 1. / 5., 0., 0., 0., 0.,
8 : 0., 3. / 40., 9. / 40., 0., 0., 0., 0., 3. / 10., -9. / 10., 6. / 5.,
9 : 0., 0., 0., -11. / 54., 5. / 2., -70. / 27., 35. / 27., 0., 0., 1631.
10 : / 55296., 175. / 512., 575. / 13824., 44275. / 110592., 253.
11 : / 4096., 0. };
12 :
13 : const double b[] = { 37. / 378., 0, 250. / 621., 125. / 594., 0., 512.
14 : / 1771. };
15 :
16 : const double bs[] = { 2825. / 27648., 0., 18575. / 48384., 13525.
17 : / 55296., 277. / 14336., 1. / 4. };
18 :
19 :
20 :
21 4 : DiffusionSDE::DiffusionSDE(ref_ptr<MagneticField> magneticField, double tolerance,
22 4 : double minStep, double maxStep, double epsilon) :
23 4 : minStep(0)
24 : {
25 4 : setMagneticField(magneticField);
26 4 : setMaximumStep(maxStep);
27 4 : setMinimumStep(minStep);
28 4 : setTolerance(tolerance);
29 4 : setEpsilon(epsilon);
30 4 : setScale(1.);
31 4 : setAlpha(1./3.);
32 4 : }
33 :
34 3 : DiffusionSDE::DiffusionSDE(ref_ptr<MagneticField> magneticField, ref_ptr<AdvectionField> advectionField, double tolerance, double minStep, double maxStep, double epsilon) :
35 3 : minStep(0)
36 : {
37 6 : setMagneticField(magneticField);
38 3 : setAdvectionField(advectionField);
39 3 : setMaximumStep(maxStep);
40 3 : setMinimumStep(minStep);
41 3 : setTolerance(tolerance);
42 3 : setEpsilon(epsilon);
43 3 : setScale(1.);
44 3 : setAlpha(1./3.);
45 3 : }
46 :
47 480003 : void DiffusionSDE::process(Candidate *candidate) const {
48 :
49 : // save the new previous particle state
50 :
51 480003 : ParticleState ¤t = candidate->current;
52 : candidate->previous = current;
53 :
54 480003 : double h = clip(candidate->getNextStep(), minStep, maxStep) / c_light;
55 480003 : Vector3d PosIn = current.getPosition();
56 480003 : Vector3d DirIn = current.getDirection();
57 :
58 480003 : double time = candidate->getTime();
59 :
60 : // rectilinear propagation for neutral particles
61 : // If an advection field is provided the drift is also included
62 480003 : if (current.getCharge() == 0) {
63 1 : Vector3d dir = current.getDirection();
64 1 : Vector3d Pos = current.getPosition();
65 :
66 : Vector3d LinProp(0.);
67 1 : if (advectionField){
68 0 : driftStep(Pos, LinProp, h, time);
69 : }
70 :
71 1 : current.setPosition(Pos + LinProp + dir*h*c_light);
72 1 : candidate->setCurrentStep(h * c_light);
73 1 : candidate->setNextStep(maxStep);
74 : return;
75 : }
76 :
77 480002 : double z = candidate->getRedshift();
78 480002 : double rig = current.getEnergy() / current.getCharge();
79 :
80 : // Calculate the Diffusion tensor
81 480002 : double BTensor[] = {0., 0., 0., 0., 0., 0., 0., 0., 0.};
82 480002 : calculateBTensor(rig, BTensor, PosIn, DirIn, z);
83 :
84 :
85 : // Generate random numbers
86 480002 : double eta[] = {0., 0., 0.};
87 1920008 : for(size_t i=0; i < 3; i++) {
88 1440006 : eta[i] = Random::instance().randNorm();
89 : }
90 :
91 480002 : double TStep = BTensor[0] * eta[0];
92 480002 : double NStep = BTensor[4] * eta[1];
93 480002 : double BStep = BTensor[8] * eta[2];
94 :
95 : Vector3d TVec(0.);
96 : Vector3d NVec(0.);
97 : Vector3d BVec(0.);
98 :
99 : Vector3d DirOut = Vector3d(0.);
100 :
101 :
102 480002 : double propTime = TStep * sqrt(h) / c_light;
103 : size_t counter = 0;
104 : double r=42.; //arbitrary number larger than one
105 :
106 : do {
107 : Vector3d PosOut = Vector3d(0.);
108 : Vector3d PosErr = Vector3d(0.);
109 480002 : tryStep(PosIn, PosOut, PosErr, z, propTime);
110 : // calculate the relative position error r and the next time step h
111 480002 : r = PosErr.getR() / tolerance;
112 480002 : propTime *= 0.5;
113 480002 : counter += 1;
114 :
115 : // Check for better break condition
116 480002 : } while (r > 1 && fabs(propTime) >= minStep/c_light);
117 :
118 :
119 480002 : size_t stepNumber = pow(2, counter-1);
120 480002 : double allowedTime = TStep * sqrt(h) / c_light / stepNumber;
121 : Vector3d Start = PosIn;
122 : Vector3d PosOut = Vector3d(0.);
123 : Vector3d PosErr = Vector3d(0.);
124 960004 : for (size_t j=0; j<stepNumber; j++) {
125 480002 : tryStep(Start, PosOut, PosErr, z, allowedTime);
126 : Start = PosOut;
127 : }
128 :
129 : // Normalize the tangent vector
130 480002 : TVec = (PosOut-PosIn).getUnitVector();
131 : // Exception: If the magnetic field vanishes: Use only advection.
132 : // If an advection field is not provided --> rectilinear propagation.
133 : double tTest = TVec.getR();
134 480002 : if (tTest != tTest) {
135 2 : Vector3d dir = current.getDirection();
136 2 : Vector3d Pos = current.getPosition();
137 : Vector3d LinProp(0.);
138 2 : if (advectionField){
139 1 : driftStep(Pos, LinProp, h, time);
140 1 : current.setPosition(Pos + LinProp);
141 1 : candidate->setCurrentStep(h*c_light);
142 1 : double newStep = 5*h*c_light;
143 : newStep = clip(newStep, minStep, maxStep);
144 1 : candidate->setNextStep(newStep);
145 : return;
146 : }
147 1 : current.setPosition(Pos + dir*h*c_light);
148 1 : candidate->setCurrentStep(h*c_light);
149 1 : double newStep = 5*h*c_light;
150 1 : newStep = clip(newStep, minStep, maxStep);
151 1 : candidate->setNextStep(newStep);
152 1 : return;
153 : }
154 :
155 : // Choose a random perpendicular vector as the Normal-vector.
156 : // Prevent 'nan's in the NVec-vector in the case of <TVec, NVec> = 0.
157 960000 : while (NVec.getR()==0.){
158 480000 : Vector3d RandomVector = Random::instance().randVector();
159 : NVec = TVec.cross( RandomVector );
160 : }
161 480000 : NVec = NVec.getUnitVector();
162 :
163 : // Calculate the Binormal-vector
164 480000 : BVec = (TVec.cross(NVec)).getUnitVector();
165 :
166 : // Calculate the advection step
167 : Vector3d LinProp(0.);
168 480000 : if (advectionField){
169 120000 : driftStep(PosIn, LinProp, h, time);
170 : }
171 :
172 : // Integration of the SDE with a Mayorama-Euler-method
173 480000 : Vector3d PO = PosOut + LinProp + (NVec * NStep + BVec * BStep) * sqrt(h) ;
174 :
175 : // Throw error message if something went wrong with propagation.
176 : // Deactivate candidate.
177 : bool NaN = std::isnan(PO.getR());
178 480000 : if (NaN == true){
179 0 : candidate->setActive(false);
180 0 : KISS_LOG_WARNING
181 0 : << "\nCandidate with 'nan'-position occured: \n"
182 0 : << "position = " << PO << "\n"
183 0 : << "PosIn = " << PosIn << "\n"
184 0 : << "TVec = " << TVec << "\n"
185 0 : << "TStep = " << std::abs(TStep) << "\n"
186 0 : << "NVec = " << NVec << "\n"
187 : << "NStep = " << NStep << "\n"
188 0 : << "BVec = " << BVec << "\n"
189 : << "BStep = " << BStep << "\n"
190 0 : << "Candidate is deactivated!\n";
191 : return;
192 : }
193 :
194 : //DirOut = (PO - PosIn - LinProp).getUnitVector(); //Advection does not change the momentum vector
195 : // Random direction around the tangential direction accounts for the pitch angle average.
196 480000 : DirOut = Random::instance().randConeVector(TVec, M_PI/2.);
197 480000 : current.setPosition(PO);
198 480000 : current.setDirection(DirOut);
199 480000 : candidate->setCurrentStep(h * c_light);
200 :
201 : double nextStep;
202 480000 : if (stepNumber>1){
203 0 : nextStep = h*pow(stepNumber, -2.)*c_light;
204 : }
205 : else {
206 480000 : nextStep = 4 * h*c_light;
207 : }
208 :
209 480000 : candidate->setNextStep(nextStep);
210 :
211 : // Debugging and Testing
212 : // Delete comments if additional information should be stored in candidate
213 : // This property "arcLength" can be interpreted as the effective arclength
214 : // of the propagation along a magnetic field line.
215 :
216 : /*
217 : const std::string AL = "arcLength";
218 : if (candidate->hasProperty(AL) == false){
219 : double arcLen = (TStep + NStep + BStep) * sqrt(h);
220 : candidate->setProperty(AL, arcLen);
221 : return;
222 : }
223 : else {
224 : double arcLen = candidate->getProperty(AL);
225 : arcLen += (TStep + NStep + BStep) * sqrt(h);
226 : candidate->setProperty(AL, arcLen);
227 : }
228 : */
229 :
230 : }
231 :
232 :
233 960004 : void DiffusionSDE::tryStep(const Vector3d &PosIn, Vector3d &POut, Vector3d &PosErr,double z, double propStep) const {
234 :
235 : Vector3d k[] = {Vector3d(0.),Vector3d(0.),Vector3d(0.),Vector3d(0.),Vector3d(0.),Vector3d(0.)};
236 : POut = PosIn;
237 : //calculate the sum k_i * b_i
238 6720028 : for (size_t i = 0; i < 6; i++) {
239 :
240 : Vector3d y_n = PosIn;
241 20160084 : for (size_t j = 0; j < i; j++)
242 14400060 : y_n += k[j] * a[i * 6 + j] * propStep;
243 :
244 : // update k_i = direction of the regular magnetic mean field
245 5760024 : Vector3d BField = getMagneticFieldAtPosition(y_n, z);
246 :
247 5760024 : k[i] = BField.getUnitVector() * c_light;
248 :
249 5760024 : POut += k[i] * b[i] * propStep;
250 5760024 : PosErr += (k[i] * (b[i] - bs[i])) * propStep / kpc;
251 :
252 : }
253 960004 : }
254 :
255 120001 : void DiffusionSDE::driftStep(const Vector3d &pos, Vector3d &linProp, double h, double t) const {
256 120001 : Vector3d advField = getAdvectionFieldAtPosition(pos, t);
257 : linProp += advField * h;
258 120001 : return;
259 : }
260 :
261 480002 : void DiffusionSDE::calculateBTensor(double r, double BTen[], Vector3d pos, Vector3d dir, double z) const {
262 :
263 480002 : double DifCoeff = scale * 6.1e24 * pow((std::abs(r) / 4.0e9), alpha);
264 480002 : BTen[0] = pow( 2 * DifCoeff, 0.5);
265 480002 : BTen[4] = pow(2 * epsilon * DifCoeff, 0.5);
266 480002 : BTen[8] = pow(2 * epsilon * DifCoeff, 0.5);
267 480002 : return;
268 :
269 : }
270 :
271 :
272 7 : void DiffusionSDE::setMinimumStep(double min) {
273 7 : if (min < 0)
274 0 : throw std::runtime_error("DiffusionSDE: minStep < 0 ");
275 7 : if (min > maxStep)
276 0 : throw std::runtime_error("DiffusionSDE: minStep > maxStep");
277 7 : minStep = min;
278 7 : }
279 :
280 7 : void DiffusionSDE::setMaximumStep(double max) {
281 7 : if (max < minStep)
282 0 : throw std::runtime_error("DiffusionSDE: maxStep < minStep");
283 7 : maxStep = max;
284 7 : }
285 :
286 :
287 7 : void DiffusionSDE::setTolerance(double tol) {
288 7 : if ((tol > 1) or (tol < 0))
289 : throw std::runtime_error(
290 0 : "DiffusionSDE: tolerance error not in range 0-1");
291 7 : tolerance = tol;
292 7 : }
293 :
294 7 : void DiffusionSDE::setEpsilon(double e) {
295 7 : if ((e > 1) or (e < 0))
296 : throw std::runtime_error(
297 0 : "DiffusionSDE: epsilon not in range 0-1");
298 7 : epsilon = e;
299 7 : }
300 :
301 :
302 7 : void DiffusionSDE::setAlpha(double a) {
303 7 : if ((a > 2.) or (a < 0))
304 : throw std::runtime_error(
305 0 : "DiffusionSDE: alpha not in range 0-2");
306 7 : alpha = a;
307 7 : }
308 :
309 7 : void DiffusionSDE::setScale(double s) {
310 7 : if (s < 0)
311 : throw std::runtime_error(
312 0 : "DiffusionSDE: Scale error: Scale < 0");
313 7 : scale = s;
314 7 : }
315 :
316 7 : void DiffusionSDE::setMagneticField(ref_ptr<MagneticField> f) {
317 7 : magneticField = f;
318 7 : }
319 :
320 3 : void DiffusionSDE::setAdvectionField(ref_ptr<AdvectionField> f) {
321 3 : advectionField = f;
322 3 : }
323 :
324 1 : double DiffusionSDE::getMinimumStep() const {
325 1 : return minStep;
326 : }
327 :
328 1 : double DiffusionSDE::getMaximumStep() const {
329 1 : return maxStep;
330 : }
331 :
332 1 : double DiffusionSDE::getTolerance() const {
333 1 : return tolerance;
334 : }
335 :
336 1 : double DiffusionSDE::getEpsilon() const {
337 1 : return epsilon;
338 : }
339 :
340 5 : double DiffusionSDE::getAlpha() const {
341 5 : return alpha;
342 : }
343 :
344 5 : double DiffusionSDE::getScale() const {
345 5 : return scale;
346 : }
347 :
348 0 : ref_ptr<MagneticField> DiffusionSDE::getMagneticField() const {
349 0 : return magneticField;
350 : }
351 :
352 5760025 : Vector3d DiffusionSDE::getMagneticFieldAtPosition(Vector3d pos, double z) const {
353 : Vector3d B(0, 0, 0);
354 : try {
355 : // check if field is valid and use the field vector at the
356 : // position pos with the redshift z
357 5760025 : if (magneticField.valid())
358 5760025 : B = magneticField->getField(pos, z);
359 : }
360 0 : catch (std::exception &e) {
361 0 : KISS_LOG_ERROR << "DiffusionSDE: Exception in DiffusionSDE::getMagneticFieldAtPosition.\n"
362 0 : << e.what();
363 0 : }
364 5760025 : return B;
365 : }
366 :
367 0 : ref_ptr<AdvectionField> DiffusionSDE::getAdvectionField() const {
368 0 : return advectionField;
369 : }
370 :
371 120002 : Vector3d DiffusionSDE::getAdvectionFieldAtPosition(Vector3d pos, double t) const {
372 : Vector3d AdvField(0.);
373 : try {
374 : // check if field is valid and use the field vector at the
375 : // position pos
376 120002 : if (advectionField.valid())
377 120002 : AdvField = advectionField->getField(pos, t);
378 : }
379 0 : catch (std::exception &e) {
380 0 : KISS_LOG_ERROR << "DiffusionSDE: Exception in DiffusionSDE::getAdvectionFieldAtPosition.\n"
381 0 : << e.what();
382 0 : }
383 120002 : return AdvField;
384 : }
385 :
386 0 : std::string DiffusionSDE::getDescription() const {
387 0 : std::stringstream s;
388 0 : s << "minStep: " << minStep / kpc << " kpc, ";
389 0 : s << "maxStep: " << maxStep / kpc << " kpc, ";
390 0 : s << "tolerance: " << tolerance << "\n";
391 :
392 0 : if (epsilon != 0.1) {
393 0 : s << "epsilon: " << epsilon << ", ";
394 : }
395 :
396 0 : if (alpha != 1./3.) {
397 0 : s << "alpha: " << alpha << "\n";
398 : }
399 :
400 0 : if (scale != 1.) {
401 0 : s << "D_0: " << scale*6.1e24 << " m^2/s" << "\n";
402 : }
403 :
404 0 : return s.str();
405 0 : }
|
