Garfield++ v2r0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Medium.cc
Go to the documentation of this file.
1#include <iostream>
2#include <iomanip>
3#include <cmath>
4
5#include "Medium.hh"
8#include "Random.hh"
9#include "Numerics.hh"
10
11namespace {
12
13void PrintNotImplemented(const std::string& cls, const std::string& fcn) {
14
15 std::cerr << cls << "::" << fcn << ": Function is not implemented.\n";
16}
17
18}
19
20namespace Garfield {
21
22int Medium::m_idCounter = -1;
23
25 : m_className("Medium"),
26 m_id(++m_idCounter),
27 m_name(""),
28 m_temperature(293.15),
29 m_pressure(760.),
30 m_epsilon(1.),
31 m_nComponents(1),
32 m_z(1.),
33 m_a(0.),
34 m_density(0.),
35 m_driftable(false),
36 m_microscopic(false),
37 m_ionisable(false),
38 m_w(0.),
39 m_fano(0.),
40 m_isChanged(true),
41 m_debug(false),
42 m_map2d(false) {
43
44 // Initialise the transport tables.
45 m_bFields.assign(1, 0.);
46 m_bAngles.assign(1, 0.);
47
56
57 m_hasHoleVelocityE = false;
58 m_hasHoleVelocityB = false;
60 m_hasHoleDiffLong = false;
61 m_hasHoleDiffTrans = false;
62 m_hasHoleTownsend = false;
63 m_hasHoleAttachment = false;
64 m_hasHoleDiffTens = false;
65
66 m_hasIonMobility = false;
67 m_hasIonDiffLong = false;
68 m_hasIonDiffTrans = false;
70
85
93
97
98 // Set the default grid.
99 SetFieldGrid(100., 100000., 20, true, 0., 0., 1, 0., 0., 1);
100}
101
103
104void Medium::SetTemperature(const double t) {
105
106 if (t <= 0.) {
107 std::cerr << m_className << "::SetTemperature:\n";
108 std::cerr << " Temperature [K] must be greater than zero.\n";
109 return;
110 }
111 m_temperature = t;
112 m_isChanged = true;
113}
114
115void Medium::SetPressure(const double p) {
116
117 if (p <= 0.) {
118 std::cerr << m_className << "::SetPressure:\n";
119 std::cerr << " Pressure [Torr] must be greater than zero.\n";
120 return;
121 }
122 m_pressure = p;
123 m_isChanged = true;
124}
125
126void Medium::SetDielectricConstant(const double eps) {
127
128 if (eps < 1.) {
129 std::cerr << m_className << "::SetDielectricConstant:\n";
130 std::cerr << " Dielectric constant must be >= 1.\n";
131 return;
132 }
133 m_epsilon = eps;
134 m_isChanged = true;
135}
136
138
139 return m_density * AtomicMassUnit * m_a;
140}
141
142void Medium::GetComponent(const unsigned int i,
143 std::string& label, double& f) {
144
145 if (i >= m_nComponents) {
146 std::cerr << m_className << "::GetComponent:\n";
147 std::cerr << " Index out of range.\n";
148 }
149
150 label = m_name;
151 f = 1.;
152}
153
154void Medium::SetAtomicNumber(const double z) {
155
156 if (z < 1.) {
157 std::cerr << m_className << "::SetAtomicNumber:\n";
158 std::cerr << " Atomic number must be >= 1.\n";
159 return;
160 }
161 m_z = z;
162 m_isChanged = true;
163}
164
165void Medium::SetAtomicWeight(const double a) {
166
167 if (a <= 0.) {
168 std::cerr << m_className << "::SetAtomicWeight:\n";
169 std::cerr << " Atomic weight must be greater than zero.\n";
170 return;
171 }
172 m_a = a;
173 m_isChanged = true;
174}
175
176void Medium::SetNumberDensity(const double n) {
177
178 if (n <= 0.) {
179 std::cerr << m_className << "::SetNumberDensity:\n";
180 std::cerr << " Density [cm-3] must be greater than zero.\n";
181 return;
182 }
183 m_density = n;
184 m_isChanged = true;
185}
186
187void Medium::SetMassDensity(const double rho) {
188
189 if (rho <= 0.) {
190 std::cerr << m_className << "::SetMassDensity:\n";
191 std::cerr << " Density [g/cm3] must be greater than zero.\n";
192 return;
193 }
194
195 if (m_a <= 0.) {
196 std::cerr << m_className << "::SetMassDensity:\n";
197 std::cerr << " Atomic weight is not defined.\n";
198 return;
199 }
200 m_density = rho / (AtomicMassUnit * m_a);
201 m_isChanged = true;
202}
203
204bool Medium::ElectronVelocity(const double ex, const double ey, const double ez,
205 const double bx, const double by, const double bz,
206 double& vx, double& vy, double& vz) {
207
208 vx = vy = vz = 0.;
209 // Make sure there is at least a table of velocities along E.
210 if (!m_hasElectronVelocityE) return false;
211
212 // Compute the magnitude of the electric field.
213 const double e = sqrt(ex * ex + ey * ey + ez * ez);
214 const double e0 = ScaleElectricField(e);
215 if (e < Small || e0 < Small) return false;
216
217 // Compute the magnitude of the magnetic field.
218 const double b = sqrt(bx * bx + by * by + bz * bz);
219
220 // Compute the angle between B field and E field.
221 const double ebang = GetAngle(ex, ey, ez, bx, by, bz, e, b);
222
223 if (b < Small) {
224 // No magnetic field.
225
226 // Calculate the velocity along E.
227 double ve = 0.;
228 if (m_map2d) {
230 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, ve,
232 std::cerr << m_className << "::ElectronVelocity:\n";
233 std::cerr << " Interpolation of velocity along E failed.\n";
234 return false;
235 }
236 } else {
239 }
240 const double q = -1.;
241 const double mu = q * ve / e;
242 vx = mu * ex;
243 vy = mu * ey;
244 vz = mu * ez;
245
247 // Magnetic field, velocities along ExB and Bt available
248
249 // Compute unit vectors along E, E x B and Bt.
250 double ue[3] = {ex / e, ey / e, ez / e};
251 double uexb[3] = {ey * bz - ez * by, ez * bx - ex * bz, ex * by - ey * bx};
252 const double exb =
253 sqrt(uexb[0] * uexb[0] + uexb[1] * uexb[1] + uexb[2] * uexb[2]);
254 if (exb > 0.) {
255 uexb[0] /= exb;
256 uexb[1] /= exb;
257 uexb[2] /= exb;
258 } else {
259 uexb[0] = ue[0];
260 uexb[1] = ue[1];
261 uexb[2] = ue[2];
262 }
263
264 double ubt[3] = {uexb[1] * ez - uexb[2] * ey, uexb[2] * ex - uexb[0] * ez,
265 uexb[0] * ey - uexb[1] * ex};
266 const double bt = sqrt(ubt[0] * ubt[0] + ubt[1] * ubt[1] + ubt[2] * ubt[2]);
267
268 if (bt > 0.) {
269 ubt[0] /= bt;
270 ubt[1] /= bt;
271 ubt[2] /= bt;
272 } else {
273 ubt[0] = ue[0];
274 ubt[1] = ue[1];
275 ubt[2] = ue[2];
276 }
277
278 if (m_debug) {
279 std::cout << std::setprecision(5);
280 std::cout << m_className << "::ElectronVelocity:\n";
281 std::cout << " unit vector along E: (" << ue[0] << ", " << ue[1]
282 << ", " << ue[2] << ")\n";
283 std::cout << " unit vector along E x B: (" << uexb[0] << ", "
284 << uexb[1] << ", " << uexb[2] << ")\n";
285 std::cout << " unit vector along Bt: (" << ubt[0] << ", " << ubt[1]
286 << ", " << ubt[2] << ")\n";
287 }
288
289 // Calculate the velocities in all directions.
290 double ve = 0., vbt = 0., vexb = 0.;
291 if (m_map2d) {
293 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, ve,
295 std::cerr << m_className << "::ElectronVelocity:\n";
296 std::cerr << " Interpolation of velocity along E failed.\n";
297 return false;
298 }
300 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, vexb,
302 std::cerr << m_className << "::ElectronVelocity:\n";
303 std::cerr << " Interpolation of velocity along ExB failed.\n";
304 return false;
305 }
307 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, vbt,
309 std::cerr << m_className << "::ElectronVelocity:\n";
310 std::cerr << " Interpolation of velocity along Bt failed.\n";
311 return false;
312 }
313 } else {
323 }
324 const double q = -1.;
325 if (ex * bx + ey * by + ez * bz > 0.) vbt = fabs(vbt);
326 else vbt = -fabs(vbt);
327 vx = q * (ve * ue[0] + q * q * vbt * ubt[0] + q * vexb * uexb[0]);
328 vy = q * (ve * ue[1] + q * q * vbt * ubt[1] + q * vexb * uexb[1]);
329 vz = q * (ve * ue[2] + q * q * vbt * ubt[2] + q * vexb * uexb[2]);
330
331 } else {
332 // Magnetic field, velocities along ExB, Bt not available
333
334 // Calculate the velocity along E.
335 double ve = 0.;
336 if (m_map2d) {
338 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, ve,
340 std::cerr << m_className << "::ElectronVelocity:\n";
341 std::cerr << " Interpolation of velocity along E failed.\n";
342 return false;
343 }
344 } else {
348 }
349
350 const double q = -1.;
351 const double mu = q * ve / e;
352 const double eb = bx * ex + by * ey + bz * ez;
353 const double nom = 1. + pow(mu * b, 2);
354 vx = mu * (ex + mu * (ey * bz - ez * by) + mu * mu * bx * eb) / nom;
355 vy = mu * (ey + mu * (ez * bx - ex * bz) + mu * mu * by * eb) / nom;
356 vz = mu * (ez + mu * (ex * by - ey * bx) + mu * mu * bz * eb) / nom;
357 }
358
359 return true;
360}
361
362bool Medium::ElectronDiffusion(const double ex, const double ey,
363 const double ez, const double bx,
364 const double by, const double bz, double& dl,
365 double& dt) {
366
367 dl = dt = 0.;
368 // Compute the magnitude of the electric field.
369 const double e = sqrt(ex * ex + ey * ey + ez * ez);
370 const double e0 = ScaleElectricField(e);
371 if (e < Small || e0 < Small) return true;
372
373 if (m_map2d) {
374 // Compute the magnitude of the magnetic field.
375 const double b = sqrt(bx * bx + by * by + bz * bz);
376 // Compute the angle between B field and E field.
377 const double ebang = GetAngle(ex, ey, ez, bx, by, bz, e, b);
378
379 // Interpolate.
382 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, dl,
384 dl = 0.;
385 }
386 }
389 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, dt,
391 dt = 0.;
392 }
393 }
394 } else {
399 }
404 }
405 }
406
407 // If no data available, calculate
408 // the diffusion coefficients using the Einstein relation
410 const double d = sqrt(2. * BoltzmannConstant * m_temperature / e);
411 if (!m_hasElectronDiffLong) dl = d;
412 if (!m_hasElectronDiffTrans) dt = d;
413 }
414 // Verify values and apply scaling.
415 if (dl < 0.) dl = 0.;
416 if (dt < 0.) dt = 0.;
417 dl = ScaleDiffusion(dl);
418 dt = ScaleDiffusion(dt);
419
420 return true;
421}
422
423bool Medium::ElectronDiffusion(const double ex, const double ey,
424 const double ez, const double bx,
425 const double by, const double bz,
426 double cov[3][3]) {
427
428 // Initialise the tensor.
429 cov[0][0] = cov[0][1] = cov[0][2] = 0.;
430 cov[1][0] = cov[1][1] = cov[1][2] = 0.;
431 cov[2][0] = cov[2][1] = cov[2][2] = 0.;
432
433 if (!m_hasElectronDiffTens) return false;
434
435 // Compute the magnitude of the electric field.
436 const double e = sqrt(ex * ex + ey * ey + ez * ez);
437 const double e0 = ScaleElectricField(e);
438 if (e < Small || e0 < Small) return true;
439
440 if (m_map2d) {
441 // Compute the magnitude of the magnetic field.
442 const double b = sqrt(bx * bx + by * by + bz * bz);
443 // Compute the angle between B field and E field.
444 const double ebang = GetAngle(ex, ey, ez, bx, by, bz, e, b);
445
446 // Interpolate.
447 double diff = 0.;
448 for (int l = 0; l < 6; ++l) {
450 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, diff,
452 diff = 0.;
453 }
454 // Apply scaling.
455 diff = ScaleDiffusionTensor(diff);
456 if (l < 3) {
457 cov[l][l] = diff;
458 } else if (l == 3) {
459 cov[0][1] = cov[1][0] = diff;
460 } else if (l == 4) {
461 cov[0][2] = cov[2][0] = diff;
462 } else if (l == 5) {
463 cov[1][2] = cov[2][1] = diff;
464 }
465 }
466 } else {
467 // Interpolate.
468 for (int l = 0; l < 6; ++l) {
469 double diff =
473 // Apply scaling.
474 diff = ScaleDiffusionTensor(diff);
475 if (l < 3) {
476 cov[l][l] = diff;
477 } else if (l == 3) {
478 cov[0][1] = cov[1][0] = diff;
479 } else if (l == 4) {
480 cov[0][2] = cov[2][0] = diff;
481 } else if (l == 5) {
482 cov[1][2] = cov[2][1] = diff;
483 }
484 }
485 }
486
487 return true;
488}
489
490bool Medium::ElectronTownsend(const double ex, const double ey, const double ez,
491 const double bx, const double by, const double bz,
492 double& alpha) {
493
494 alpha = 0.;
495 if (tabElectronTownsend.empty()) return false;
496 // Compute the magnitude of the electric field.
497 const double e = sqrt(ex * ex + ey * ey + ez * ez);
498 const double e0 = ScaleElectricField(e);
499 if (e < Small || e0 < Small) return true;
500
501 if (m_map2d) {
502 // Compute the magnitude of the magnetic field.
503 const double b = sqrt(bx * bx + by * by + bz * bz);
504 // Compute the angle between B field and E field.
505 const double ebang = GetAngle(ex, ey, ez, bx, by, bz, e, b);
506
507 // Interpolate.
508 if (e0 < m_eFields[thrElectronTownsend]) {
510 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, alpha,
511 1)) {
512 alpha = -30.;
513 }
514 } else {
516 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, alpha,
518 alpha = -30.;
519 }
520 }
521 } else {
522 // Interpolate.
523 if (e0 < m_eFields[thrElectronTownsend]) {
524 alpha = Interpolate1D(e0, tabElectronTownsend[0][0], m_eFields, 1,
526 } else {
530 }
531 }
532
533 if (alpha < -20.) {
534 alpha = 0.;
535 } else {
536 alpha = exp(alpha);
537 }
538
539 // Apply scaling.
540 alpha = ScaleTownsend(alpha);
541 return true;
542}
543
544bool Medium::ElectronAttachment(const double ex, const double ey,
545 const double ez, const double bx,
546 const double by, const double bz, double& eta) {
547
548 eta = 0.;
549 if (!m_hasElectronAttachment) return false;
550 // Compute the magnitude of the electric field.
551 const double e = sqrt(ex * ex + ey * ey + ez * ez);
552 const double e0 = ScaleElectricField(e);
553 if (e < Small || e0 < Small) return true;
554
555 if (m_map2d) {
556 // Compute the magnitude of the magnetic field.
557 const double b = sqrt(bx * bx + by * by + bz * bz);
558 // Compute the angle between B field and E field.
559 const double ebang = GetAngle(ex, ey, ez, bx, by, bz, e, b);
560
561 // Interpolate.
564 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, eta,
565 1)) {
566 eta = -30.;
567 }
568 } else {
570 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, eta,
572 eta = -30.;
573 }
574 }
575 } else {
576 // Interpolate.
580 } else {
581 eta =
585 }
586 }
587
588 if (eta < -20.) {
589 eta = 0.;
590 } else {
591 eta = exp(eta);
592 }
593
594 // Apply scaling.
595 eta = ScaleAttachment(eta);
596 return true;
597}
598
599bool Medium::ElectronLorentzAngle(const double ex, const double ey,
600 const double ez, const double bx,
601 const double by, const double bz,
602 double& lor) {
603
604 lor = 0.;
605 if (!m_hasElectronLorentzAngle) return false;
606 // Compute the magnitude of the electric field.
607 const double e = sqrt(ex * ex + ey * ey + ez * ez);
608 const double e0 = ScaleElectricField(e);
609 if (e < Small || e0 < Small) return true;
610
611 if (m_map2d) {
612 // Compute the magnitude of the magnetic field.
613 const double b = sqrt(bx * bx + by * by + bz * bz);
614 // Compute the angle between B field and E field.
615 const double ebang = GetAngle(ex, ey, ez, bx, by, bz, e, b);
616
617 // Interpolate.
619 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, lor,
621 lor = 0.;
622 }
623 } else {
624 // Interpolate.
625 lor =
629 }
630 // Apply scaling.
631 lor = ScaleLorentzAngle(lor);
632 return true;
633}
634
635double Medium::GetElectronEnergy(const double px, const double py,
636 const double pz, double& vx, double& vy,
637 double& vz, const int band) {
638
639 if (band > 0) {
640 std::cerr << m_className << "::GetElectronEnergy:\n";
641 std::cerr << " Unknown band index.\n";
642 }
643
644 vx = SpeedOfLight * px / ElectronMass;
645 vy = SpeedOfLight * py / ElectronMass;
646 vz = SpeedOfLight * pz / ElectronMass;
647
648 return 0.5 * (px * px + py * py + pz * pz) / ElectronMass;
649}
650
651void Medium::GetElectronMomentum(const double e, double& px, double& py,
652 double& pz, int& band) {
653
654 const double p = sqrt(2. * ElectronMass * e) / SpeedOfLight;
655 RndmDirection(px, py, pz, p);
656 band = -1;
657}
658
659double Medium::GetElectronNullCollisionRate(const int /*band*/) {
660
661 if (m_debug) PrintNotImplemented(m_className, "GetElectronNullCollisionRate");
662 return 0.;
663}
664
665double Medium::GetElectronCollisionRate(const double /*e*/,
666 const int /*band*/) {
667
668 if (m_debug) PrintNotImplemented(m_className, "GetElectronCollisionRate");
669 return 0.;
670}
671
672bool Medium::GetElectronCollision(const double e, int& type, int& level,
673 double& e1, double& dx, double& dy,
674 double& dz, int& nion, int& ndxc, int& band) {
675
676 type = level = -1;
677 e1 = e;
678 nion = ndxc = band = 0;
679 RndmDirection(dx, dy, dz);
680
681 if (m_debug) PrintNotImplemented(m_className, "GetElectronCollision");
682 return false;
683}
684
685bool Medium::GetIonisationProduct(const unsigned int /*i*/,
686 int& type, double& energy) const {
687
688 if (m_debug) PrintNotImplemented(m_className, "GetIonisationProduct");
689 type = 0;
690 energy = 0.;
691 return false;
692}
693
694bool Medium::GetDeexcitationProduct(const unsigned int /*i*/,
695 double& t, double& s,
696 int& type, double& energy) const {
697
698 if (m_debug) PrintNotImplemented(m_className, "GetDeexcitationProduct");
699 t = s = energy = 0.;
700 type = 0;
701 return false;
702}
703
704bool Medium::HoleVelocity(const double ex, const double ey, const double ez,
705 const double bx, const double by, const double bz,
706 double& vx, double& vy, double& vz) {
707
708 vx = vy = vz = 0.;
709 // Make sure there is at least a table of velocities along E.
710 if (!m_hasHoleVelocityE) return false;
711
712 // Compute the magnitude of the electric field.
713 const double e = sqrt(ex * ex + ey * ey + ez * ez);
714 const double e0 = ScaleElectricField(e);
715 if (e < Small || e0 < Small) return true;
716
717 // Compute the magnitude of the magnetic field.
718 const double b = sqrt(bx * bx + by * by + bz * bz);
719
720 // Compute the angle between B field and E field.
721 const double ebang = GetAngle(ex, ey, ez, bx, by, bz, e, b);
722
723 if (b < Small) {
724 // No magnetic field.
725 // Calculate the velocity along E.
726 double ve = 0.;
727 if (m_map2d) {
729 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, ve,
731 std::cerr << m_className << "::HoleVelocity:\n";
732 std::cerr << " Interpolation of velocity along E failed.\n";
733 return false;
734 }
735 } else {
738 }
739 const double q = 1.;
740 const double mu = q * ve / e;
741 vx = mu * ex;
742 vy = mu * ey;
743 vz = mu * ez;
744
746 // Magnetic field, velocities along ExB and Bt available
747
748 // Compute unit vectors along E, E x B and Bt.
749 double ue[3] = {ex / e, ey / e, ez / e};
750 double uexb[3] = {ey * bz - ez * by, ez * bx - ex * bz, ex * by - ey * bx};
751 const double exb =
752 sqrt(uexb[0] * uexb[0] + uexb[1] * uexb[1] + uexb[2] * uexb[2]);
753 if (exb > 0.) {
754 uexb[0] /= exb;
755 uexb[1] /= exb;
756 uexb[2] /= exb;
757 } else {
758 uexb[0] = ue[0];
759 uexb[1] = ue[1];
760 uexb[2] = ue[2];
761 }
762
763 double ubt[3] = {uexb[1] * ez - uexb[2] * ey, uexb[2] * ex - uexb[0] * ez,
764 uexb[0] * ey - uexb[1] * ex};
765 const double bt = sqrt(ubt[0] * ubt[0] + ubt[1] * ubt[1] + ubt[2] * ubt[2]);
766
767 if (bt > 0.) {
768 ubt[0] /= bt;
769 ubt[1] /= bt;
770 ubt[2] /= bt;
771 } else {
772 ubt[0] = ue[0];
773 ubt[1] = ue[1];
774 ubt[2] = ue[2];
775 }
776
777 // Calculate the velocities in all directions.
778 double ve = 0., vbt = 0., vexb = 0.;
779 if (m_map2d) {
781 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, ve,
783 std::cerr << m_className << "::HoleVelocity:\n";
784 std::cerr << " Interpolation of velocity along E failed.\n";
785 return false;
786 }
788 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, vexb,
790 std::cerr << m_className << "::HoleVelocity:\n";
791 std::cerr << " Interpolation of velocity along ExB failed.\n";
792 return false;
793 }
795 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, vbt,
797 std::cerr << m_className << "::HoleVelocity:\n";
798 std::cerr << " Interpolation of velocity along Bt failed.\n";
799 return false;
800 }
801 } else {
808 }
809 const double q = 1.;
810 if (ex * bx + ey * by + ez * bz > 0.) vbt = fabs(vbt);
811 else vbt = -fabs(vbt);
812 vx = q * (ve * ue[0] + q * q * vbt * ubt[0] + q * vexb * uexb[0]);
813 vy = q * (ve * ue[1] + q * q * vbt * ubt[1] + q * vexb * uexb[1]);
814 vz = q * (ve * ue[2] + q * q * vbt * ubt[2] + q * vexb * uexb[2]);
815
816 } else {
817 // Magnetic field, velocities along ExB, Bt not available
818
819 // Calculate the velocity along E.
820 double ve = 0.;
821 if (m_map2d) {
823 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, ve,
825 std::cerr << m_className << "::HoleVelocity:\n";
826 std::cerr << " Interpolation of velocity along E failed.\n";
827 return false;
828 }
829 } else {
832 }
833
834 const double q = 1.;
835 const double mu = q * ve / e;
836 const double eb = bx * ex + by * ey + bz * ez;
837 const double nom = 1. + pow(mu * b, 2);
838 vx = mu * (ex + mu * (ey * bz - ez * by) + mu * mu * bx * eb) / nom;
839 vy = mu * (ey + mu * (ez * bx - ex * bz) + mu * mu * by * eb) / nom;
840 vz = mu * (ez + mu * (ex * by - ey * bx) + mu * mu * bz * eb) / nom;
841 }
842
843 return true;
844}
845
846bool Medium::HoleDiffusion(const double ex, const double ey, const double ez,
847 const double bx, const double by, const double bz,
848 double& dl, double& dt) {
849
850 dl = dt = 0.;
851 // Compute the magnitude of the electric field.
852 const double e = sqrt(ex * ex + ey * ey + ez * ez);
853 const double e0 = ScaleElectricField(e);
854 if (e < Small || e0 < Small) return true;
855
856 if (m_map2d) {
857 // Compute the magnitude of the magnetic field.
858 const double b = sqrt(bx * bx + by * by + bz * bz);
859 // Compute the angle between B field and E field.
860 const double ebang = GetAngle(ex, ey, ez, bx, by, bz, e, b);
861
862 // Interpolate.
863 if (m_hasHoleDiffLong) {
865 m_bFields.size(), m_eFields.size(), ebang, b, e0, dl,
867 dl = 0.;
868 }
869 }
870 if (m_hasHoleDiffTrans) {
872 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, dt,
874 dt = 0.;
875 }
876 }
877 } else {
878 if (m_hasHoleDiffLong) {
881 }
882 if (m_hasHoleDiffTrans) {
885 }
886 }
887
888 // If no data available, calculate
889 // the diffusion coefficients using the Einstein relation
890 if (!m_hasHoleDiffLong) {
891 dl = sqrt(2. * BoltzmannConstant * m_temperature / e);
892 }
893 if (!m_hasHoleDiffTrans) {
894 dt = sqrt(2. * BoltzmannConstant * m_temperature / e);
895 }
896
897 // Verify values and apply scaling.
898 if (dl < 0.) dl = 0.;
899 if (dt < 0.) dt = 0.;
900 dl = ScaleDiffusion(dl);
901 dt = ScaleDiffusion(dt);
902
903 return true;
904}
905
906bool Medium::HoleDiffusion(const double ex, const double ey, const double ez,
907 const double bx, const double by, const double bz,
908 double cov[3][3]) {
909
910 // Initialise the tensor.
911 cov[0][0] = cov[0][1] = cov[0][2] = 0.;
912 cov[1][0] = cov[1][1] = cov[1][2] = 0.;
913 cov[2][0] = cov[2][1] = cov[2][2] = 0.;
914
915 if (!m_hasHoleDiffTens) return false;
916
917 // Compute the magnitude of the electric field.
918 const double e = sqrt(ex * ex + ey * ey + ez * ez);
919 const double e0 = ScaleElectricField(e);
920 if (e < Small || e0 < Small) return true;
921
922 if (m_map2d) {
923 // Compute the magnitude of the magnetic field.
924 const double b = sqrt(bx * bx + by * by + bz * bz);
925 // Compute the angle between B field and E field.
926 const double ebang = GetAngle(ex, ey, ez, bx, by, bz, e, b);
927
928 // Interpolate.
929 double diff = 0.;
930 for (int l = 0; l < 6; ++l) {
932 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, diff,
934 diff = 0.;
935 }
936 // Apply scaling.
937 diff = ScaleDiffusionTensor(diff);
938 if (l < 3) {
939 cov[l][l] = diff;
940 } else if (l == 3) {
941 cov[0][1] = cov[1][0] = diff;
942 } else if (l == 4) {
943 cov[0][2] = cov[2][0] = diff;
944 } else if (l == 5) {
945 cov[1][2] = cov[2][1] = diff;
946 }
947 }
948 } else {
949 // Interpolate.
950 for (int l = 0; l < 6; ++l) {
951 double diff =
954 // Apply scaling.
955 diff = ScaleDiffusionTensor(diff);
956 if (l < 3) {
957 cov[l][l] = diff;
958 } else if (l == 3) {
959 cov[0][1] = cov[1][0] = diff;
960 } else if (l == 4) {
961 cov[0][2] = cov[2][0] = diff;
962 } else if (l == 5) {
963 cov[1][2] = cov[2][1] = diff;
964 }
965 }
966 }
967
968 return true;
969}
970
971bool Medium::HoleTownsend(const double ex, const double ey, const double ez,
972 const double bx, const double by, const double bz,
973 double& alpha) {
974
975 alpha = 0.;
976 if (!m_hasHoleTownsend) return false;
977 // Compute the magnitude of the electric field.
978 const double e = sqrt(ex * ex + ey * ey + ez * ez);
979 const double e0 = ScaleElectricField(e);
980 if (e < Small || e0 < Small) return true;
981
982 if (m_map2d) {
983 // Compute the magnitude of the magnetic field.
984 const double b = sqrt(bx * bx + by * by + bz * bz);
985 // Compute the angle between B field and E field.
986 const double ebang = GetAngle(ex, ey, ez, bx, by, bz, e, b);
987
988 // Interpolate.
989 if (e0 < m_eFields[thrHoleTownsend]) {
991 m_bFields.size(), m_eFields.size(), ebang, b, e0, alpha, 1)) {
992 alpha = -30.;
993 }
994 } else {
996 m_bFields.size(), m_eFields.size(), ebang, b, e0, alpha,
998 alpha = -30.;
999 }
1000 }
1001 } else {
1002 // Interpolate.
1003 if (e0 < m_eFields[thrHoleTownsend]) {
1004 alpha = Interpolate1D(e0, tabHoleTownsend[0][0], m_eFields, 1,
1006 } else {
1009 }
1010 }
1011
1012 if (alpha < -20.) {
1013 alpha = 0.;
1014 } else {
1015 alpha = exp(alpha);
1016 }
1017
1018 // Apply scaling.
1019 alpha = ScaleTownsend(alpha);
1020 return true;
1021}
1022
1023bool Medium::HoleAttachment(const double ex, const double ey, const double ez,
1024 const double bx, const double by, const double bz,
1025 double& eta) {
1026
1027 eta = 0.;
1028 if (!m_hasHoleAttachment) return false;
1029 // Compute the magnitude of the electric field.
1030 const double e = sqrt(ex * ex + ey * ey + ez * ez);
1031 const double e0 = ScaleElectricField(e);
1032 if (e < Small || e0 < Small) return true;
1033
1034 if (m_map2d) {
1035 // Compute the magnitude of the magnetic field.
1036 const double b = sqrt(bx * bx + by * by + bz * bz);
1037 // Compute the angle between B field and E field.
1038 const double ebang = GetAngle(ex, ey, ez, bx, by, bz, e, b);
1039
1040 // Interpolate.
1041 if (e0 < m_eFields[thrHoleAttachment]) {
1043 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, eta,
1044 1)) {
1045 eta = -30.;
1046 }
1047 } else {
1049 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, eta,
1051 eta = -30.;
1052 }
1053 }
1054 } else {
1055 // Interpolate.
1056 if (e0 < m_eFields[thrHoleAttachment]) {
1057 eta = Interpolate1D(e0, tabHoleAttachment[0][0], m_eFields, 1,
1059 } else {
1060 eta = Interpolate1D(e0, tabHoleAttachment[0][0], m_eFields,
1063 }
1064 }
1065
1066 if (eta < -20.) {
1067 eta = 0.;
1068 } else {
1069 eta = exp(eta);
1070 }
1071
1072 // Apply scaling.
1073 eta = ScaleAttachment(eta);
1074 return true;
1075}
1076
1077bool Medium::IonVelocity(const double ex, const double ey, const double ez,
1078 const double bx, const double by, const double bz,
1079 double& vx, double& vy, double& vz) {
1080
1081 vx = vy = vz = 0.;
1082 if (!m_hasIonMobility) return false;
1083 // Compute the magnitude of the electric field.
1084 const double e = sqrt(ex * ex + ey * ey + ez * ez);
1085 const double e0 = ScaleElectricField(e);
1086 if (e < Small || e0 < Small) return true;
1087 // Compute the magnitude of the electric field.
1088 const double b = sqrt(bx * bx + by * by + bz * bz);
1089
1090 double mu = 0.;
1091 if (m_map2d) {
1092 // Compute the angle between B field and E field.
1093 const double ebang = GetAngle(ex, ey, ez, bx, by, bz, e, b);
1094
1096 m_bFields.size(), m_eFields.size(), ebang, b, e0, mu,
1097 m_intpMobility)) {
1098 mu = 0.;
1099 }
1100 } else {
1103 }
1104
1105 const double q = 1.;
1106 mu *= q;
1107 if (b < Small) {
1108 vx = mu * ex;
1109 vy = mu * ey;
1110 vz = mu * ez;
1111 } else {
1112 const double eb = bx * ex + by * ey + bz * ez;
1113 const double nom = 1. + pow(mu * b, 2);
1114 vx = mu * (ex + mu * (ey * bz - ez * by) + mu * mu * bx * eb) / nom;
1115 vy = mu * (ey + mu * (ez * bx - ex * bz) + mu * mu * by * eb) / nom;
1116 vz = mu * (ez + mu * (ex * by - ey * bx) + mu * mu * bz * eb) / nom;
1117 }
1118
1119 return true;
1120}
1121
1122bool Medium::IonDiffusion(const double ex, const double ey, const double ez,
1123 const double bx, const double by, const double bz,
1124 double& dl, double& dt) {
1125
1126 dl = dt = 0.;
1127 // Compute the magnitude of the electric field.
1128 const double e = sqrt(ex * ex + ey * ey + ez * ez);
1129 const double e0 = ScaleElectricField(e);
1130 if (e < Small || e0 < Small) return true;
1131
1132 if (m_map2d) {
1133 // Compute the magnitude of the magnetic field.
1134 const double b = sqrt(bx * bx + by * by + bz * bz);
1135 // Compute the angle between B field and E field.
1136 const double ebang = GetAngle(ex, ey, ez, bx, by, bz, e, b);
1137
1138 // Interpolate.
1139 if (m_hasIonDiffLong) {
1141 m_bFields.size(), m_eFields.size(), ebang, b, e0, dl,
1142 m_intpDiffusion)) {
1143 dl = 0.;
1144 }
1145 }
1146 if (m_hasIonDiffTrans) {
1148 m_bFields.size(), m_eFields.size(), ebang, b, e0, dt,
1149 m_intpDiffusion)) {
1150 dt = 0.;
1151 }
1152 }
1153 } else {
1154 if (m_hasIonDiffLong) {
1157 }
1158 if (m_hasIonDiffTrans) {
1161 }
1162 }
1163
1164 // If no data available, calculate
1165 // the diffusion coefficients using the Einstein relation
1166 if (!m_hasIonDiffLong) {
1167 dl = sqrt(2. * BoltzmannConstant * m_temperature / e);
1168 }
1169 if (!m_hasIonDiffTrans) {
1170 dt = sqrt(2. * BoltzmannConstant * m_temperature / e);
1171 }
1172
1173 return true;
1174}
1175
1176bool Medium::IonDissociation(const double ex, const double ey, const double ez,
1177 const double bx, const double by, const double bz,
1178 double& diss) {
1179
1180 diss = 0.;
1181 if (!m_hasIonDissociation) return false;
1182 // Compute the magnitude of the electric field.
1183 const double e = sqrt(ex * ex + ey * ey + ez * ez);
1184 const double e0 = ScaleElectricField(e);
1185 if (e < Small || e0 < Small) return true;
1186
1187 if (m_map2d) {
1188 // Compute the magnitude of the magnetic field.
1189 const double b = sqrt(bx * bx + by * by + bz * bz);
1190 // Compute the angle between B field and E field.
1191 const double ebang = GetAngle(ex, ey, ez, bx, by, bz, e, b);
1192 // Interpolate.
1193 if (e0 < m_eFields[thrIonDissociation]) {
1195 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, diss,
1196 1)) {
1197 diss = -30.;
1198 }
1199 } else {
1201 m_bAngles.size(), m_bFields.size(), m_eFields.size(), ebang, b, e0, diss,
1203 diss = -30.;
1204 }
1205 }
1206 } else {
1207 // Interpolate.
1208 if (e0 < m_eFields[thrIonDissociation]) {
1209 diss = Interpolate1D(e0, tabIonDissociation[0][0], m_eFields, 1,
1211 } else {
1212 diss = Interpolate1D(e0, tabHoleTownsend[0][0], m_eFields,
1215 }
1216 }
1217
1218 if (diss < -20.) {
1219 diss = 0.;
1220 } else {
1221 diss = exp(diss);
1222 }
1223
1224 // Apply scaling.
1225 diss = ScaleDissociation(diss);
1226 return true;
1227}
1228
1229bool Medium::GetOpticalDataRange(double& emin, double& emax,
1230 const unsigned int i) {
1231
1232 if (i >= m_nComponents) {
1233 std::cerr << m_className << "::GetOpticalDataRange:\n";
1234 std::cerr << " Component " << i << " does not exist.\n";
1235 return false;
1236 }
1237
1238 if (m_debug) PrintNotImplemented(m_className, "GetOpticalDataRange");
1239 emin = emax = 0.;
1240 return false;
1241}
1242
1243bool Medium::GetDielectricFunction(const double e, double& eps1, double& eps2,
1244 const unsigned int i) {
1245
1246 if (i >= m_nComponents) {
1247 std::cerr << m_className << "::GetDielectricFunction:\n";
1248 std::cerr << " Component " << i << " does not exist.\n";
1249 return false;
1250 }
1251
1252 if (e < 0.) {
1253 std::cerr << m_className << "::GetDielectricFunction:\n";
1254 std::cerr << " Energy must be > 0.\n";
1255 return false;
1256 }
1257
1258 if (m_debug) PrintNotImplemented(m_className, "GetDielectricFunction");
1259 eps1 = 1.;
1260 eps2 = 0.;
1261 return false;
1262}
1263
1264bool Medium::GetPhotoAbsorptionCrossSection(const double e, double& sigma,
1265 const unsigned int i) {
1266
1267 if (i >= m_nComponents) {
1268 std::cerr << m_className << "::GetPhotoAbsorptionCrossSection:\n";
1269 std::cerr << " Component " << i << " does not exist.\n";
1270 return false;
1271 }
1272
1273 if (e < 0.) {
1274 std::cerr << m_className << "::GetPhotoAbsorptionCrossSection:\n";
1275 std::cerr << " Energy must be > 0.\n";
1276 return false;
1277 }
1278
1279 if (m_debug) {
1280 PrintNotImplemented(m_className, "GetPhotoAbsorptionCrossSection");
1281 }
1282 sigma = 0.;
1283 return false;
1284}
1285
1286double Medium::GetPhotonCollisionRate(const double e) {
1287
1288 double sigma = 0.;
1289 if (!GetPhotoAbsorptionCrossSection(e, sigma)) return 0.;
1290
1291 return sigma * m_density * SpeedOfLight;
1292}
1293
1294bool Medium::GetPhotonCollision(const double e, int& type, int& level,
1295 double& e1, double& ctheta, int& nsec,
1296 double& esec) {
1297
1298 type = level = -1;
1299 e1 = e;
1300 ctheta = 1.;
1301 nsec = 0;
1302 esec = 0.;
1303 return false;
1304}
1305
1307
1308 tabElectronVelocityE.clear();
1309 tabElectronVelocityB.clear();
1310 tabElectronVelocityExB.clear();
1311 m_hasElectronVelocityE = false;
1312 m_hasElectronVelocityB = false;
1314}
1315
1317
1318 tabElectronDiffLong.clear();
1319 tabElectronDiffTrans.clear();
1320 tabElectronDiffTens.clear();
1321 m_hasElectronDiffLong = false;
1322 m_hasElectronDiffTrans = false;
1323 m_hasElectronDiffTens = false;
1324}
1325
1327
1328 tabElectronTownsend.clear();
1329}
1330
1332
1333 tabElectronAttachment.clear();
1335}
1336
1338
1341}
1342
1344
1345 tabHoleVelocityE.clear();
1346 tabHoleVelocityB.clear();
1347 tabHoleVelocityExB.clear();
1348 m_hasHoleVelocityE = false;
1349 m_hasHoleVelocityB = false;
1350 m_hasHoleVelocityExB = false;
1351}
1352
1354
1355 tabHoleDiffLong.clear();
1356 tabHoleDiffTrans.clear();
1357 tabHoleDiffTens.clear();
1358 m_hasHoleDiffLong = false;
1359 m_hasHoleDiffTrans = false;
1360 m_hasHoleDiffTens = false;
1361}
1362
1364
1365 tabHoleTownsend.clear();
1366 m_hasHoleTownsend = false;
1367}
1368
1370
1371 tabHoleAttachment.clear();
1372 m_hasHoleAttachment = false;
1373}
1374
1376
1377 tabIonMobility.clear();
1378 m_hasIonMobility = false;
1379}
1380
1382
1383 tabIonDiffLong.clear();
1384 tabIonDiffTrans.clear();
1385 m_hasIonDiffLong = false;
1386 m_hasIonDiffTrans = false;
1387}
1388
1390
1391 tabIonDissociation.clear();
1392 m_hasIonDissociation = false;
1393}
1394
1395void Medium::SetFieldGrid(double emin, double emax, int ne, bool logE,
1396 double bmin, double bmax, int nb, double amin,
1397 double amax, int na) {
1398
1399 // Check if the requested E-field range makes sense.
1400 if (ne <= 0) {
1401 std::cerr << m_className << "::SetFieldGrid:\n";
1402 std::cerr << " Number of E-fields must be > 0.\n";
1403 return;
1404 }
1405
1406 if (emin < 0. || emax < 0.) {
1407 std::cerr << m_className << "::SetFieldGrid:\n";
1408 std::cerr << " Electric fields must be positive.\n";
1409 return;
1410 }
1411
1412 if (emax < emin) {
1413 std::cerr << m_className << "::SetFieldGrid:\n";
1414 std::cerr << " Swapping min./max. E-field.\n";
1415 const double etemp = emin;
1416 emin = emax;
1417 emax = etemp;
1418 }
1419
1420 double estep = 0.;
1421 if (logE) {
1422 // Logarithmic scale
1423 if (emin < Small) {
1424 std::cerr << m_className << "::SetFieldGrid:\n";
1425 std::cerr << " Min. E-field must be non-zero for log. scale.\n";
1426 return;
1427 }
1428 if (ne == 1) {
1429 std::cerr << m_className << "::SetFieldGrid:\n";
1430 std::cerr << " Number of E-fields must be > 1 for log. scale.\n";
1431 return;
1432 }
1433 estep = pow(emax / emin, 1. / (ne - 1.));
1434 } else {
1435 // Linear scale
1436 if (ne > 1) estep = (emax - emin) / (ne - 1.);
1437 }
1438
1439 // Check if the requested B-field range makes sense.
1440 if (nb <= 0) {
1441 std::cerr << m_className << "::SetFieldGrid:\n";
1442 std::cerr << " Number of B-fields must be > 0.\n";
1443 return;
1444 }
1445 if (bmax < 0. || bmin < 0.) {
1446 std::cerr << m_className << "::SetFieldGrid:\n";
1447 std::cerr << " Magnetic fields must be positive.\n";
1448 return;
1449 }
1450 if (bmax < bmin) {
1451 std::cerr << m_className << "::SetFieldGrid:\n";
1452 std::cerr << " Swapping min./max. B-field.\n";
1453 const double btemp = bmin;
1454 bmin = bmax;
1455 bmax = btemp;
1456 }
1457
1458 double bstep = 0.;
1459 if (nb > 1) bstep = (bmax - bmin) / (nb - 1.);
1460
1461 // Check if the requested angular range makes sense.
1462 if (na <= 0) {
1463 std::cerr << m_className << "::SetFieldGrid:\n";
1464 std::cerr << " Number of angles must be > 0.\n";
1465 return;
1466 }
1467 if (amax < 0. || amin < 0.) {
1468 std::cerr << m_className << "::SetFieldGrid:\n";
1469 std::cerr << " Angles must be positive.\n";
1470 return;
1471 }
1472 if (amax < amin) {
1473 std::cerr << m_className << "::SetFieldGrid:\n";
1474 std::cerr << " Swapping min./max. angle.\n";
1475 const double atemp = amin;
1476 amin = amax;
1477 amax = atemp;
1478 }
1479 double astep = 0.;
1480 if (na > 1) astep = (amax - amin) / (na - 1.);
1481
1482 // Setup the field grids.
1483 std::vector<double> eFieldsNew(ne);
1484 std::vector<double> bFieldsNew(nb);
1485 std::vector<double> bAnglesNew(na);
1486 for (int i = 0; i < ne; ++i) {
1487 if (logE) {
1488 eFieldsNew[i] = emin * pow(estep, i);
1489 } else {
1490 eFieldsNew[i] = emin + i * estep;
1491 }
1492 }
1493 for (int i = 0; i < nb; ++i) {
1494 bFieldsNew[i] = bmin + i * bstep;
1495 }
1496 if (na == 1 && nb == 1 && fabs(bmin) < 1.e-4) {
1497 bAnglesNew[0] = HalfPi;
1498 } else {
1499 for (int i = 0; i < na; ++i) {
1500 bAnglesNew[i] = amin + i * astep;
1501 }
1502 }
1503 SetFieldGrid(eFieldsNew, bFieldsNew, bAnglesNew);
1504}
1505
1506void Medium::SetFieldGrid(const std::vector<double>& efields,
1507 const std::vector<double>& bfields,
1508 const std::vector<double>& angles) {
1509
1510 if (efields.empty()) {
1511 std::cerr << m_className << "::SetFieldGrid:\n";
1512 std::cerr << " Number of E-fields must be > 0.\n";
1513 return;
1514 }
1515 if (bfields.empty()) {
1516 std::cerr << m_className << "::SetFieldGrid:\n";
1517 std::cerr << " Number of B-fields must be > 0.\n";
1518 return;
1519 }
1520 if (angles.empty()) {
1521 std::cerr << m_className << "::SetFieldGrid:\n";
1522 std::cerr << " Number of angles must be > 0.\n";
1523 return;
1524 }
1525
1526 // Make sure the values are not negative.
1527 if (efields[0] < 0.) {
1528 std::cerr << m_className << "::SetFieldGrid:\n";
1529 std::cerr << " E-fields must be >= 0.\n";
1530 }
1531 if (bfields[0] < 0.) {
1532 std::cerr << m_className << "::SetFieldGrid:\n";
1533 std::cerr << " B-fields must be >= 0.\n";
1534 }
1535 if (angles[0] < 0.) {
1536 std::cerr << m_className << "::SetFieldGrid:\n";
1537 std::cerr << " Angles must be >= 0.\n";
1538 }
1539
1540 const unsigned int nEfieldsNew = efields.size();
1541 const unsigned int nBfieldsNew = bfields.size();
1542 const unsigned int nAnglesNew = angles.size();
1543 // Make sure the values are in strictly monotonic, ascending order.
1544 if (nEfieldsNew > 1) {
1545 for (unsigned int i = 1; i < nEfieldsNew; ++i) {
1546 if (efields[i] <= efields[i - 1]) {
1547 std::cerr << m_className << "::SetFieldGrid:\n";
1548 std::cerr << " E-fields are not in ascending order.\n";
1549 return;
1550 }
1551 }
1552 }
1553 if (nBfieldsNew > 1) {
1554 for (unsigned int i = 1; i < nBfieldsNew; ++i) {
1555 if (bfields[i] <= bfields[i - 1]) {
1556 std::cerr << m_className << "::SetFieldGrid:\n";
1557 std::cerr << " B-fields are not in ascending order.\n";
1558 return;
1559 }
1560 }
1561 }
1562 if (nAnglesNew > 1) {
1563 for (unsigned int i = 1; i < nAnglesNew; ++i) {
1564 if (angles[i] <= angles[i - 1]) {
1565 std::cerr << m_className << "::SetFieldGrid:\n";
1566 std::cerr << " Angles are not in ascending order.\n";
1567 return;
1568 }
1569 }
1570 }
1571
1572 if (m_debug) {
1573 std::cout << m_className << "::SetFieldGrid:\n";
1574 std::cout << " E-fields:\n";
1575 for (unsigned int i = 0; i < nEfieldsNew; ++i) {
1576 std::cout << " " << efields[i] << "\n";
1577 }
1578 std::cout << " B-fields:\n";
1579 for (unsigned int i = 0; i < nBfieldsNew; ++i) {
1580 std::cout << " " << bfields[i] << "\n";
1581 }
1582 std::cout << " Angles:\n";
1583 for (unsigned int i = 0; i < nAnglesNew; ++i) {
1584 std::cout << " " << angles[i] << "\n";
1585 }
1586 }
1587
1588 // Clone the existing tables.
1589 // Electrons
1590 CloneTable(tabElectronVelocityE, efields, bfields, angles, m_intpVelocity,
1592 "electron velocity along E");
1593 CloneTable(tabElectronVelocityB, efields, bfields, angles, m_intpVelocity,
1595 "electron velocity along Bt");
1596 CloneTable(tabElectronVelocityExB, efields, bfields, angles, m_intpVelocity,
1598 "electron velocity along ExB");
1599 CloneTable(tabElectronDiffLong, efields, bfields, angles, m_intpDiffusion,
1601 "electron longitudinal diffusion");
1602 CloneTable(tabElectronDiffTrans, efields, bfields, angles, m_intpDiffusion,
1604 "electron transverse diffusion");
1605 CloneTable(tabElectronTownsend, efields, bfields, angles, m_intpTownsend,
1607 "electron Townsend coefficient");
1608 CloneTable(tabElectronAttachment, efields, bfields, angles, m_intpAttachment,
1610 "electron attachment coefficient");
1611 CloneTable(tabElectronLorentzAngle, efields, bfields, angles, m_intpLorentzAngle,
1613 "electron attachment coefficient");
1615 CloneTensor(tabElectronDiffTens, 6, efields, bfields, angles, m_intpDiffusion,
1617 "electron diffusion tensor");
1618 }
1619
1620 // Holes
1621 CloneTable(tabHoleVelocityE, efields, bfields, angles, m_intpVelocity,
1623 "hole velocity along E");
1624 CloneTable(tabHoleVelocityB, efields, bfields, angles, m_intpVelocity,
1626 "hole velocity along Bt");
1627 CloneTable(tabHoleVelocityExB, efields, bfields, angles, m_intpVelocity,
1629 "hole velocity along ExB");
1630 CloneTable(tabHoleDiffLong, efields, bfields, angles, m_intpDiffusion,
1632 "hole longitudinal diffusion");
1633 CloneTable(tabHoleDiffTrans, efields, bfields, angles, m_intpDiffusion,
1635 "hole transverse diffusion");
1636 CloneTable(tabHoleTownsend, efields, bfields, angles, m_intpTownsend,
1638 "hole Townsend coefficient");
1639 CloneTable(tabHoleAttachment, efields, bfields, angles, m_intpAttachment,
1641 "hole attachment coefficient");
1642 if (m_hasHoleDiffTens) {
1643 CloneTensor(tabHoleDiffTens, 6, efields, bfields, angles, m_intpDiffusion,
1645 "hole diffusion tensor");
1646 }
1647
1648 // Ions
1649 CloneTable(tabIonMobility, efields, bfields, angles, m_intpMobility,
1651 "ion mobility");
1652 CloneTable(tabIonDiffLong, efields, bfields, angles, m_intpDiffusion,
1654 "ion longitudinal diffusion");
1655 CloneTable(tabIonDiffTrans, efields, bfields, angles, m_intpDiffusion,
1657 "ion transverse diffusion");
1658 CloneTable(tabIonDissociation, efields, bfields, angles, m_intpDissociation,
1660 "ion dissociation");
1661
1662 if (nBfieldsNew > 1 || nAnglesNew > 1) m_map2d = true;
1663 m_eFields.resize(nEfieldsNew);
1664 m_bFields.resize(nBfieldsNew);
1665 m_bAngles.resize(nAnglesNew);
1666 m_eFields = efields;
1667 m_bFields = bfields;
1668 m_bAngles = angles;
1669}
1670
1671void Medium::GetFieldGrid(std::vector<double>& efields,
1672 std::vector<double>& bfields,
1673 std::vector<double>& angles) {
1674
1675 efields = m_eFields;
1676 bfields = m_bFields;
1677 angles = m_bAngles;
1678}
1679
1680bool Medium::GetElectronVelocityE(const unsigned int ie,
1681 const unsigned int ib,
1682 const unsigned int ia, double& v) {
1683
1684 if (ie >= m_eFields.size() || ib >= m_bFields.size() || ia >= m_bAngles.size()) {
1685 std::cerr << m_className << "::GetElectronVelocityE:\n";
1686 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
1687 << ") out of range.\n";
1688 v = 0.;
1689 return false;
1690 }
1692 if (m_debug) {
1693 std::cerr << m_className << "::GetElectronVelocityE:\n";
1694 std::cerr << " Data not available.\n";
1695 }
1696 v = 0.;
1697 return false;
1698 }
1699
1700 v = tabElectronVelocityE[ia][ib][ie];
1701 return true;
1702}
1703
1704bool Medium::GetElectronVelocityExB(const unsigned int ie,
1705 const unsigned int ib,
1706 const unsigned int ia, double& v) {
1707
1708 if (ie >= m_eFields.size() || ib >= m_bFields.size() || ia >= m_bAngles.size()) {
1709 std::cerr << m_className << "::GetElectronVelocityExB:\n";
1710 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
1711 << ") out of range.\n";
1712 v = 0.;
1713 return false;
1714 }
1716 if (m_debug) {
1717 std::cerr << m_className << "::GetElectronVelocityExB:\n";
1718 std::cerr << " Data not available.\n";
1719 }
1720 v = 0.;
1721 return false;
1722 }
1723
1724 v = tabElectronVelocityExB[ia][ib][ie];
1725 return true;
1726}
1727
1728bool Medium::GetElectronVelocityB(const unsigned int ie,
1729 const unsigned int ib,
1730 const unsigned int ia, double& v) {
1731
1732 if (ie >= m_eFields.size() || ib >= m_bFields.size() || ia >= m_bAngles.size()) {
1733 std::cerr << m_className << "::GetElectronVelocityB:\n";
1734 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
1735 << ") out of range.\n";
1736 v = 0.;
1737 return false;
1738 }
1740 if (m_debug) {
1741 std::cerr << m_className << "::GetElectronVelocityB:\n";
1742 std::cerr << " Data not available.\n";
1743 }
1744 v = 0.;
1745 return false;
1746 }
1747
1748 v = tabElectronVelocityB[ia][ib][ie];
1749 return true;
1750}
1751
1753 const unsigned int ib,
1754 const unsigned int ia,
1755 double& dl) {
1756
1757 if (ie >= m_eFields.size() || ib >= m_bFields.size() || ia >= m_bAngles.size()) {
1758 std::cerr << m_className << "::GetElectronLongitudinalDiffusion:\n";
1759 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
1760 << ") out of range.\n";
1761 dl = 0.;
1762 return false;
1763 }
1764 if (!m_hasElectronDiffLong) {
1765 if (m_debug) {
1766 std::cerr << m_className << "::GetElectronLongitudinalDiffusion:\n";
1767 std::cerr << " Data not available.\n";
1768 }
1769 dl = 0.;
1770 return false;
1771 }
1772
1773 dl = tabElectronDiffLong[ia][ib][ie];
1774 return true;
1775}
1776
1777bool Medium::GetElectronTransverseDiffusion(const unsigned int ie,
1778 const unsigned int ib,
1779 const unsigned int ia,
1780 double& dt) {
1781
1782 if (ie >= m_eFields.size() || ib >= m_bFields.size() || ia >= m_bAngles.size()) {
1783 std::cerr << m_className << "::GetElectronTransverseDiffusion:\n";
1784 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
1785 << ") out of range.\n";
1786 dt = 0.;
1787 return false;
1788 }
1790 if (m_debug) {
1791 std::cerr << m_className << "::GetElectronTransverseDiffusion:\n";
1792 std::cerr << " Data not available.\n";
1793 }
1794 dt = 0.;
1795 return false;
1796 }
1797
1798 dt = tabElectronDiffTrans[ia][ib][ie];
1799 return true;
1800}
1801
1802bool Medium::GetElectronTownsend(const unsigned int ie,
1803 const unsigned int ib,
1804 const unsigned int ia, double& alpha) {
1805
1806 if (ie >= m_eFields.size() || ib >= m_bFields.size() || ia >= m_bAngles.size()) {
1807 std::cerr << m_className << "::GetElectronTownsend:\n";
1808 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
1809 << ") out of range.\n";
1810 alpha = 0.;
1811 return false;
1812 }
1813 if (tabElectronTownsend.empty()) {
1814 if (m_debug) {
1815 std::cerr << m_className << "::GetElectronTownsend:\n";
1816 std::cerr << " Data not available.\n";
1817 }
1818 alpha = 0.;
1819 return false;
1820 }
1821
1822 alpha = tabElectronTownsend[ia][ib][ie];
1823 return true;
1824}
1825
1826bool Medium::GetElectronAttachment(const unsigned int ie,
1827 const unsigned int ib,
1828 const unsigned int ia, double& eta) {
1829
1830 if (ie >= m_eFields.size() || ib >= m_bFields.size() || ia >= m_bAngles.size()) {
1831 std::cerr << m_className << "::GetElectronAttachment:\n";
1832 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
1833 << ") out of range.\n";
1834 eta = 0.;
1835 return false;
1836 }
1838 if (m_debug) {
1839 std::cerr << m_className << "::GetElectronAttachment:\n";
1840 std::cerr << " Data not available.\n";
1841 }
1842 eta = 0.;
1843 return false;
1844 }
1845
1846 eta = tabElectronAttachment[ia][ib][ie];
1847 return true;
1848}
1849
1850bool Medium::GetElectronLorentzAngle(const unsigned int ie,
1851 const unsigned int ib,
1852 const unsigned int ia, double& lor) {
1853
1854 if (ie >= m_eFields.size() || ib >= m_bFields.size() || ia >= m_bAngles.size()) {
1855 std::cerr << m_className << "::GetElectronLorentzAngle:\n";
1856 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
1857 << ") out of range.\n";
1858 lor = 0.;
1859 return false;
1860 }
1862 if (m_debug) {
1863 std::cerr << m_className << "::GetElectronLorentzAngle:\n";
1864 std::cerr << " Data not available.\n";
1865 }
1866 lor = 0.;
1867 return false;
1868 }
1869
1870 lor = tabElectronLorentzAngle[ia][ib][ie];
1871 return true;
1872}
1873
1874bool Medium::GetHoleVelocityE(const unsigned int ie,
1875 const unsigned int ib,
1876 const unsigned int ia, double& v) {
1877
1878 if (ie >= m_eFields.size() || ib >= m_bFields.size() || ia >= m_bAngles.size()) {
1879 std::cerr << m_className << "::GetHoleVelocityE:\n";
1880 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
1881 << ") out of range.\n";
1882 v = 0.;
1883 return false;
1884 }
1885 if (!m_hasHoleVelocityE) {
1886 if (m_debug) {
1887 std::cerr << m_className << "::GetHoleVelocityE:\n";
1888 std::cerr << " Data not available.\n";
1889 }
1890 v = 0.;
1891 return false;
1892 }
1893
1894 v = tabHoleVelocityE[ia][ib][ie];
1895 return true;
1896}
1897
1898bool Medium::GetHoleVelocityExB(const unsigned int ie,
1899 const unsigned int ib,
1900 const unsigned int ia, double& v) {
1901
1902 if (ie >= m_eFields.size() || ib >= m_bFields.size() || ia >= m_bAngles.size()) {
1903 std::cerr << m_className << "::GetHoleVelocityExB:\n";
1904 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
1905 << ") out of range.\n";
1906 v = 0.;
1907 return false;
1908 }
1909 if (!m_hasHoleVelocityExB) {
1910 if (m_debug) {
1911 std::cerr << m_className << "::GetHoleVelocityExB:\n";
1912 std::cerr << " Data not available.\n";
1913 }
1914 v = 0.;
1915 return false;
1916 }
1917
1918 v = tabHoleVelocityExB[ia][ib][ie];
1919 return true;
1920}
1921
1922bool Medium::GetHoleVelocityB(const unsigned int ie,
1923 const unsigned int ib,
1924 const unsigned int ia, double& v) {
1925
1926 if (ie >= m_eFields.size() || ib >= m_bFields.size() || ia >= m_bAngles.size()) {
1927 std::cerr << m_className << "::GetHoleVelocityB:\n";
1928 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
1929 << ") out of range.\n";
1930 v = 0.;
1931 return false;
1932 }
1933 if (!m_hasHoleVelocityB) {
1934 if (m_debug) {
1935 std::cerr << m_className << "::GetHoleVelocityB:\n";
1936 std::cerr << " Data not available.\n";
1937 }
1938 v = 0.;
1939 return false;
1940 }
1941
1942 v = tabHoleVelocityB[ia][ib][ie];
1943 return true;
1944}
1945
1946bool Medium::GetHoleLongitudinalDiffusion(const unsigned int ie,
1947 const unsigned int ib,
1948 const unsigned int ia, double& dl) {
1949
1950 if (ie >= m_eFields.size() || ib >= m_bFields.size() || ia >= m_bAngles.size()) {
1951 std::cerr << m_className << "::GetHoleLongitudinalDiffusion:\n";
1952 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
1953 << ") out of range.\n";
1954 dl = 0.;
1955 return false;
1956 }
1957 if (!m_hasHoleDiffLong) {
1958 if (m_debug) {
1959 std::cerr << m_className << "::GetHoleLongitudinalDiffusion:\n";
1960 std::cerr << " Data not available.\n";
1961 }
1962 dl = 0.;
1963 return false;
1964 }
1965
1966 dl = tabHoleDiffLong[ia][ib][ie];
1967 return true;
1968}
1969
1970bool Medium::GetHoleTransverseDiffusion(const unsigned int ie,
1971 const unsigned int ib,
1972 const unsigned int ia, double& dt) {
1973
1974 if (ie >= m_eFields.size() || ib >= m_bFields.size() || ia >= m_bAngles.size()) {
1975 std::cerr << m_className << "::GetHoleTransverseDiffusion:\n";
1976 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
1977 << ") out of range.\n";
1978 dt = 0.;
1979 return false;
1980 }
1981 if (!m_hasHoleDiffTrans) {
1982 if (m_debug) {
1983 std::cerr << m_className << "::GetHoleTransverseDiffusion:\n";
1984 std::cerr << " Data not available.\n";
1985 }
1986 dt = 0.;
1987 return false;
1988 }
1989
1990 dt = tabHoleDiffTrans[ia][ib][ie];
1991 return true;
1992}
1993
1994bool Medium::GetHoleTownsend(const unsigned int ie,
1995 const unsigned int ib,
1996 const unsigned int ia, double& alpha) {
1997
1998 if (ie >= m_eFields.size() || ib >= m_bFields.size() || ia >= m_bAngles.size()) {
1999 std::cerr << m_className << "::GetHoleTownsend:\n";
2000 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
2001 << ") out of range.\n";
2002 alpha = 0.;
2003 return false;
2004 }
2005 if (!m_hasHoleTownsend) {
2006 if (m_debug) {
2007 std::cerr << m_className << "::GetHoleTownsend:\n";
2008 std::cerr << " Data not available.\n";
2009 }
2010 alpha = 0.;
2011 return false;
2012 }
2013
2014 alpha = tabHoleTownsend[ia][ib][ie];
2015 return true;
2016}
2017
2018bool Medium::GetHoleAttachment(const unsigned int ie,
2019 const unsigned int ib,
2020 const unsigned int ia, double& eta) {
2021
2022 if (ie >= m_eFields.size() || ib >= m_bFields.size() || ia >= m_bAngles.size()) {
2023 std::cerr << m_className << "::GetHoleAttachment:\n";
2024 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
2025 << ") out of range.\n";
2026 eta = 0.;
2027 return false;
2028 }
2029 if (!m_hasHoleAttachment) {
2030 if (m_debug) {
2031 std::cerr << m_className << "::GetHoleAttachment:\n";
2032 std::cerr << " Data not available.\n";
2033 }
2034 eta = 0.;
2035 return false;
2036 }
2037
2038 eta = tabHoleAttachment[ia][ib][ie];
2039 return true;
2040}
2041
2042bool Medium::GetIonMobility(const unsigned int ie, const unsigned int ib,
2043 const unsigned int ia, double& mu) {
2044
2045 if (ie >= m_eFields.size() || ib >= m_bFields.size() || ia >= m_bAngles.size()) {
2046 std::cerr << m_className << "::GetIonMobility:\n";
2047 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
2048 << ") out of range.\n";
2049 mu = 0.;
2050 return false;
2051 }
2052 if (!m_hasIonMobility) {
2053 if (m_debug) {
2054 std::cerr << m_className << "::GetIonMobility:\n";
2055 std::cerr << " Data not available.\n";
2056 }
2057 mu = 0.;
2058 return false;
2059 }
2060
2061 mu = tabIonMobility[ia][ib][ie];
2062 return true;
2063}
2064
2065bool Medium::GetIonLongitudinalDiffusion(const unsigned int ie,
2066 const unsigned int ib,
2067 const unsigned int ia, double& dl) {
2068
2069 if (ie >= m_eFields.size() || ib >= m_bFields.size() || ia >= m_bAngles.size()) {
2070 std::cerr << m_className << "::GetIonLongitudinalDiffusion:\n";
2071 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
2072 << ") out of range.\n";
2073 dl = 0.;
2074 return false;
2075 }
2076 if (!m_hasIonDiffLong) {
2077 if (m_debug) {
2078 std::cerr << m_className << "::GetIonLongitudinalDiffusion:\n";
2079 std::cerr << " Data not available.\n";
2080 }
2081 dl = 0.;
2082 return false;
2083 }
2084
2085 dl = tabIonDiffLong[ia][ib][ie];
2086 return true;
2087}
2088
2089bool Medium::GetIonTransverseDiffusion(const unsigned int ie,
2090 const unsigned int ib,
2091 const unsigned int ia, double& dt) {
2092
2093 if (ie >= m_eFields.size() || ib >= m_bFields.size() || ia >= m_bAngles.size()) {
2094 std::cerr << m_className << "::GetIonTransverseDiffusion:\n";
2095 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
2096 << ") out of range.\n";
2097 dt = 0.;
2098 return false;
2099 }
2100 if (!m_hasIonDiffTrans) {
2101 if (m_debug) {
2102 std::cerr << m_className << "::GetIonTransverseDiffusion:\n";
2103 std::cerr << " Data not available.\n";
2104 }
2105 dt = 0.;
2106 return false;
2107 }
2108
2109 dt = tabIonDiffTrans[ia][ib][ie];
2110 return true;
2111}
2112
2113bool Medium::GetIonDissociation(const unsigned int ie,
2114 const unsigned int ib,
2115 const unsigned int ia, double& diss) {
2116
2117 if (ie >= m_eFields.size() || ib >= m_bFields.size() || ia >= m_bAngles.size()) {
2118 std::cerr << m_className << "::GetIonDissociation:\n";
2119 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
2120 << ") out of range.\n";
2121 diss = 0.;
2122 return false;
2123 }
2124 if (!m_hasIonDissociation) {
2125 if (m_debug) {
2126 std::cerr << m_className << "::GetIonDissociation:\n";
2127 std::cerr << " Data not available.\n";
2128 }
2129 diss = 0.;
2130 return false;
2131 }
2132
2133 diss = tabIonDissociation[ia][ib][ie];
2134 return true;
2135}
2136
2137void Medium::CloneTable(std::vector<std::vector<std::vector<double> > >& tab,
2138 const std::vector<double>& efields,
2139 const std::vector<double>& bfields,
2140 const std::vector<double>& angles,
2141 const unsigned int intp,
2142 const unsigned int extrLow,
2143 const unsigned int extrHigh,
2144 const double init, const std::string& label) {
2145
2146 if (m_debug) {
2147 std::cout << m_className << "::CloneTable:\n";
2148 std::cout << " Copying values of " << label << " to new grid.\n";
2149 }
2150
2151 if (tab.empty()) {
2152 if (m_debug) std::cout << " Table is empty.\n";
2153 return;
2154 }
2155 // Get the dimensions of the new grid.
2156 const int nEfieldsNew = efields.size();
2157 const int nBfieldsNew = bfields.size();
2158 const int nAnglesNew = angles.size();
2159
2160 // Create a temporary table to store the values at the new grid points.
2161 std::vector<std::vector<std::vector<double> > > tabClone;
2162 InitParamArrays(nEfieldsNew, nBfieldsNew, nAnglesNew, tabClone, init);
2163
2164 // Fill the temporary table.
2165 for (int i = 0; i < nEfieldsNew; ++i) {
2166 for (int j = 0; j < nBfieldsNew; ++j) {
2167 for (int k = 0; k < nAnglesNew; ++k) {
2168 double val = 0.;
2169 if (m_map2d) {
2171 m_bFields.size(), m_eFields.size(), angles[k], bfields[j],
2172 efields[i], val, intp)) {
2173 std::cerr << m_className << "::SetFieldGrid:\n";
2174 std::cerr << " Interpolation of " << label << " failed.\n";
2175 std::cerr << " Cannot copy value to new grid at: \n";
2176 std::cerr << " E = " << efields[i] << "\n";
2177 std::cerr << " B = " << bfields[j] << "\n";
2178 std::cerr << " angle: " << angles[k] << "\n";
2179 } else {
2180 tabClone[k][j][i] = val;
2181 }
2182 } else {
2183 val = Interpolate1D(efields[i], tab[0][0], m_eFields, intp, extrLow,
2184 extrHigh);
2185 tabClone[k][j][i] = val;
2186 }
2187 }
2188 }
2189 }
2190 // Re-dimension the original table.
2191 InitParamArrays(nEfieldsNew, nBfieldsNew, nAnglesNew, tab, init);
2192 // Copy the values to the original table.
2193 for (int i = 0; i < nEfieldsNew; ++i) {
2194 for (int j = 0; j < nBfieldsNew; ++j) {
2195 for (int k = 0; k < nAnglesNew; ++k) {
2196 tab[k][j][i] = tabClone[k][j][i];
2197 }
2198 }
2199 }
2200 tabClone.clear();
2201}
2202
2204 std::vector<std::vector<std::vector<std::vector<double> > > >& tab,
2205 const unsigned int n,
2206 const std::vector<double>& efields, const std::vector<double>& bfields,
2207 const std::vector<double>& angles,
2208 const unsigned int intp,
2209 const unsigned int extrLow, const unsigned int extrHigh,
2210 const double init,
2211 const std::string& label) {
2212
2213 // If the table does not exist, do nothing.
2214 if (tab.empty()) return;
2215
2216 // Get the dimensions of the new grid.
2217 const unsigned int nEfieldsNew = efields.size();
2218 const unsigned int nBfieldsNew = bfields.size();
2219 const unsigned int nAnglesNew = angles.size();
2220
2221 // Create a temporary table to store the values at the new grid points.
2222 std::vector<std::vector<std::vector<std::vector<double> > > > tabClone;
2223 tabClone.clear();
2224 InitParamTensor(nEfieldsNew, nBfieldsNew, nAnglesNew, n, tabClone, init);
2225
2226 // Fill the temporary table.
2227 for (unsigned int l = 0; l < n; ++l) {
2228 for (unsigned int i = 0; i < nEfieldsNew; ++i) {
2229 for (unsigned int j = 0; j < nBfieldsNew; ++j) {
2230 for (unsigned int k = 0; k < nAnglesNew; ++k) {
2231 double val = 0.;
2232 if (m_map2d) {
2234 m_bFields.size(), m_eFields.size(), angles[k], bfields[j],
2235 efields[i], val, intp)) {
2236 std::cerr << m_className << "::SetFieldGrid:\n";
2237 std::cerr << " Interpolation of " << label << " failed.\n";
2238 std::cerr << " Cannot copy value to new grid at: \n";
2239 std::cerr << " Index: " << l << "\n";
2240 std::cerr << " E = " << efields[i] << "\n";
2241 std::cerr << " B = " << bfields[j] << "\n";
2242 std::cerr << " angle: " << angles[k] << "\n";
2243 } else {
2244 tabClone[l][k][j][i] = val;
2245 }
2246 } else {
2247 val = Interpolate1D(efields[i], tab[l][0][0], m_eFields, intp,
2248 extrLow, extrHigh);
2249 tabClone[l][k][j][i] = val;
2250 }
2251 }
2252 }
2253 }
2254 }
2255 // Re-dimension the original table.
2256 InitParamTensor(nEfieldsNew, nBfieldsNew, nAnglesNew, n, tab, 0.);
2257 // Copy the values to the original table.
2258 for (unsigned int l = 0; l < n; ++l) {
2259 for (unsigned int i = 0; i < nEfieldsNew; ++i) {
2260 for (unsigned int j = 0; j < nBfieldsNew; ++j) {
2261 for (unsigned int k = 0; k < nAnglesNew; ++k) {
2262 tab[l][k][j][i] = tabClone[l][k][j][i];
2263 }
2264 }
2265 }
2266 }
2267}
2268
2269bool Medium::SetIonMobility(const unsigned int ie, const unsigned int ib,
2270 const unsigned int ia, const double mu) {
2271
2272 // Check the index.
2273 if (ie >= m_eFields.size() || ib >= m_bFields.size() ||
2274 ia >= m_bAngles.size()) {
2275 std::cerr << m_className << "::SetIonMobility:\n";
2276 std::cerr << " Index (" << ie << ", " << ib << ", " << ia
2277 << ") out of range.\n";
2278 return false;
2279 }
2280
2281 if (!m_hasIonMobility) {
2282 std::cerr << m_className << "::SetIonMobility:\n";
2283 std::cerr << " Ion mobility table not initialised.\n";
2284 return false;
2285 }
2286
2287 if (mu == 0.) {
2288 std::cerr << m_className << "::SetIonMobility:\n";
2289 std::cerr << " Zero value not permitted.\n";
2290 return false;
2291 }
2292
2293 tabIonMobility[ia][ib][ie] = mu;
2294 if (m_debug) {
2295 std::cout << m_className << "::SetIonMobility:\n";
2296 std::cout << " Ion mobility at E = " << m_eFields[ie]
2297 << " V/cm, B = " << m_bFields[ib] << " T, angle " << m_bAngles[ia]
2298 << " set to " << mu << " cm2/(V ns).\n";
2299 }
2300 return true;
2301}
2302
2303bool Medium::SetIonMobility(const std::vector<double>& efields,
2304 const std::vector<double>& mobilities) {
2305
2306 const int ne = efields.size();
2307 const int nm = mobilities.size();
2308 if (ne != nm) {
2309 std::cerr << m_className << "::SetIonMobility:\n";
2310 std::cerr << " E-field and mobility arrays have different sizes.\n";
2311 return false;
2312 }
2313
2315 const unsigned int nEfields = m_eFields.size();
2316 const unsigned int nBfields = m_bFields.size();
2317 const unsigned int nAngles = m_bAngles.size();
2318 InitParamArrays(nEfields, nBfields, nAngles, tabIonMobility, 0.);
2319 for (unsigned int i = 0; i < nEfields; ++i) {
2320 const double e = m_eFields[i];
2321 const double mu = Interpolate1D(e, mobilities, efields, m_intpMobility,
2323 tabIonMobility[0][0][i] = mu;
2324 }
2325
2326 if (m_map2d) {
2327 for (unsigned int i = 0; i < nAngles; ++i) {
2328 for (unsigned int j = 0; j< nBfields; ++j) {
2329 for (unsigned int k = 0; k < nEfields; ++k) {
2330 tabIonMobility[i][j][k] = tabIonMobility[0][0][k];
2331 }
2332 }
2333 }
2334 }
2335 m_hasIonMobility = true;
2336
2337 return true;
2338}
2339
2340void Medium::SetExtrapolationMethodVelocity(const std::string& extrLow,
2341 const std::string& extrHigh) {
2342
2343 unsigned int iExtr = 0;
2344 if (GetExtrapolationIndex(extrLow, iExtr)) {
2345 m_extrLowVelocity = iExtr;
2346 } else {
2347 std::cerr << m_className << "::SetExtrapolationMethodVelocity:\n";
2348 std::cerr << " Unknown extrapolation method (" << extrLow << ")\n";
2349 }
2350
2351 if (GetExtrapolationIndex(extrHigh, iExtr)) {
2352 m_extrHighVelocity = iExtr;
2353 } else {
2354 std::cerr << m_className << "::SetExtrapolationMethodVelocity:\n";
2355 std::cerr << " Unknown extrapolation method (" << extrHigh << ")\n";
2356 }
2357}
2358
2359void Medium::SetExtrapolationMethodDiffusion(const std::string& extrLow,
2360 const std::string& extrHigh) {
2361
2362 unsigned int iExtr = 0;
2363 if (GetExtrapolationIndex(extrLow, iExtr)) {
2364 m_extrLowDiffusion = iExtr;
2365 } else {
2366 std::cerr << m_className << "::SetExtrapolationMethodDiffusion:\n";
2367 std::cerr << " Unknown extrapolation method (" << extrLow << ")\n";
2368 }
2369
2370 if (GetExtrapolationIndex(extrHigh, iExtr)) {
2371 m_extrHighDiffusion = iExtr;
2372 } else {
2373 std::cerr << m_className << "::SetExtrapolationMethodDiffusion:\n";
2374 std::cerr << " Unknown extrapolation method (" << extrHigh << ")\n";
2375 }
2376}
2377
2378void Medium::SetExtrapolationMethodTownsend(const std::string& extrLow,
2379 const std::string& extrHigh) {
2380
2381 unsigned int iExtr = 0;
2382 if (GetExtrapolationIndex(extrLow, iExtr)) {
2383 m_extrLowTownsend = iExtr;
2384 } else {
2385 std::cerr << m_className << "::SetExtrapolationMethodTownsend:\n";
2386 std::cerr << " Unknown extrapolation method (" << extrLow << ")\n";
2387 }
2388
2389 if (GetExtrapolationIndex(extrHigh, iExtr)) {
2390 m_extrHighTownsend = iExtr;
2391 } else {
2392 std::cerr << m_className << "::SetExtrapolationMethodTownsend:\n";
2393 std::cerr << " Unknown extrapolation method (" << extrHigh << ")\n";
2394 }
2395}
2396
2397void Medium::SetExtrapolationMethodAttachment(const std::string& extrLow,
2398 const std::string& extrHigh) {
2399
2400 unsigned int iExtr = 0;
2401 if (GetExtrapolationIndex(extrLow, iExtr)) {
2402 m_extrLowAttachment = iExtr;
2403 } else {
2404 std::cerr << m_className << "::SetExtrapolationMethodAttachment:\n";
2405 std::cerr << " Unknown extrapolation method (" << extrLow << ")\n";
2406 }
2407
2408 if (GetExtrapolationIndex(extrHigh, iExtr)) {
2409 m_extrHighAttachment = iExtr;
2410 } else {
2411 std::cerr << m_className << "::SetExtrapolationMethodAttachment:\n";
2412 std::cerr << " Unknown extrapolation method (" << extrHigh << ")\n";
2413 }
2414}
2415
2416void Medium::SetExtrapolationMethodIonMobility(const std::string& extrLow,
2417 const std::string& extrHigh) {
2418
2419 unsigned int iExtr = 0;
2420 if (GetExtrapolationIndex(extrLow, iExtr)) {
2421 m_extrLowMobility = iExtr;
2422 } else {
2423 std::cerr << m_className << "::SetExtrapolationMethodIonMobility:\n";
2424 std::cerr << " Unknown extrapolation method (" << extrLow << ")\n";
2425 }
2426 if (GetExtrapolationIndex(extrHigh, iExtr)) {
2427 m_extrHighMobility = iExtr;
2428 } else {
2429 std::cerr << m_className << "::SetExtrapolationMethodIonMobility:\n";
2430 std::cerr << " Unknown extrapolation method (" << extrHigh << ")\n";
2431 }
2432}
2433
2435 const std::string& extrHigh) {
2436
2437 unsigned int iExtr = 0;
2438 if (GetExtrapolationIndex(extrLow, iExtr)) {
2439 m_extrLowDissociation = iExtr;
2440 } else {
2441 std::cerr << m_className << "::SetExtrapolationMethodIonDissociation:\n";
2442 std::cerr << " Unknown extrapolation method (" << extrLow << ")\n";
2443 }
2444
2445 if (GetExtrapolationIndex(extrHigh, iExtr)) {
2446 m_extrHighDissociation = iExtr;
2447 } else {
2448 std::cerr << m_className << "::SetExtrapolationMethodIonDissociation:\n";
2449 std::cerr << " Unknown extrapolation method (" << extrHigh << ")\n";
2450 }
2451}
2452
2453bool Medium::GetExtrapolationIndex(std::string extrStr, unsigned int& extrNb) {
2454
2455 // Convert to upper-case
2456 for (unsigned int i = 0; i < extrStr.length(); ++i) {
2457 extrStr[i] = toupper(extrStr[i]);
2458 }
2459
2460 if (extrStr == "CONST" || extrStr == "CONSTANT") {
2461 extrNb = 0;
2462 } else if (extrStr == "LIN" || extrStr == "LINEAR") {
2463 extrNb = 1;
2464 } else if (extrStr == "EXP" || extrStr == "EXPONENTIAL") {
2465 extrNb = 2;
2466 } else {
2467 return false;
2468 }
2469
2470 return true;
2471}
2472
2473void Medium::SetInterpolationMethodVelocity(const unsigned int intrp) {
2474
2475 if (intrp > 0) m_intpVelocity = intrp;
2476}
2477
2478void Medium::SetInterpolationMethodDiffusion(const unsigned int intrp) {
2479
2480 if (intrp > 0) m_intpDiffusion = intrp;
2481}
2482
2483void Medium::SetInterpolationMethodTownsend(const unsigned int intrp) {
2484
2485 if (intrp > 0) m_intpTownsend = intrp;
2486}
2487
2488void Medium::SetInterpolationMethodAttachment(const unsigned int intrp) {
2489
2490 if (intrp > 0) m_intpAttachment = intrp;
2491}
2492
2493void Medium::SetInterpolationMethodIonMobility(const unsigned int intrp) {
2494
2495 if (intrp > 0) m_intpMobility = intrp;
2496}
2497
2499
2500 if (intrp > 0) m_intpDissociation = intrp;
2501}
2502
2503double Medium::GetAngle(const double ex, const double ey, const double ez,
2504 const double bx, const double by, const double bz,
2505 const double e, const double b) const {
2506
2507 // Ion
2508 /*
2509 if (e * b > 0.) {
2510 const double eb = fabs(ex * bx + ey * by + ez * bz);
2511 if (eb > 0.2 * e * b) {
2512 ebang = asin(std::min(
2513 1., sqrt(pow(ex * by - ey * bx, 2) + pow(ex * bz - ez * bx, 2) +
2514 pow(ez * by - ey * bz, 2)) /
2515 (e * b)));
2516 } else {
2517 ebang = acos(std::min(1., eb / (e * b)));
2518 }
2519 } else {
2520 ebang = m_bAngles[0];
2521 }
2522 */
2523 // Hole
2524 /*
2525 if (e * b > 0.) {
2526 const double eb = fabs(ex * bx + ey * by + ez * bz);
2527 if (eb > 0.2 * e * b) {
2528 ebang = asin(std::min(
2529 1., sqrt(pow(ex * by - ey * bx, 2) + pow(ex * bz - ez * bx, 2) +
2530 pow(ez * by - ey * bz, 2)) /
2531 (e * b)));
2532 } else {
2533 ebang = acos(std::min(1., eb / (e * b)));
2534 }
2535 } else {
2536 ebang = m_bAngles[0];
2537 }
2538 */
2539 /*
2540 if (e * b > 0.) {
2541 const double eb = fabs(ex * bx + ey * by + ez * bz);
2542 if (eb > 0.2 * e * b) {
2543 ebang = asin(std::min(
2544 1., sqrt(pow(ex * by - ey * bx, 2) + pow(ex * bz - ez * bx, 2) +
2545 pow(ez * by - ey * bz, 2)) /
2546 (e * b)));
2547 } else {
2548 ebang = acos(std::min(1., eb / (e * b)));
2549 }
2550 } else {
2551 ebang = m_bAngles[0];
2552 }
2553 */
2554 if (e * b <= 0.) return m_bAngles[0];
2555 const double eb = fabs(ex * bx + ey * by + ez * bz);
2556 if (eb > 0.2 * e * b) {
2557 const double ebxy = ex * by - ey * bx;
2558 const double ebxz = ex * bz - ez * bx;
2559 const double ebzy = ez * by - ey * bz;
2560 return asin(std::min(1., sqrt(ebxy * ebxy + ebxz * ebxz + ebzy * ebzy) / (e * b)));
2561 }
2562 return acos(std::min(1., eb / (e * b)));
2563}
2564
2565double Medium::Interpolate1D(const double e, const std::vector<double>& table,
2566 const std::vector<double>& fields,
2567 const unsigned int intpMeth, const int extrLow,
2568 const int extrHigh) {
2569
2570 // This function is a generalized version of the Fortran functions
2571 // GASVEL, GASVT1, GASVT2, GASLOR, GASMOB, GASDFT, and GASDFL
2572 // for the case of a 1D table. All variables are generic.
2573
2574 const int nSizeTable = fields.size();
2575
2576 if (e < 0. || nSizeTable < 1) return 0.;
2577
2578 double result = 0.;
2579
2580 if (nSizeTable == 1) {
2581 // Only one point
2582 result = table[0];
2583 } else if (e < fields[0]) {
2584 // Extrapolation towards small fields
2585 if (fields[0] >= fields[1]) {
2586 if (m_debug) {
2587 std::cerr << m_className << "::Interpolate1D:\n";
2588 std::cerr << " First two field values coincide.\n";
2589 std::cerr << " No extrapolation to lower fields.\n";
2590 }
2591 result = table[0];
2592 } else if (extrLow == 1) {
2593 // Linear extrapolation
2594 const double extr4 = (table[1] - table[0]) / (fields[1] - fields[0]);
2595 const double extr3 = table[0] - extr4 * fields[0];
2596 result = extr3 + extr4 * e;
2597 } else if (extrLow == 2) {
2598 // Logarithmic extrapolation
2599 const double extr4 = log(table[1] / table[0]) / (fields[1] - fields[0]);
2600 const double extr3 = log(table[0] - extr4 * fields[0]);
2601 result = std::exp(std::min(50., extr3 + extr4 * e));
2602 } else {
2603 result = table[0];
2604 }
2605 } else if (e > fields[nSizeTable - 1]) {
2606 // Extrapolation towards large fields
2607 if (fields[nSizeTable - 1] <= fields[nSizeTable - 2]) {
2608 if (m_debug) {
2609 std::cerr << m_className << "::Interpolate1D:\n";
2610 std::cerr << " Last two field values coincide.\n";
2611 std::cerr << " No extrapolation to higher fields.\n";
2612 }
2613 result = table[nSizeTable - 1];
2614 } else if (extrHigh == 1) {
2615 // Linear extrapolation
2616 const double extr2 = (table[nSizeTable - 1] - table[nSizeTable - 2]) /
2617 (fields[nSizeTable - 1] - fields[nSizeTable - 2]);
2618 const double extr1 =
2619 table[nSizeTable - 1] - extr2 * fields[nSizeTable - 1];
2620 result = extr1 + extr2 * e;
2621 } else if (extrHigh == 2) {
2622 // Logarithmic extrapolation
2623 const double extr2 = log(table[nSizeTable - 1] / table[nSizeTable - 2]) /
2624 (fields[nSizeTable - 1] - fields[nSizeTable - 2]);
2625 const double extr1 =
2626 log(table[nSizeTable - 1]) - extr2 * fields[nSizeTable - 1];
2627 result = exp(std::min(50., extr1 + extr2 * e));
2628 } else {
2629 result = table[nSizeTable - 1];
2630 }
2631 } else {
2632 // Intermediate points, spline interpolation (not implemented).
2633 // Intermediate points, Newtonian interpolation
2634 result = Numerics::Divdif(table, fields, nSizeTable, e, intpMeth);
2635 }
2636
2637 return result;
2638}
2639
2641 const unsigned int eRes, const unsigned int bRes,
2642 const unsigned int aRes,
2643 std::vector<std::vector<std::vector<double> > >& tab, const double val) {
2644
2645 if (eRes == 0 || bRes == 0 || aRes == 0) {
2646 std::cerr << m_className << "::InitParamArrays:\n";
2647 std::cerr << " Invalid grid.\n";
2648 return;
2649 }
2650
2651 tab.assign(aRes, std::vector<std::vector<double> >(bRes, std::vector<double>(eRes, val)));
2652 /*
2653 tab.resize(aRes);
2654 for (unsigned int i = 0; i < aRes; ++i) {
2655 tab[i].resize(bRes);
2656 for (unsigned int j = 0; j < bRes; ++j) {
2657 tab[i][j].assign(eRes, val);
2658 }
2659 }
2660 */
2661}
2662
2664 const unsigned int eRes, const unsigned int bRes,
2665 const unsigned int aRes, const unsigned int tRes,
2666 std::vector<std::vector<std::vector<std::vector<double> > > >& tab,
2667 const double val) {
2668
2669 if (eRes == 0 || bRes == 0 || aRes == 0 || tRes == 0) {
2670 std::cerr << m_className << "::InitParamArrays: Invalid grid.\n";
2671 return;
2672 }
2673
2674 tab.resize(tRes);
2675 for (unsigned int l = 0; l < tRes; ++l) {
2676 tab[l].resize(aRes);
2677 for (unsigned int i = 0; i < aRes; ++i) {
2678 tab[l][i].resize(bRes);
2679 for (unsigned int j = 0; j < bRes; ++j) {
2680 tab[l][i][j].assign(eRes, val);
2681 }
2682 }
2683 }
2684}
2685}
void ResetHoleAttachment()
Definition: Medium.cc:1369
virtual void GetComponent(const unsigned int i, std::string &label, double &f)
Definition: Medium.cc:142
double m_density
Definition: Medium.hh:315
bool m_hasElectronDiffTens
Definition: Medium.hh:340
void ResetIonMobility()
Definition: Medium.cc:1375
virtual bool GetPhotonCollision(const double e, int &type, int &level, double &e1, double &ctheta, int &nsec, double &esec)
Definition: Medium.cc:1294
bool SetIonMobility(const unsigned int ie, const unsigned int ib, const unsigned int ia, const double mu)
Definition: Medium.cc:2269
void ResetElectronAttachment()
Definition: Medium.cc:1331
virtual double GetMassDensity() const
Definition: Medium.cc:137
unsigned int m_extrLowMobility
Definition: Medium.hh:392
std::vector< double > m_bFields
Definition: Medium.hh:333
void ResetIonDiffusion()
Definition: Medium.cc:1381
bool GetIonDissociation(const unsigned int ie, const unsigned int ib, const unsigned int ia, double &diss)
Definition: Medium.cc:2113
virtual void GetElectronMomentum(const double e, double &px, double &py, double &pz, int &band)
Definition: Medium.cc:651
void SetTemperature(const double t)
Definition: Medium.cc:104
unsigned int m_intpDiffusion
Definition: Medium.hh:397
bool GetElectronVelocityB(const unsigned int ie, const unsigned int ib, const unsigned int ia, double &v)
Definition: Medium.cc:1728
std::vector< std::vector< std::vector< double > > > tabHoleVelocityE
Definition: Medium.hh:359
unsigned int m_extrLowDiffusion
Definition: Medium.hh:388
bool GetElectronVelocityE(const unsigned int ie, const unsigned int ib, const unsigned int ia, double &v)
Definition: Medium.cc:1680
void ResetElectronLorentzAngle()
Definition: Medium.cc:1337
double m_pressure
Definition: Medium.hh:305
void CloneTensor(std::vector< std::vector< std::vector< std::vector< double > > > > &tab, const unsigned int n, const std::vector< double > &efields, const std::vector< double > &bfields, const std::vector< double > &angles, const unsigned int intp, const unsigned int extrLow, const unsigned int extrHigh, const double init, const std::string &label)
Definition: Medium.cc:2203
virtual bool HoleTownsend(const double ex, const double ey, const double ez, const double bx, const double by, const double bz, double &alpha)
Definition: Medium.cc:971
bool m_hasElectronDiffTrans
Definition: Medium.hh:340
virtual double ScaleAttachment(const double eta) const
Definition: Medium.hh:270
unsigned int m_extrLowLorentzAngle
Definition: Medium.hh:391
virtual bool GetElectronCollision(const double e, int &type, int &level, double &e1, double &dx, double &dy, double &dz, int &nion, int &ndxc, int &band)
Definition: Medium.cc:672
void SetExtrapolationMethodIonDissociation(const std::string &extrLow, const std::string &extrHigh)
Definition: Medium.cc:2434
void SetInterpolationMethodIonDissociation(const unsigned int intrp)
Definition: Medium.cc:2498
double GetAngle(const double ex, const double ey, const double ez, const double bx, const double by, const double bz, const double e, const double b) const
Definition: Medium.cc:2503
std::vector< std::vector< std::vector< double > > > tabElectronVelocityB
Definition: Medium.hh:345
unsigned int m_extrHighAttachment
Definition: Medium.hh:390
virtual double ScaleDiffusion(const double d) const
Definition: Medium.hh:267
virtual bool HoleVelocity(const double ex, const double ey, const double ez, const double bx, const double by, const double bz, double &vx, double &vy, double &vz)
Definition: Medium.cc:704
unsigned int m_intpMobility
Definition: Medium.hh:401
virtual double ScaleLorentzAngle(const double lor) const
Definition: Medium.hh:271
unsigned int m_intpAttachment
Definition: Medium.hh:399
virtual bool ElectronLorentzAngle(const double ex, const double ey, const double ez, const double bx, const double by, const double bz, double &lor)
Definition: Medium.cc:599
bool GetIonLongitudinalDiffusion(const unsigned int ie, const unsigned int ib, const unsigned int ia, double &dl)
Definition: Medium.cc:2065
bool m_hasElectronVelocityB
Definition: Medium.hh:339
bool GetHoleAttachment(const unsigned int ie, const unsigned int ib, const unsigned int ia, double &eta)
Definition: Medium.cc:2018
bool GetHoleVelocityExB(const unsigned int ie, const unsigned int ib, const unsigned int ia, double &v)
Definition: Medium.cc:1898
void ResetHoleVelocity()
Definition: Medium.cc:1343
bool GetExtrapolationIndex(std::string extrStr, unsigned int &extrNb)
Definition: Medium.cc:2453
void ResetIonDissociation()
Definition: Medium.cc:1389
bool GetElectronAttachment(const unsigned int ie, const unsigned int ib, const unsigned int ia, double &eta)
Definition: Medium.cc:1826
unsigned int m_extrHighDiffusion
Definition: Medium.hh:388
virtual bool GetDielectricFunction(const double e, double &eps1, double &eps2, const unsigned int i=0)
Definition: Medium.cc:1243
std::vector< std::vector< std::vector< double > > > tabIonDissociation
Definition: Medium.hh:376
unsigned int m_extrHighMobility
Definition: Medium.hh:392
std::vector< std::vector< std::vector< double > > > tabElectronVelocityE
Definition: Medium.hh:343
void SetFieldGrid(double emin, double emax, int ne, bool logE, double bmin=0., double bmax=0., int nb=1, double amin=0., double amax=0., int na=1)
Definition: Medium.cc:1395
virtual void SetNumberDensity(const double n)
Definition: Medium.cc:176
std::string m_name
Definition: Medium.hh:301
unsigned int m_intpTownsend
Definition: Medium.hh:398
unsigned int m_extrLowAttachment
Definition: Medium.hh:390
unsigned int m_extrHighDissociation
Definition: Medium.hh:393
void InitParamArrays(const unsigned int eRes, const unsigned int bRes, const unsigned int aRes, std::vector< std::vector< std::vector< double > > > &tab, const double val)
Definition: Medium.cc:2640
double Interpolate1D(const double e, const std::vector< double > &table, const std::vector< double > &fields, const unsigned int intpMeth, const int jExtr, const int iExtr)
Definition: Medium.cc:2565
virtual bool ElectronVelocity(const double ex, const double ey, const double ez, const double bx, const double by, const double bz, double &vx, double &vy, double &vz)
Definition: Medium.cc:204
virtual bool IonVelocity(const double ex, const double ey, const double ez, const double bx, const double by, const double bz, double &vx, double &vy, double &vz)
Definition: Medium.cc:1077
bool GetIonMobility(const unsigned int ie, const unsigned int ib, const unsigned int ia, double &mu)
Definition: Medium.cc:2042
unsigned int m_intpVelocity
Definition: Medium.hh:396
void SetExtrapolationMethodTownsend(const std::string &extrLow, const std::string &extrHigh)
Definition: Medium.cc:2378
bool m_hasHoleTownsend
Definition: Medium.hh:358
std::vector< std::vector< std::vector< double > > > tabElectronDiffLong
Definition: Medium.hh:346
virtual void SetAtomicNumber(const double z)
Definition: Medium.cc:154
std::vector< std::vector< std::vector< double > > > tabElectronLorentzAngle
Definition: Medium.hh:350
unsigned int m_extrHighLorentzAngle
Definition: Medium.hh:391
void SetExtrapolationMethodIonMobility(const std::string &extrLow, const std::string &extrHigh)
Definition: Medium.cc:2416
virtual double ScaleElectricField(const double e) const
Definition: Medium.hh:264
std::vector< std::vector< std::vector< double > > > tabElectronAttachment
Definition: Medium.hh:349
int thrIonDissociation
Definition: Medium.hh:384
virtual double ScaleDiffusionTensor(const double d) const
Definition: Medium.hh:268
virtual double GetElectronNullCollisionRate(const int band=0)
Definition: Medium.cc:659
bool GetElectronTownsend(const unsigned int ie, const unsigned int ib, const unsigned int ia, double &alpha)
Definition: Medium.cc:1802
bool GetElectronLorentzAngle(const unsigned int ie, const unsigned int ib, const unsigned int ia, double &lor)
Definition: Medium.cc:1850
virtual bool GetIonisationProduct(const unsigned int i, int &type, double &energy) const
Definition: Medium.cc:685
unsigned int m_extrLowVelocity
Definition: Medium.hh:387
virtual double ScaleTownsend(const double alpha) const
Definition: Medium.hh:269
virtual double GetPhotonCollisionRate(const double e)
Definition: Medium.cc:1286
void SetDielectricConstant(const double eps)
Definition: Medium.cc:126
unsigned int m_extrHighTownsend
Definition: Medium.hh:389
virtual bool GetDeexcitationProduct(const unsigned int i, double &t, double &s, int &type, double &energy) const
Definition: Medium.cc:694
virtual double GetElectronCollisionRate(const double e, const int band=0)
Definition: Medium.cc:665
bool GetHoleTownsend(const unsigned int ie, const unsigned int ib, const unsigned int ia, double &alpha)
Definition: Medium.cc:1994
void SetInterpolationMethodVelocity(const unsigned int intrp)
Definition: Medium.cc:2473
virtual bool ElectronDiffusion(const double ex, const double ey, const double ez, const double bx, const double by, const double bz, double &dl, double &dt)
Definition: Medium.cc:362
unsigned int m_extrHighVelocity
Definition: Medium.hh:387
void SetInterpolationMethodDiffusion(const unsigned int intrp)
Definition: Medium.cc:2478
std::vector< std::vector< std::vector< double > > > tabHoleTownsend
Definition: Medium.hh:364
double m_epsilon
Definition: Medium.hh:307
unsigned int m_extrLowTownsend
Definition: Medium.hh:389
void SetExtrapolationMethodDiffusion(const std::string &extrLow, const std::string &extrHigh)
Definition: Medium.cc:2359
std::vector< double > m_eFields
Definition: Medium.hh:332
void GetFieldGrid(std::vector< double > &efields, std::vector< double > &bfields, std::vector< double > &angles)
Definition: Medium.cc:1671
static int m_idCounter
Definition: Medium.hh:296
std::vector< std::vector< std::vector< std::vector< double > > > > tabHoleDiffTens
Definition: Medium.hh:367
virtual bool ElectronTownsend(const double ex, const double ey, const double ez, const double bx, const double by, const double bz, double &alpha)
Definition: Medium.cc:490
void ResetElectronDiffusion()
Definition: Medium.cc:1316
std::vector< std::vector< std::vector< double > > > tabIonDiffTrans
Definition: Medium.hh:375
virtual bool GetOpticalDataRange(double &emin, double &emax, const unsigned int i=0)
Definition: Medium.cc:1229
virtual void SetMassDensity(const double rho)
Definition: Medium.cc:187
void CloneTable(std::vector< std::vector< std::vector< double > > > &tab, const std::vector< double > &efields, const std::vector< double > &bfields, const std::vector< double > &angles, const unsigned int intp, const unsigned int extrLow, const unsigned int extrHigh, const double init, const std::string &label)
Definition: Medium.cc:2137
virtual double ScaleDissociation(const double diss) const
Definition: Medium.hh:272
bool GetElectronVelocityExB(const unsigned int ie, const unsigned int ib, const unsigned int ia, double &v)
Definition: Medium.cc:1704
unsigned int m_extrLowDissociation
Definition: Medium.hh:393
bool m_hasHoleVelocityB
Definition: Medium.hh:356
std::vector< std::vector< std::vector< std::vector< double > > > > tabElectronDiffTens
Definition: Medium.hh:353
std::vector< double > m_bAngles
Definition: Medium.hh:334
void InitParamTensor(const unsigned int eRes, const unsigned int bRes, const unsigned int aRes, const unsigned int tRes, std::vector< std::vector< std::vector< std::vector< double > > > > &tab, const double val)
Definition: Medium.cc:2663
bool GetElectronTransverseDiffusion(const unsigned int ie, const unsigned int ib, const unsigned int ia, double &dt)
Definition: Medium.cc:1777
virtual double GetElectronEnergy(const double px, const double py, const double pz, double &vx, double &vy, double &vz, const int band=0)
Definition: Medium.cc:635
void SetInterpolationMethodIonMobility(const unsigned int intrp)
Definition: Medium.cc:2493
std::vector< std::vector< std::vector< double > > > tabHoleVelocityExB
Definition: Medium.hh:360
void SetExtrapolationMethodVelocity(const std::string &extrLow, const std::string &extrHigh)
Definition: Medium.cc:2340
virtual bool ElectronAttachment(const double ex, const double ey, const double ez, const double bx, const double by, const double bz, double &eta)
Definition: Medium.cc:544
bool m_hasElectronVelocityExB
Definition: Medium.hh:339
std::vector< std::vector< std::vector< double > > > tabIonMobility
Definition: Medium.hh:373
bool m_hasElectronLorentzAngle
Definition: Medium.hh:342
void SetPressure(const double p)
Definition: Medium.cc:115
bool m_hasHoleDiffTrans
Definition: Medium.hh:357
bool m_hasIonDiffTrans
Definition: Medium.hh:371
bool GetHoleVelocityE(const unsigned int ie, const unsigned int ib, const unsigned int ia, double &v)
Definition: Medium.cc:1874
int thrElectronAttachment
Definition: Medium.hh:380
bool m_hasIonMobility
Definition: Medium.hh:370
bool GetHoleVelocityB(const unsigned int ie, const unsigned int ib, const unsigned int ia, double &v)
Definition: Medium.cc:1922
bool m_hasElectronDiffLong
Definition: Medium.hh:340
std::vector< std::vector< std::vector< double > > > tabHoleDiffTrans
Definition: Medium.hh:363
virtual void SetAtomicWeight(const double a)
Definition: Medium.cc:165
void SetExtrapolationMethodAttachment(const std::string &extrLow, const std::string &extrHigh)
Definition: Medium.cc:2397
unsigned int m_intpLorentzAngle
Definition: Medium.hh:400
std::vector< std::vector< std::vector< double > > > tabHoleVelocityB
Definition: Medium.hh:361
bool m_hasHoleDiffLong
Definition: Medium.hh:357
std::vector< std::vector< std::vector< double > > > tabHoleDiffLong
Definition: Medium.hh:362
unsigned int m_nComponents
Definition: Medium.hh:309
int thrHoleAttachment
Definition: Medium.hh:383
std::vector< std::vector< std::vector< double > > > tabIonDiffLong
Definition: Medium.hh:374
std::string m_className
Definition: Medium.hh:294
bool GetHoleTransverseDiffusion(const unsigned int ie, const unsigned int ib, const unsigned int ia, double &dt)
Definition: Medium.cc:1970
bool GetElectronLongitudinalDiffusion(const unsigned int ie, const unsigned int ib, const unsigned int ia, double &dl)
Definition: Medium.cc:1752
void SetInterpolationMethodAttachment(const unsigned int intrp)
Definition: Medium.cc:2488
virtual bool HoleDiffusion(const double ex, const double ey, const double ez, const double bx, const double by, const double bz, double &dl, double &dt)
Definition: Medium.cc:846
bool GetHoleLongitudinalDiffusion(const unsigned int ie, const unsigned int ib, const unsigned int ia, double &dl)
Definition: Medium.cc:1946
std::vector< std::vector< std::vector< double > > > tabHoleAttachment
Definition: Medium.hh:365
virtual bool IonDissociation(const double ex, const double ey, const double ez, const double bx, const double by, const double bz, double &diss)
Definition: Medium.cc:1176
bool m_hasIonDissociation
Definition: Medium.hh:372
bool m_hasIonDiffLong
Definition: Medium.hh:371
int thrHoleTownsend
Definition: Medium.hh:382
unsigned int m_intpDissociation
Definition: Medium.hh:402
virtual ~Medium()
Definition: Medium.cc:102
bool m_hasHoleAttachment
Definition: Medium.hh:358
void ResetHoleDiffusion()
Definition: Medium.cc:1353
bool m_hasHoleVelocityExB
Definition: Medium.hh:356
bool m_isChanged
Definition: Medium.hh:326
bool m_hasElectronVelocityE
Definition: Medium.hh:339
std::vector< std::vector< std::vector< double > > > tabElectronTownsend
Definition: Medium.hh:348
bool m_hasHoleVelocityE
Definition: Medium.hh:356
int thrElectronTownsend
Definition: Medium.hh:379
bool m_hasElectronAttachment
Definition: Medium.hh:341
void ResetElectronTownsend()
Definition: Medium.cc:1326
bool m_hasHoleDiffTens
Definition: Medium.hh:357
virtual bool HoleAttachment(const double ex, const double ey, const double ez, const double bx, const double by, const double bz, double &eta)
Definition: Medium.cc:1023
void ResetHoleTownsend()
Definition: Medium.cc:1363
void ResetElectronVelocity()
Definition: Medium.cc:1306
void SetInterpolationMethodTownsend(const unsigned int intrp)
Definition: Medium.cc:2483
virtual bool GetPhotoAbsorptionCrossSection(const double e, double &sigma, const unsigned int i=0)
Definition: Medium.cc:1264
virtual bool IonDiffusion(const double ex, const double ey, const double ez, const double bx, const double by, const double bz, double &dl, double &dt)
Definition: Medium.cc:1122
double m_temperature
Definition: Medium.hh:303
std::vector< std::vector< std::vector< double > > > tabElectronDiffTrans
Definition: Medium.hh:347
std::vector< std::vector< std::vector< double > > > tabElectronVelocityExB
Definition: Medium.hh:344
bool GetIonTransverseDiffusion(const unsigned int ie, const unsigned int ib, const unsigned int ia, double &dt)
Definition: Medium.cc:2089
bool Boxin3(const std::vector< std::vector< std::vector< double > > > &value, const std::vector< double > &xAxis, const std::vector< double > &yAxis, const std::vector< double > &zAxis, const int nx, const int ny, const int nz, const double xx, const double yy, const double zz, double &f, const int iOrder)
Definition: Numerics.cc:770
double Divdif(const std::vector< double > &f, const std::vector< double > &a, int nn, double x, int mm)
Definition: Numerics.cc:650
void RndmDirection(double &dx, double &dy, double &dz, const double length=1.)
Draw a random (isotropic) direction vector.
Definition: Random.hh:106