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