OPAL (Object Oriented Parallel Accelerator Library)  2.2.0
OPAL
SecondaryEmissionPhysics.cpp
Go to the documentation of this file.
1 #include "Structure/SecondaryEmissionPhysics.h"
2 #include <vector>
3 #include <cmath>
4 #include <algorithm>
5 #include <sys/stat.h>
6 #include "Utilities/Options.h"
7 #include "Utility/IpplTimings.h"
8 
9 using namespace myeps;
10 using namespace Physics;
11 SecondaryEmissionPhysics::SecondaryEmissionPhysics() {
12 
13  TPnSec_m = IpplTimings::getTimer("Secondary emission");
14 }
20 SecondaryEmissionPhysics::~SecondaryEmissionPhysics() {
21 
22 }
23 
24 void SecondaryEmissionPhysics::nSec(const double &incEnergy,
25  const double &cosTheta,
26  const int &matNumber,
27  int &seNum,
28  int &seType,
29  const double &incQ,
30  const Vector_t &TriNorm,
31  const Vector_t &inteCoords,
32  const Vector_t &localX,
33  PartBunchBase<double, 3> *itsBunch,
34  double &seyNum,
35  const double &ppVw,
36  const double &vVThermal,
37  const bool nEmissionMode) {
38 
39  IpplTimings::startTimer(TPnSec_m);
40  double prob[11] = {0};
41  std::vector<Vector_t> se_P;
42  setSeMaterial(matNumber);//set material based parameters
43  seyNum=calcProb(incEnergy, cosTheta, prob);//calculate probability
44  calcEmiNum(incEnergy, cosTheta, prob, seNum);//calculate emitted number
45  PAssert_LT(seNum, 11);//debug
46  PAssert_GE(seNum, 0);//debug
47  double Eemit[10];
48  double emiTheta[10];
49  double emiPhi[10];
50  Vector_t interCoords_l = inteCoords;
51  Vector_t TriNorm_l = TriNorm;
52  double incQ_l = incQ;
53 
54 
55  /*===========================Definitions for benchmark===================================*/
56  double vw=ppVw; //1.6*1e-19*1200/9.10938188*1e-31/(2*3.1415926*2.0*1e8)/0.03;//benchmark
57  double vt=vVThermal;//7.268929821*1e5;//1.5eV//benchmark
58  double f_max=vw/vt*exp(-0.5);//benchmark
59  double test_a=vt/vw;//benchmark
60  double test_asq=test_a*test_a;//benchmark
61  /*---------------------------------------------------------------------------------------*/
62  if( Options::ppdebug ) {
63 
64  } else {
65 
66  if(seNum != 0) {
67 
68  for(int i = 0; i < seNum; i++) {
69 
70  double tmp1 = IpplRandom();
71  double tmp2 = IpplRandom();
72  double temp = 1.0 / (1.0 + seAlpha_m);
73  emiTheta[i] = acos(pow(tmp1, temp));
74  emiPhi[i] = Physics::two_pi * tmp2;
75 
76  }
77  }
78  }
79 
80 
81  if(seNum == 0) {
82 
83  // The incident particle will be marked for deletion
84  if (!nEmissionMode) {
85  if( Options::ppdebug ) {
86  /*=========(Velocity with Maxwellian Distribution For Parallel Plate Benchmark)========*/
87  double test_s=1;
88  double f_x=0;
89  double test_x=0;
90  while (test_s>f_x) {
91  test_s=IpplRandom();
92  test_s*=f_max;
93  test_x=IpplRandom();
94  test_x*=10*test_a;// range for normalized emission speed(0,10*test_a);
95  f_x=test_x/test_asq*exp(-test_x*test_x/2/test_asq);
96  }
97  double v_emi=test_x*vw;
98  Eemit[0]=(1.0/sqrt(1.0-v_emi*v_emi/9.0/1e16)-1)*Physics::m_e*1.0e9;
99  // cout<<"Single Eemit[0]: "<<Eemit[0]<<endl;
100  /*---------------------End Of Maxwellian Distribution(For Benchmark)--------------------*/
101  } else {
102 
103  // For absorption case in constant simulation particle mode, just use true secondary emission with 1 particle energy distribution for Furman Pivi model.
104  double u2 = IpplRandom();
105 
106  double temp = incEnergy / seEpsn_m[0];
107  double p0 = gammp(sePn_m[0], temp);
108  temp = p0 * u2;
109  Eemit[0] = seEpsn_m[0] * invgammp(temp, sePn_m[0]) ;
110 
111 
112  }
113 
114  }
115 
116 
117  } else if(seNum == 1) {
118 
119  if( Options::ppdebug ) {
120  /*=========(Velocity with Maxwellian Distribution For Parallel Plate Benchmark)========*/
121  double test_s=1;
122  double f_x=0;
123  double test_x=0;
124  while (test_s>f_x) {
125  test_s=IpplRandom();
126  test_s*=f_max;
127  test_x=IpplRandom();
128  test_x*=10*test_a;//range for normalized emission speed(0,10*test_a);
129  f_x=test_x/test_asq*exp(-test_x*test_x/2/test_asq);
130  }
131  double v_emi=test_x*vw;
132  Eemit[0]=(1.0/sqrt(1.0-v_emi*v_emi/9.0/1e16)-1)*Physics::m_e*1.0e9;
133  //cout<<"Single Eemit[0]: "<<Eemit[0]<<endl;
134  /*---------------------End Of Maxwellian Distribution(For Benchmark)--------------------*/
135  } else {
136 
137  //double delta_e = calcDeltae(incEnergy, cosTheta);
138  //double delta_r = calcDeltar(incEnergy, cosTheta);
139  double tmp = prob[1] + deltae_m + deltar_m;
140  //double ae = delta_e / tmp;
141  double ae = deltae_m / tmp;
142  //double ar = delta_r / tmp;
143  double ar = deltar_m / tmp;
144  double a_re = ae + ar;
145  double urand = IpplRandom();
146 
147 
148  if(urand < ae) {
149  int t_count = 0;
150  do {
151  Eemit[0] = incEnergy - seDelta_m * fabs(gaussRand()) ;
152  t_count++;
153  } while(Eemit[0] < 0&&t_count<200);
154  if(Eemit[0]<0)// if the above do - while loops over 200 times, the loop will break out, and Eemit will be its mean value, i.e., incident energy.
155  Eemit[0]=incEnergy;
156  seType = 0;
157 
158  } else if(urand >= ae && urand < a_re) {
159  double u1 = IpplRandom();
160 
161  double powArg = 1.0 / (1.0 + seQ_m);
162  Eemit[0] = incEnergy * pow(u1, powArg);
163  seType = 1;
164 
165 
166  } else {
167 
168  double u2 = IpplRandom();
169 
170  double temp = incEnergy / seEpsn_m[0];
171  double p0 = gammp(sePn_m[0], temp);
172  temp = p0 * u2;
173  Eemit[0] = seEpsn_m[0] * invgammp(temp, sePn_m[0]) ;
174  seType = 2;
175 
176  }
177  }
178 
179  } else {
180  seType = 2;
181 
182  if( Options::ppdebug ) {
183  /*==========(Velocity with Maxwellian Distribution For Parallel Plate Benchmark)========*/
184  if (!nEmissionMode) {
185  /*double Eemit_mean = 0.0;
186  for(int i = 0; i < seNum; i++) {
187  Eemit_mean += Eemit[i];
188  }*/
189  //Eemit[0] = Eemit_mean/seNum;
190  double test_s=1.0;
191  double f_x=0;
192  double test_x=0;
193  while (test_s>f_x) {
194  test_s=IpplRandom();
195  test_s*=f_max;
196  test_x=IpplRandom();
197  test_x*=10*test_a;//range for normalized emission speed(0,10*test_a);
198  f_x=test_x/test_asq*exp(-test_x*test_x/2/test_asq);
199  }
200  double v_emi=test_x*vw;
201  Eemit[0]=(1.0/sqrt(1.0-v_emi*v_emi/9.0/1e16)-1)*Physics::m_e*1.0e9;
202 
203 
204  } else {
205  for(int i = 0; i < seNum; i++) {
206  double test_s=1.0;
207  double f_x=0;
208  double test_x=0;
209  while (test_s>f_x) {
210  test_s=IpplRandom();
211  test_s*=f_max;
212  test_x=IpplRandom();
213  test_x*=10*test_a;//range for normalized emission speed(0,10*test_a);
214  f_x=test_x/test_asq*exp(-test_x*test_x/2/test_asq);
215  }
216  double v_emi=test_x*vw;
217  Eemit[i]=(1.0/sqrt(1.0-v_emi*v_emi/9.0/1e16)-1)*Physics::m_e*1.0e9;
218 
219  }
220 
221 
222  }
223  /*---------------------End Of Maxwellian Distribution(For Benchmark)--------------------*/
224  } else {
225  /*=======================3D velocity for Furman-Pivi's Model====================*/
226  double x0 = incEnergy / seEpsn_m[seNum-1];
227  double parg = seNum * sePn_m[seNum-1];
228  double p0 = gammp(parg, x0);
229  double sin2theta_n[seNum];
230  double cos2theta_n[seNum];
231  double rand_y = IpplRandom();
232  double invarg = rand_y * p0;
233  double y2 = invgammp(invarg, parg);
234  double y2_n[seNum];
235  double multisin = 1.0;
236 
237  if (!nEmissionMode) {// only emit 1 particle
238 
239  double Eemisum = 0.0;
240  for(int i = 0; i < seNum - 1; i++) {
241  double mu = sePn_m[seNum-1] * (seNum - i);
242  double nu = sePn_m[seNum-1];
243  double rand_n = IpplRandom();
244  sin2theta_n[i] = invbetai(rand_n, mu, nu);
245  cos2theta_n[i] = 1.0 - sin2theta_n[i];
246  if(i != 0) {
247 
248  multisin *= sin2theta_n[i-1];
249  }
250  y2_n[i] = y2 * multisin * cos2theta_n[i];
251  Eemit[i] = seEpsn_m[seNum-1] * y2_n[i];
252  Eemisum += Eemit[i];
253 
254 
255  }
256 
257  Eemit[seNum-1] = Eemit[seNum-2] / cos2theta_n[seNum-2] * sin2theta_n[seNum-2];
258  Eemisum += Eemit[seNum-1];
259 
260  Eemit[0] = Eemisum/seNum;
261 
262  } else {// emit seNum particles
263  for(int i = 0; i < seNum - 1; i++) {
264  double mu = sePn_m[seNum-1] * (seNum - i);
265  double nu = sePn_m[seNum-1];
266  double rand_n = IpplRandom();
267  sin2theta_n[i] = invbetai(rand_n, mu, nu);
268  cos2theta_n[i] = 1.0 - sin2theta_n[i];
269  if(i != 0) {
270 
271  multisin *= sin2theta_n[i-1];
272  }
273  y2_n[i] = y2 * multisin * cos2theta_n[i];
274  Eemit[i] = seEpsn_m[seNum-1] * y2_n[i];
275 
276 
277  }
278 
279  Eemit[seNum-1] = Eemit[seNum-2] / cos2theta_n[seNum-2] * sin2theta_n[seNum-2];
280 
281  /*=====================================================================================================*/
282  }
283  }
284 
285  }
286 
287  Vector_t z_unit = TriNorm_l;
288  Vector_t x_unit = localX - interCoords_l;
289  double tmp = sqrt(dot(x_unit,x_unit));
290  x_unit /= tmp;
291  Vector_t y_unit = cross(z_unit,x_unit);
292 
293  size_t lowMark = itsBunch->getLocalNum();
294  if (!nEmissionMode) {// emit only 1 particle with larger charge instead of emit seNum secondaries.
295  //if (seNum>0) {// we dont delete particles even when seNum==0.
296  double gamma_const = Eemit[0] / Physics::m_e/1.0e9 + 1.0;
297  double beta_const = sqrt(1.0 - 1.0 / pow(gamma_const, 2.0));
298  double P_emitted = gamma_const * beta_const;
299 
300  Vector_t P_local;
301  Vector_t P_global = (0.0);
302 
303  if( Options::ppdebug ) {
304 
305  P_global = P_emitted*TriNorm_l;// 1D for parallel plate benchmark
306 
307  } else {
308  /*==================3D for Furman-Pivi's Model====================*/
309  P_local[2] = P_emitted * cos(emiTheta[0]);
310  P_local[1] = P_emitted * sin(emiTheta[0]) * sin(emiPhi[0]);
311  P_local[0] = P_emitted * sin(emiTheta[0]) * cos(emiPhi[0]);
312  P_global = P_local[0]*x_unit + P_local[1]*y_unit + P_local[2]*z_unit;//Pivi's model
313  /*================================================================*/
314 
315 
316  }
317  itsBunch->create(1);
318  itsBunch->R[lowMark] = interCoords_l;
319 
320  itsBunch->P[lowMark] = P_global;
321  itsBunch->Bin[lowMark] = 0;
322  itsBunch->PType[lowMark] = ParticleType::NEWSECONDARY;
323  itsBunch->TriID[lowMark] = 0;
324  //itsBunch->Q[lowMark] = incQ_l*seNum;// charge of simulation particle will be sum of secondaies
325  itsBunch->Q[lowMark] = incQ_l*seyNum;// charge of simulation particle will be multiplied by SEY.
326  itsBunch->Ef[lowMark] = Vector_t(0.0);
327  itsBunch->Bf[lowMark] = Vector_t(0.0);
328  itsBunch->dt[lowMark] = itsBunch->getdT();
329  //}
330 
331 
332 
333  } else {
334  for(size_t i = 0; i < (size_t) seNum; i++) {
335 
336  double gamma_const = Eemit[i] / Physics::m_e/1.0e9 + 1.0;
337  double beta_const = sqrt(1.0 - 1.0 / pow(gamma_const, 2.0));
338  double P_emitted = gamma_const * beta_const;
339 
340  Vector_t P_local;
341  Vector_t P_global = (0.0);
342 
343  if( Options::ppdebug ) {
344 
345  P_global = P_emitted*TriNorm_l;// 1D for parallel plate benchmark
346 
347  } else {
348  /*==================3D for Furman-Pivi's Model====================*/
349  P_local[2] = P_emitted * cos(emiTheta[i]);
350  P_local[1] = P_emitted * sin(emiTheta[i]) * sin(emiPhi[i]);
351  P_local[0] = P_emitted * sin(emiTheta[i]) * cos(emiPhi[i]);
352  P_global = P_local[0]*x_unit + P_local[1]*y_unit + P_local[2]*z_unit;//Pivi's model
353  /*================================================================*/
354 
355 
356  }
357  itsBunch->create(1);
358  itsBunch->R[lowMark+i] = interCoords_l;
359 
360  itsBunch->P[lowMark+i] = P_global;
361  itsBunch->Bin[lowMark+i] = 0;
362  itsBunch->PType[lowMark+i] = ParticleType::NEWSECONDARY;
363  itsBunch->TriID[lowMark+i] = 0;
364  itsBunch->Q[lowMark+i] = incQ_l;
365  itsBunch->Ef[lowMark+i] = Vector_t(0.0);
366  itsBunch->Bf[lowMark+i] = Vector_t(0.0);
367  itsBunch->dt[lowMark+i] = itsBunch->getdT();
368 
369  }
370  }
371 
372  IpplTimings::stopTimer(TPnSec_m);
373 }
374 
375 
376 //Vaughan's secondary emission model.
377 void SecondaryEmissionPhysics::nSec(const double &incEnergy,
378  const double &cosTheta,
379  int &seNum, int &seType,
380  const double &incQ,
381  const Vector_t &TriNorm,
382  const Vector_t &inteCoords,
383  const Vector_t &localX,
384  PartBunchBase<double, 3> *itsBunch,
385  double &seyNum, const double &ppVw,
386  const double &vSeyZero,
387  const double &vEzero,
388  const double &vSeyMax,
389  const double &vEmax,
390  const double &vKenergy,
391  const double &vKtheta,
392  const double &vVThermal,
393  const bool nEmissionMode)
394 {
395 
396 
397  IpplTimings::startTimer(TPnSec_m);
398 
399  std::vector<Vector_t> se_P;
400  calcEmiNum(incEnergy, cosTheta, seNum, vSeyZero, vEzero, vSeyMax, vEmax, vKenergy, vKtheta, seyNum);//calculate emitted number & SEY factor
401  double Eemit[seNum];
402  double emiTheta[seNum];
403  double emiPhi[seNum];
404  Vector_t interCoords_l = inteCoords;
405  Vector_t TriNorm_l = TriNorm;//simpler model
406  double vw=ppVw; //1.6*1e-19*1200/9.10938188*1e-31/(2*3.1415926*2.0*1e8)/0.03;//benchmark
407  double vt=vVThermal;//7.268929821*1e5 1.5eV//benchmark
408  double f_max;//benchmark
409  double test_a;//benchmark
410  double test_asq;//benchmark
411  if( Options::ppdebug ) {
412  f_max=vw/vt*exp(-0.5);// velocity a Maxwellian distribution. See Anza et al.,Phys. Plasmas 17, 062110 (2010)
413 
414  test_a=vt/vw;
415  test_asq=test_a*test_a;
416  } else {// Energy Maxwell-Boltzmann distribution f(E)=E/a^2*exp(-E/a), a= kinetic energy w.r.t thermal speed vt. See www.nuc.berkeley.edu/dept/Courses/NE-255/minicourse.pdf p.20
417  test_a = Physics::m_e*(1.0/sqrt(1-vt*vt/Physics::c/Physics::c)-1.0)*1.0e9;// m_e GeV change it to eV
418  test_asq=test_a*test_a;
419  f_max= 1.0/test_a*exp(-1.0);
420 
421  }
422 
423  double incQ_l = incQ;
424  if( Options::ppdebug ) {
425  // 1D emission angle along the surface triangle normal
426 
427  } else {
428  if (!nEmissionMode) {
429 
430  double tmp1 = IpplRandom();
431  double tmp2 = IpplRandom();
432  double seAlpha = 1.0;
433  double temp = 1.0 / (1.0 + seAlpha);// pow(cosine(theta),seAlpha) distribution. Here seAlpha=1.0
434  emiTheta[0] = acos(pow(tmp1, temp));
435  emiPhi[0] = Physics::two_pi * tmp2;
436 
437 
438  } else {
439 
440  for(int i = 0; i < seNum; i++) {
441 
442  double tmp1 = IpplRandom();
443  double tmp2 = IpplRandom();
444  double seAlpha = 1.0;
445  double temp = 1.0 / (1.0 + seAlpha);// pow(cosine(theta),seAlpha) distribution. Here seAlpha=1.0
446  emiTheta[i] = acos(pow(tmp1, temp));
447  emiPhi[i] = Physics::two_pi * tmp2;
448 
449  }
450 
451  }
452  }
453 
454  if(seNum == 0) {
455  if (!nEmissionMode) {
456  if( Options::ppdebug ) {
457  double test_s=1.0;
458  double f_x=0;
459  double test_x=0;
460  while (test_s>f_x) {
461  test_s=IpplRandom();
462  test_s*=f_max;
463  test_x=IpplRandom();
464  test_x*=10*test_a;// truncation range for emission speed(0,10*test_a);
465  f_x=test_x/test_asq*exp(-test_x*test_x/2/test_asq);
466  }
467  double v_emi=test_x*vw;
468  Eemit[0]=(1.0/sqrt(1.0-v_emi*v_emi/Physics::c/Physics::c)-1)*Physics::m_e*1.0e9;
469  } else {
470  double test_s=1.0;
471  double f_x=0;
472  double test_x=0;
473  while (test_s>f_x) {// rejection & acceptance method
474  test_s=IpplRandom();
475  test_s*=f_max;
476  test_x=IpplRandom();
477  test_x*=10*test_a;// truncation range for emission energy(0,10*test_a);
478  f_x=test_x/test_asq*exp(-test_x/test_a);// thermal distribution for energy
479  }
480 
481  Eemit[0]=test_x;
482 
483 
484  }
485 
486  }
487  // else The incident particle will be marked for deletion
488 
489 
490  } else {
491  seType = 2;
492  if (!nEmissionMode) {
493  if( Options::ppdebug ) {
494  double test_s=1.0;
495  double f_x=0;
496  double test_x=0;
497  while (test_s>f_x) {
498  test_s=IpplRandom();
499  test_s*=f_max;
500  test_x=IpplRandom();
501  test_x*=10*test_a;// truncation range for emission speed(0,10*test_a);
502  f_x=test_x/test_asq*exp(-test_x*test_x/2/test_asq);
503  }
504  double v_emi=test_x*vw;
505  Eemit[0]=(1.0/sqrt(1.0-v_emi*v_emi/Physics::c/Physics::c)-1)*Physics::m_e*1.0e9;
506  } else {
507  double test_s=1.0;
508  double f_x=0;
509  double test_x=0;
510  while (test_s>f_x) {
511  test_s=IpplRandom();
512  test_s*=f_max;
513  test_x=IpplRandom();
514  test_x*=10*test_a;// truncation range for emission energy(0,10*test_a);
515  f_x=test_x/test_asq*exp(-test_x/test_a);// thermal distribution for energy
516  }
517 
518  Eemit[0]=test_x;
519 
520  }
521 
522  } else {
523  if( Options::ppdebug ) {
524  /*=======================Maxwellian Distribution====================*/
525  // the velocity distribution for Vaughan's model:fu=u*vw*vw/vt/vt*exp(-u*u*vw*vw/vt/vt) valid for benchmarking;
526 
527  for(int i = 0; i < seNum; i++) {
528  double test_s=1.0;
529  double f_x=0;
530  double test_x=0;
531  while (test_s>f_x) {
532  test_s=IpplRandom();
533  test_s*=f_max;
534  test_x=IpplRandom();
535  test_x*=10*test_a;//range for normalized emission speed(0,10*test_a);
536  f_x=test_x/test_asq*exp(-test_x*test_x/2/test_asq);
537  }
538  double v_emi=test_x*vw;
539  Eemit[i]=(1.0/sqrt(1.0-v_emi*v_emi/9.0/1e16)-1)*Physics::m_e*1.0e9;
540 
541  }
542 
543  /*---------------------End Of Maxwellian Distribution--------------------*/
544  } else {// energy thermal distribution for Vaughan model
545  for(int i = 0; i < seNum; i++) {
546  double test_s=1.0;
547  double f_x=0;
548  double test_x=0;
549  while (test_s>f_x) {
550  test_s=IpplRandom();
551  test_s*=f_max;
552  test_x=IpplRandom();
553  test_x*=10*test_a;// truncation range for emission energy(0,10*test_a);
554  f_x=test_x/test_asq*exp(-test_x/test_a);// thermal distribution for energy
555  }
556 
557  Eemit[i]=test_x;
558  }
559  }
560  }
561  }
562 
563  Vector_t z_unit = TriNorm;
564  Vector_t x_unit = localX - interCoords_l;
565  double tmp = sqrt(dot(x_unit,x_unit));
566  x_unit /= tmp;
567  Vector_t y_unit = cross(z_unit,x_unit);
568 
569  size_t lowMark = itsBunch->getLocalNum();
570  if (!nEmissionMode) {
571  double gamma_const = Eemit[0] / Physics::m_e/1.0e9 + 1.0;
572  double beta_const = sqrt(1.0 - 1.0 / pow(gamma_const, 2.0));
573  double P_emitted = gamma_const * beta_const;
574 
575 
576  Vector_t P_local;
577  Vector_t P_global = (0.0);
578  if( Options::ppdebug ) {
579 
580  P_global = P_emitted*TriNorm_l;//1D for parallel plate benchmark
581 
582  } else {
583  /*==================3D Vaughan' Model==========================================*/
584  P_local[2] = P_emitted * cos(emiTheta[0]);
585  P_local[1] = P_emitted * sin(emiTheta[0]) * sin(emiPhi[0]);
586  P_local[0] = P_emitted * sin(emiTheta[0]) * cos(emiPhi[0]);
587  P_global = P_local[0]*x_unit + P_local[1]*y_unit + P_local[2]*z_unit;// the same cosine distribution as Pivi's model
588  /*=============================================================================*/
589  }
590 
591  itsBunch->create(1);
592  itsBunch->R[lowMark] = interCoords_l;
593  itsBunch->P[lowMark] = P_global;
594  itsBunch->Bin[lowMark] = 0;
595  itsBunch->PType[lowMark] = ParticleType::NEWSECONDARY;
596  itsBunch->TriID[lowMark] = 0;
597  itsBunch->Q[lowMark] = incQ_l*seyNum;
598  itsBunch->Ef[lowMark] = Vector_t(0.0);
599  itsBunch->Bf[lowMark] = Vector_t(0.0);
600  itsBunch->dt[lowMark] = itsBunch->getdT();
601 
602  } else {
603  for(size_t i = 0; i < (size_t) seNum; i++) {
604 
605  double gamma_const = Eemit[i] / Physics::m_e/1.0e9 + 1.0;
606  double beta_const = sqrt(1.0 - 1.0 / pow(gamma_const, 2.0));
607  double P_emitted = gamma_const * beta_const;
608 
609 
610  Vector_t P_local;
611  Vector_t P_global = (0.0);
612  if( Options::ppdebug ) {
613 
614  P_global = P_emitted*TriNorm_l;//1D for parallel plate benchmark
615 
616  } else {
617  /*==================3D Vaughan' Model==========================================*/
618  P_local[2] = P_emitted * cos(emiTheta[i]);
619  P_local[1] = P_emitted * sin(emiTheta[i]) * sin(emiPhi[i]);
620  P_local[0] = P_emitted * sin(emiTheta[i]) * cos(emiPhi[i]);
621  P_global = P_local[0]*x_unit + P_local[1]*y_unit + P_local[2]*z_unit;//the same cosine distribution as Pivi's model
622  /*=============================================================================*/
623  }
624 
625  itsBunch->create(1);
626  itsBunch->R[lowMark+i] = interCoords_l;
627 
628  itsBunch->P[lowMark+i] = P_global;
629  itsBunch->Bin[lowMark+i] = 0;
630  itsBunch->PType[lowMark+i] = ParticleType::NEWSECONDARY;
631  itsBunch->TriID[lowMark+i] = 0;
632  itsBunch->Q[lowMark+i] = incQ_l;
633  itsBunch->Ef[lowMark+i] = Vector_t(0.0);
634  itsBunch->Bf[lowMark+i] = Vector_t(0.0);
635  itsBunch->dt[lowMark+i] = itsBunch->getdT();
636 
637 
638  }
639  }
640 
641  IpplTimings::stopTimer(TPnSec_m);
642 
643 }
644 
645 void SecondaryEmissionPhysics::calcEmiNum(const double incEnergy,
646  const double cosTheta,
647  int &seNum,
648  const double &vSeyZero,
649  const double &vEzero,
650  const double &vSeyMax,
651  const double &vEmax,
652  const double &vKenergy,
653  const double &vKtheta,
654  double &seyNum) {// For Vaughan's model.
655 
656  double vSEY = 0;
657 
658  if (incEnergy<vEzero) {
659  vSEY = vSeyZero;
660  } else {
661 
662  double theta = acos(cosTheta);
663  double delta_max = vSeyMax*(1.0+vKtheta*theta*theta/Physics::two_pi);//here the symbols k are different with reference.
664  double E_max = vEmax*(1.0+vKenergy*theta*theta/Physics::two_pi);
665  PAssert_GT(E_max - vEzero, 0);
666  double v = (incEnergy-vEzero)/(E_max-vEzero);
667 
668  if (v<=3.6) {
669 
670  if (v<1.0) {
671  vSEY = delta_max*pow(v*exp(1.0-v),0.56);
672  }else {
673  vSEY = delta_max*pow(v*exp(1.0-v),0.25);
674  }
675 
676  }else {
677  vSEY = delta_max*1.125/pow(v,0.35);
678  }
679 
680  }
681  double L = exp(-vSEY);// poisson distribution: Knuth's algorithm.
682  int k = 0;
683  double p = 1.0;
684  do {
685  k++;
686  double u = IpplRandom();
687  p*=u;
688  }
689  while (p>L);
690  seNum = k-1;
691  seyNum = vSEY;
692 
693 }
694 
695 void SecondaryEmissionPhysics::calcEmiNum(const double incEnergy, const double cosTheta, const double *prob, int &seNum) {// For Furman-Pivi's model
696 
697  double prob_max = 0.0;
698  // Acceptance-rejection methods to generate random number with specified distribution.
699 
700  for(int i = 0; i < 11; i++) {
701 
702  if(prob[i] > prob_max) {
703  prob_max = prob[i];
704  }
705 
706  }
707  double pY = 1.0;
708  double pX = 0.0;
709  while(pY > pX) {
710 
711  double rand1 = IpplRandom();
712  //double rand1 = (*rand)(rng);
713  seNum = (int)(rand1 * 11.0);//fix me
714  pX = prob[seNum];
715  double rand2 = IpplRandom();
716  //double rand2 = (*rand)(rng);
717  pY = prob_max * rand2;
718  }
719 
720 }
721 
722 
723 double SecondaryEmissionPhysics::calcDeltats(const double incEnergy, const double cosTheta) {
724 
725  double seypeak = seYPeakTS_m * (1 + seTOneTS_m * (1.0 - pow(cosTheta, seTTwoTS_m))); //formula III.E (48a)
726  double seepeak = seEPeakTS_m * (1 + seTThreeTS_m * (1.0 - pow(cosTheta, seTFourTS_m))); //formula III.E (48b)
727  double tmpx = incEnergy / seepeak;
728  double tmpD = seSTS_m * tmpx / (seSTS_m - 1 + pow(tmpx, seSTS_m)); //formula III.D (32)
729  double ret = seypeak * tmpD; //formula III.D (31)
730  return ret;
731 
732 }
733 
734 double SecondaryEmissionPhysics::calcDeltar(const double incEnergy, const double cosTheta) {
735 
736  double tmp = pow(incEnergy / seERed_m, seR_m);
737  double ret = sePRed_m * (1.0 - exp(-1 * tmp)); //formula III.D (28)
738  ret = ret * (1.0 + seROne_m * (1.0 - pow(cosTheta, seRTwo_m))); //formula III.E (47b)
739  return ret;
740 
741 }
742 
743 
744 double SecondaryEmissionPhysics::calcDeltae(const double incEnergy, const double cosTheta) {
745 
746  double tmp = pow(std::abs(incEnergy - seEScatPeak_m) / seW_m, seP_m) / seP_m;
747  double ret = sePScat_m + (sePScatPeak_m - sePScat_m) * exp(-1 * tmp); //formula III.D (25)
748  ret = ret * (1.0 + seEOne_m * (1.0 - pow(cosTheta, seETwo_m))); //formula III.E (47a)
749  return ret;
750 
751 }
752 
753 
754 double SecondaryEmissionPhysics::calcProb(const double incEnergy, const double cosTheta, double *prob) {
755 
756  deltae_m = calcDeltae(incEnergy, cosTheta);
757  deltar_m = calcDeltar(incEnergy, cosTheta);
758  deltats_m = calcDeltats(incEnergy, cosTheta);
759 
760  double tmp = 1.0 - deltae_m - deltar_m;
761  double p = deltats_m / tmp / 10.0;
762  double q = 1.0 - p;
763  double b[11];
764  b[0] = 1.0;
765  b[1] = 10.0;
766  b[2] = 45.0;
767  b[3] = 120.0;
768  b[4] = 210.0;
769  b[5] = 252.0;
770  b[6] = 210.0;
771  b[7] = 120.0;
772  b[8] = 45.0;
773  b[9] = 10.0;
774  b[10] = 1.0;
775  for(int i = 0; i < 11; i++) {
776  prob[i] = tmp * b[i] * pow(p, i) * pow(q, (10 - i));
777  }
778  prob[1] = prob[1] + deltae_m + deltar_m;
779 
780  /*==============================================*/
781 
782  return (deltae_m+deltar_m+deltats_m);
783  //cout << "sum prob: " << sum << endl;
784  /*==============================================*/
785 }
786 
787 void SecondaryEmissionPhysics::setSeMaterial(int material_num) {
788 
789  if(material_num == 0) {
790  seAlpha_m = 1.0;
791  sePScat_m = 0.02;
792  sePScatPeak_m = 0.496;
793  seEScatPeak_m = 0;
794  seW_m = 60.86;
795  seP_m = 1.0;
796  seDelta_m = 2.0;
797  seEOne_m = 0.26;
798  seETwo_m = 2;
799 
800  sePRed_m = 0.2;
801  seERed_m = 0.041;
802  seR_m = 0.104;
803  seQ_m = 0.5;
804  seROne_m = 0.26;
805  seRTwo_m = 2;
806 
807  seYPeakTS_m = 1.8848;
808  seEPeakTS_m = 276.8;
809  seSTS_m = 1.54;
810  seTOneTS_m = 0.66;
811  seTTwoTS_m = 0.8;
812  seTThreeTS_m = 0.7;
813  seTFourTS_m = 1.0;
814  seEPeakTot_m = 271;
815  seYPeakTot_m = 2.1;
816  sePn_m[0] = 2.5;
817  sePn_m[1] = 3.3;
818  sePn_m[2] = 2.5;
819  sePn_m[3] = 2.5;
820  sePn_m[4] = 2.8;
821  sePn_m[5] = 1.3;
822  sePn_m[6] = 1.5;
823  sePn_m[7] = 1.5;
824  sePn_m[8] = 1.5;
825  sePn_m[9] = 1.5;
826  seEpsn_m[0] = 1.5;
827  seEpsn_m[1] = 1.75;
828  seEpsn_m[2] = 1.0;
829  seEpsn_m[3] = 3.75;
830  seEpsn_m[4] = 8.5;
831  seEpsn_m[5] = 11.5;
832  seEpsn_m[6] = 2.5;
833  seEpsn_m[7] = 3.0;
834  seEpsn_m[8] = 2.5;
835  seEpsn_m[9] = 3.0;
836  }
837  if(material_num == 1) {
838  seAlpha_m = 1.0;
839  sePScat_m = 0.07;
840  sePScatPeak_m = 0.5;
841  seEScatPeak_m = 0;
842  seW_m = 100;
843  seP_m = 0.9;
844  seDelta_m = 1.9;
845  seEOne_m = 0.26;
846  seETwo_m = 2;
847 
848  sePRed_m = 0.74;
849  seERed_m = 40;
850  seR_m = 1;
851  seQ_m = 0.4;
852  seROne_m = 0.26;
853  seRTwo_m = 2;
854 
855  seYPeakTS_m = 1.22;
856  seEPeakTS_m = 310;
857  seSTS_m = 1.813;
858  seTOneTS_m = 0.66;
859  seTTwoTS_m = 0.8;
860  seTThreeTS_m = 0.7;
861  seTFourTS_m = 1.0;
862  seEPeakTot_m = 292;
863  seYPeakTot_m = 2.05;
864 
865  sePn_m[0] = 1.6;
866  sePn_m[1] = 2.0;
867  sePn_m[2] = 1.8;
868  sePn_m[3] = 4.7;
869  sePn_m[4] = 1.8;
870  sePn_m[5] = 2.4;
871  sePn_m[6] = 1.8;
872  sePn_m[7] = 1.8;
873  sePn_m[8] = 2.3;
874  sePn_m[9] = 1.8;
875  seEpsn_m[0] = 3.9;
876  seEpsn_m[1] = 6.2;
877  seEpsn_m[2] = 13.0;
878  seEpsn_m[3] = 8.8;
879  seEpsn_m[4] = 6.25;
880  seEpsn_m[5] = 2.25;
881  seEpsn_m[6] = 9.2;
882  seEpsn_m[7] = 5.3;
883  seEpsn_m[8] = 17.8;
884  seEpsn_m[9] = 10;
885  }
886 }
887 
888 
889 
890 /*==========================================
891  *http://www.taygeta.com/random/gaussian.html
892  *return a gaussian distributed random number
893  *==========================================*/
894 
895 double SecondaryEmissionPhysics::gaussRand() {
896 
897  double x1;
898  double x2;
899  double w;
900  do {
901  x1 = 2.0 * IpplRandom() - 1.0;
902  //x1 = 2.0 * (*rand)(rng) - 1.0;
903  x2 = 2.0 * IpplRandom() - 1.0;
904  //x2 = 2.0 * (*rand)(rng) - 1.0;
905  w = x1 * x1 + x2 * x2;
906  } while(w >= 1.0);
907 
908  w = sqrt((-2.0 * log(w)) / w);
909  double ret = x1 * w;
910  return ret;
911 
912 }
913 
914 double SecondaryEmissionPhysics::gammp(const double a, const double x) {
915  // Returns the incomplete gamma function P .a; x/.
916  static const int ASWITCH = 100; //When to switch to quadrature method.
917  if(x < 0.0 || a <= 0.0)
918  throw("bad args in gammp");
919  if(x == 0.0)
920  return 0.0;
921  else if((int)a >= ASWITCH)
922  return gammpapprox(a, x, 1); //Quadrature.
923  else if(x < a + 1.0)
924  return gser(a, x); //Use the series representation.
925  else
926  return 1.0 - gcf(a, x); //Use the continued fraction representation.
927 }
928 
929 double SecondaryEmissionPhysics::gser(const double a, const double x) {
930  // Returns the incomplete gamma function P .a; x/ evaluated by its series representation.
931  // Also sets ln .a/ as gln. User should not call directly.
932  double sum, del, ap;
933  double gln = gammln(a);
934  ap = a;
935  del = sum = 1.0 / a;
936  for(;;) {
937  ++ap;
938  del *= x / ap;
939  sum += del;
940  if(std::abs(del) < std::abs(sum)*myeps::EPS) {
941  return sum * exp(-x + a * log(x) - gln);
942  }
943  }
944 }
945 
946 double SecondaryEmissionPhysics::gcf(const double a, const double x) {
947  // Returns the incomplete gamma function Q.a; x/ evaluated by its continued fraction representation. Also sets ln .a/ as gln. User should not call directly.
948  int i;
949  double an, b, c, d, del, h;
950  double gln = gammln(a);
951  b = x + 1.0 - a; // Set up for evaluating continued fraction
952  c = 1.0 / myeps::FPMIN; // by modified Lentz‘¡¯s method (5.2)
953  d = 1.0 / b; // with b0 D 0.
954  h = d;
955  for(i = 1;; i++) {
956  //Iterate to convergence.
957  an = -i * (i - a);
958  b += 2.0;
959  d = an * d + b;
960  if(std::abs(d) < myeps::FPMIN) d = myeps::FPMIN;
961  c = b + an / c;
962  if(std::abs(c) < myeps::FPMIN) c = myeps::FPMIN;
963  d = 1.0 / d;
964  del = d * c;
965  h *= del;
966  if(std::abs(del - 1.0) <= myeps::EPS) break;
967  }
968  return exp(-x + a * log(x) - gln) * h; //Put factors in front.
969 }
970 
971 double SecondaryEmissionPhysics::gammpapprox(double a, double x, int psig) {
972  // Incomplete gamma by quadrature. Returns P .a; x/ or Q.a; x/, when psig is 1 or 0,respectively. User should not call directly.
973 
974  const double y[18] = {0.0021695375159141994,
975  0.011413521097787704, 0.027972308950302116, 0.051727015600492421,
976  0.082502225484340941, 0.12007019910960293, 0.16415283300752470,
977  0.21442376986779355, 0.27051082840644336, 0.33199876341447887,
978  0.39843234186401943, 0.46931971407375483, 0.54413605556657973,
979  0.62232745288031077, 0.70331500465597174, 0.78649910768313447,
980  0.87126389619061517, 0.95698180152629142
981  };
982  const double w[18] = {0.0055657196642445571,
983  0.012915947284065419, 0.020181515297735382, 0.027298621498568734,
984  0.034213810770299537, 0.040875750923643261, 0.047235083490265582,
985  0.053244713977759692, 0.058860144245324798, 0.064039797355015485,
986  0.068745323835736408, 0.072941885005653087, 0.076598410645870640,
987  0.079687828912071670, 0.082187266704339706, 0.084078218979661945,
988  0.085346685739338721, 0.085983275670394821
989  };
990  int j;
991  double xu, t, sum, ans;
992  double a1 = a - 1.0, lna1 = log(a1), sqrta1 = sqrt(a1);
993  int ngau = 18;
994  double gln = gammln(a);// Set how far to integrate into the tail:
995  if(x > a1) xu = max_Gamma(a1 + 11.5 * sqrta1, x + 6.0 * sqrta1);
996  else xu = max_Gamma(0., min_Gamma(a1 - 7.5 * sqrta1, x - 5.0 * sqrta1));
997  sum = 0;
998  for(j = 0; j < ngau; j++) {
999  //Gauss-Legendre.
1000  t = x + (xu - x) * y[j];
1001  sum += w[j] * exp(-(t - a1) + a1 * (log(t) - lna1));
1002  }
1003  ans = sum * (xu - x) * exp(a1 * (lna1 - 1.) - gln);
1004  return (psig ? (ans > 0.0 ? 1.0 - ans : -ans) : (ans >= 0.0 ? ans : 1.0 + ans));
1005 }
1006 
1007 double SecondaryEmissionPhysics::gammln(const double xx) {
1008  // Returns the value ln%GÅ’%@.xx/ for xx > 0.
1009  int j;
1010  double x, tmp, y, ser;
1011  static const double cof[14] = {57.1562356658629235, -59.5979603554754912,
1012  14.1360979747417471, -0.491913816097620199, .339946499848118887e-4,
1013  .465236289270485756e-4, -.983744753048795646e-4, .158088703224912494e-3,
1014  -.210264441724104883e-3, .217439618115212643e-3, -.164318106536763890e-3,
1015  .844182239838527433e-4, -.261908384015814087e-4, .368991826595316234e-5
1016  };
1017  if(xx <= 0) throw("bad arg in gammln");
1018  y = x = xx;
1019  tmp = x + 5.24218750000000000; // Rational 671/128.
1020  tmp = (x + 0.5) * log(tmp) - tmp;
1021  ser = 0.999999999999997092;
1022  for(j = 0; j < 14; j++) ser += cof[j] / ++y;
1023  return tmp + log(2.5066282746310005 * ser / x);
1024 }
1025 
1026 double SecondaryEmissionPhysics::invgammp(double p, double a) {
1027  //Returns x such that P .a; x/ D p for an argument p between 0 and 1.
1028  int j;
1029  double x, err, t, u, pp, lna1 = 0.0, afac = 0.0, a1 = a - 1;
1030  const double EPS = 1.e-8; //Accuracy is the square of EPS.
1031  double gln = gammln(a);
1032  if(a <= 0.) throw("a must be pos in invgammap");
1033  if(p >= 1.) return max_Gamma(100., a + 100.*sqrt(a));
1034  if(p <= 0.) return 0.0;
1035  if(a > 1.) {
1036  //Initial guess based on reference [1].
1037  lna1 = log(a1);
1038  afac = exp(a1 * (lna1 - 1.) - gln);
1039  pp = (p < 0.5) ? p : 1. - p;
1040  t = sqrt(-2.*log(pp));
1041  x = (2.30753 + t * 0.27061) / (1. + t * (0.99229 + t * 0.04481)) - t;
1042  if(p < 0.5) x = -x;
1043  x = max_Gamma(1.e-3, a * pow(1. - 1. / (9.*a) - x / (3.*sqrt(a)), 3));
1044  } else {
1045  //Initial guess based on equations (6.2.8) and (6.2.9).
1046  t = 1.0 - a * (0.253 + a * 0.12);
1047 
1048  if(p < t) x = pow(p / t, 1. / a);
1049  else x = 1. - log(1. - (p - t) / (1. - t));
1050  }
1051  for(j = 0; j < 12; j++) {
1052  if(x <= 0.0) return 0.0;
1053  //x too small to compute accurately.
1054  err = gammp(a, x) - p;
1055  if(a > 1.) t = afac * exp(-(x - a1) + a1 * (log(x) - lna1));
1056  else t = exp(-x + a1 * log(x) - gln);
1057  u = err / t;
1058  x -= (t = u / (1. - 0.5 * min_Gamma(1., u * ((a - 1.) / x - 1))));
1059  //Halley‘¡¯s method.
1060  if(x <= 0.) x = 0.5 * (x + t);
1061  //Halve old value if x tries to go negative.
1062  if(std::abs(t) < EPS * x) break;
1063  }
1064  return x;
1065 }
1066 
1067 double SecondaryEmissionPhysics::betai(const double x, const double a, const double b) {
1068  //cout<<"betai called"<<endl;
1069  //Returns incomplete beta function Ix.a; b/ for positive a and b, and x between 0 and 1.
1070  static const int SWITCH = 3000;
1071  double bt;
1072  /*=========================debug code=====================
1073  if (a==b) {
1074  cout<< " x in betai "<<x<<endl;
1075  }
1076  =========================debug code=====================*/
1077  if(a <= 0.0 || b <= 0.0) {
1078  //cout<<"betai 1"<<endl;
1079  throw("Bad a or b in routine betai");
1080  }
1081 
1082  if(x < 0.0 || x > 1.0) {
1083  //cout<<"betai 2"<<endl;
1084  throw("Bad x in routine betai");
1085 
1086  }
1087  if(x == 0.0 || x == 1.0) {
1088  //cout<<"betai 3"<<endl;
1089  return x;
1090  }
1091  if(a > SWITCH && b > SWITCH) {
1092  //cout<<"betai 4"<<endl;
1093  return betaiapprox(a, b, x);
1094  }
1095  bt = exp(gammln(a + b) - gammln(a) - gammln(b) + a * log(x) + b * log(1.0 - x));
1096  if(x < (a + 1.0) / (a + b + 2.0)) {
1097  //cout<<"betai 5"<<endl;
1098  return bt * betacf(a, b, x) / a;
1099  } else {
1100  //cout<<"betai 6"<<endl;
1101  /*=========================debug code=====================
1102  if (a==b) {
1103  cout<<" betacf(b, a, 1.0 - x) = "<< betacf(b, a, 1.0 - x)<<" a = "<<a<<" b = "<<b<<endl;
1104  }
1105  =========================================================*/
1106  return 1.0 - bt * betacf(b, a, 1.0 - x) / b;
1107  }
1108 }
1109 
1110 double SecondaryEmissionPhysics::betacf(const double a, const double b, const double x) {
1111  //Evaluates continued fraction for incomplete beta function by modified Lentzâs method (5.2). User should not call directly.
1112  int m, m2;
1113  double aa, c, d, del, h, qab, qam, qap;
1114  qab = a + b; //These qab will be used in factors that occur in the coefficients (6.4.6).
1115  qap = a + 1.0;
1116  qam = a - 1.0;
1117  c = 1.0; //First step of Lentzâs method.
1118 
1119  d = 1.0 - qab * x / qap;
1120  /*========debug====================
1121  if (a == b) {
1122  cout<<"d = "<<d<<" qap = "<<qap<<" qab = "<<qab<<" x = "<<x<<endl;
1123  }
1124  ========debug====================*/
1125  if(std::abs(d) < myeps::FPMIN) d = myeps::FPMIN;
1126  d = 1.0 / d;
1127  h = d;
1128  // cout<<"h = "<<h<<" myeps::FPMIN = "<<myeps::FPMIN<<endl;
1129  for(m = 1; m < 10000; m++) {
1130  m2 = 2 * m;
1131  aa = m * (b - m) * x / ((qam + m2) * (a + m2));
1132  d = 1.0 + aa * d; //One step (the even one) of the recurrence.
1133  if(std::abs(d) < myeps::FPMIN)
1134  d = myeps::FPMIN;
1135  c = 1.0 + aa / c;
1136  if(std::abs(c) < myeps::FPMIN)
1137  c = myeps::FPMIN;
1138  d = 1.0 / d;
1139  h *= d * c;
1140  aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
1141  d = 1.0 + aa * d; //Next step of the recurrence (the odd one).
1142  if(std::abs(d) < myeps::FPMIN)
1143  d = myeps::FPMIN;
1144  c = 1.0 + aa / c;
1145  if(std::abs(c) < myeps::FPMIN)
1146  c = myeps::FPMIN;
1147  d = 1.0 / d;
1148  del = d * c;
1149  h *= del;
1150  if(std::abs(del - 1.0) <= myeps::EPS)
1151  break;
1152  }
1153  return h;
1154 }
1155 
1156 double SecondaryEmissionPhysics::betaiapprox(double a, double b, double x) {
1157  //Incomplete beta by quadrature. Returns Ix.a; b/. User should not call directly.
1158  const double y[18] = {0.0021695375159141994,
1159  0.011413521097787704, 0.027972308950302116, 0.051727015600492421,
1160  0.082502225484340941, 0.12007019910960293, 0.16415283300752470,
1161  0.21442376986779355, 0.27051082840644336, 0.33199876341447887,
1162  0.39843234186401943, 0.46931971407375483, 0.54413605556657973,
1163  0.62232745288031077, 0.70331500465597174, 0.78649910768313447,
1164  0.87126389619061517, 0.95698180152629142
1165  };
1166  const double w[18] = {0.0055657196642445571,
1167  0.012915947284065419, 0.020181515297735382, 0.027298621498568734,
1168  0.034213810770299537, 0.040875750923643261, 0.047235083490265582,
1169  0.053244713977759692, 0.058860144245324798, 0.064039797355015485,
1170  0.068745323835736408, 0.072941885005653087, 0.076598410645870640,
1171  0.079687828912071670, 0.082187266704339706, 0.084078218979661945,
1172  0.085346685739338721, 0.085983275670394821
1173  };
1174  int j;
1175  double xu, t, sum, ans;
1176  double a1 = a - 1.0, b1 = b - 1.0, mu = a / (a + b);
1177  double lnmu = log(mu), lnmuc = log(1. - mu);
1178  t = sqrt(a * b / (sqrt(a + b) * (a + b + 1.0)));
1179  if(x > a / (a + b)) { //Set how far to integrate into the tail:
1180  if(x >= 1.0)
1181  return 1.0;
1182  xu = min_Gamma(1., max_Gamma(mu + 10.*t, x + 5.0 * t));
1183  } else {
1184 
1185  if(x <= 0.0)
1186  return 0.0;
1187  xu = max_Gamma(0., min_Gamma(mu - 10.*t, x - 5.0 * t));
1188  }
1189  sum = 0;
1190  for(j = 0; j < 18; j++) { //Gauss-Legendre.
1191  t = x + (xu - x) * y[j];
1192  sum += w[j] * exp(a1 * (log(t) - lnmu) + b1 * (log(1 - t) - lnmuc));
1193  }
1194  ans = sum * (xu - x) * exp(a1 * lnmu - gammln(a) + b1 * lnmuc - gammln(b) + gammln(a + b));
1195  return ans > 0.0 ? 1.0 - ans : -ans;
1196 }
1197 
1198 double SecondaryEmissionPhysics::invbetai(double p, double a, double b) {
1199  // Inverse of incomplete beta function. Returns x such that Ix.a; b/ D p for argument p between 0 and 1.
1200  const double EPS = 1.e-8;
1201  double pp, t, u, err, x, al, h, w, afac, a1 = a - 1., b1 = b - 1.;
1202  int j;
1203 
1204  if(p <= 0.)
1205  return 0.;
1206  else if(p >= 1.)
1207  return 1.;
1208  else if(a >= 1. && b >= 1.) { // Set initial guess. See text.
1209  pp = (p < 0.5) ? p : 1. - p;
1210  t = sqrt(-2.*log(pp));
1211  x = (2.30753 + t * 0.27061) / (1. + t * (0.99229 + t * 0.04481)) - t;
1212  /*========debug====================
1213  if (a==b) {
1214 
1215  cout<<"x = "<<x<<" p = "<<p<<" w = "<<w<<" t = "<<t<<" al = "<<al<<" h = "<<h<<endl;
1216 
1217  }
1218  ========debug====================*/
1219  // if(p < 0.5)//origin code from numerical ricipes.
1220  if(x < 0.0)//fixed bug.
1221  x = -x;
1222  al = (sqrt(x) - 3.) / 6.;
1223  h = 2. / (1. / (2.*a - 1.) + 1. / (2.*b - 1.));
1224  w = (x * sqrt(al + h) / h) - (1. / (2.*b - 1) - 1. / (2.*a - 1.)) * (al + 5. / 6. - 2. / (3.*h));
1225  x = a / (a + b * exp(2.*w));
1226 
1227  } else {
1228  double lna = log(a / (a + b)), lnb = log(b / (a + b));
1229  t = exp(a * lna) / a;
1230  u = exp(b * lnb) / b;
1231  w = t + u;
1232  if(p < t / w)
1233  x = pow(a * w * p, 1. / a);
1234  else
1235  x = 1. - pow(b * w * (1. - p), 1. / b);
1236  }
1237  afac = -gammln(a) - gammln(b) + gammln(a + b);
1238  for(j = 0; j < 10; j++) {
1239  if(x == 0. || x == 1.) // a or b too small for accurate calculation.
1240  return x;
1241  err = betai(x, a, b) - p;
1242  /*=========================debug code=====================
1243  if (a==b) {
1244  cout<<"p = "<<p<<" err = "<<err<<" t = "<<t<<" u = "<<u<<" x = "<<x<<endl;
1245  }
1246  =========================================================*/
1247  t = exp(a1 * log(x) + b1 * log(1. - x) + afac);
1248  u = err / t; //Halley:
1249  x -= (t = u / (1. - 0.5 * min_Gamma(1., u * (a1 / x - b1 / (1. - x)))));
1250  /*=========================debug code=====================
1251  if (a==b) {
1252  cout<<"err = "<<err<<" t = "<<t<<" u = "<<u<<" x = "<<x<<endl;
1253  }
1254  =========================================================*/
1255  if(x <= 0.)
1256  x = 0.5 * (x + t); // Bisect if x tries to go neg or > 1.
1257  if(x >= 1.)
1258  x = 0.5 * (x + t + 1.);
1259  if(std::abs(t) < EPS * x && j > 0)
1260  break;
1261  }
1262  return x;
1263 }
SecondaryEmissionPhysics::invbetai
double invbetai(double p, double a, double b)
Definition: SecondaryEmissionPhysics.cpp:1198
PartBunchBase::P
ParticleAttrib< Vector_t > P
Definition: PartBunchBase.h:484
abs
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
Definition: IpplExpressions.h:62
acos
PETE_TUTree< FnArcCos, typename T::PETE_Expr_t > acos(const PETE_Expr< T > &l)
Definition: PETE.h:808
SecondaryEmissionPhysics::calcProb
double calcProb(const double incEnergy, const double cosTheta, double *prob)
Definition: SecondaryEmissionPhysics.cpp:754
Physics::e
constexpr double e
The value of .
Definition: Physics.h:40
PartBunchBase::Ef
ParticleAttrib< Vector_t > Ef
Definition: PartBunchBase.h:488
PartBunchBase::create
void create(size_t M)
Definition: PartBunchBase.hpp:2319
SecondaryEmissionPhysics::gammpapprox
double gammpapprox(double a, double x, int psig)
Definition: SecondaryEmissionPhysics.cpp:971
fabs
PETE_TUTree< FnFabs, typename T::PETE_Expr_t > fabs(const PETE_Expr< T > &l)
Definition: PETE.h:815
sin
Tps< T > sin(const Tps< T > &x)
Sine.
Definition: TpsMath.h:111
PartBunchBase::Q
ParticleAttrib< double > Q
Definition: PartBunchBase.h:485
Physics::two_pi
constexpr double two_pi
The value of .
Definition: Physics.h:34
exp
Tps< T > exp(const Tps< T > &x)
Exponential.
Definition: TpsMath.h:165
PartBunchBase::getdT
double getdT() const
Definition: PartBunchBase.hpp:1076
IpplRandom
RandomNumberGen IpplRandom
Definition: RandomNumberGen.cpp:30
sum
T::PETE_Expr_t::PETE_Return_t sum(const PETE_Expr< T > &expr)
Definition: PETE.h:1213
log
Tps< T > log(const Tps< T > &x)
Natural logarithm.
Definition: TpsMath.h:182
dot
double dot(const Vector3D &lhs, const Vector3D &rhs)
Vector dot product.
Definition: Vector3D.cpp:118
Physics::m_e
constexpr double m_e
The electron rest mass in GeV.
Definition: Physics.h:85
PAssert_LT
#define PAssert_LT(a, b)
Definition: PAssert.h:121
PAssert_GE
#define PAssert_GE(a, b)
Definition: PAssert.h:124
PartBunchBase::PType
ParticleAttrib< short > PType
Definition: PartBunchBase.h:495
SecondaryEmissionPhysics::gammp
double gammp(const double a, const double x)
Definition: SecondaryEmissionPhysics.cpp:914
SecondaryEmissionPhysics::calcDeltae
double calcDeltae(const double incEnergy, const double cosTheta)
Definition: SecondaryEmissionPhysics.cpp:744
PartBunchBase::TriID
ParticleAttrib< int > TriID
Definition: PartBunchBase.h:496
Physics::c
constexpr double c
The velocity of light in m/s.
Definition: Physics.h:52
cross
Vector3D cross(const Vector3D &lhs, const Vector3D &rhs)
Vector cross product.
Definition: Vector3D.cpp:111
IpplTimings::startTimer
static void startTimer(TimerRef t)
Definition: IpplTimings.h:187
pow
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
Definition: TpsMath.h:76
Vector_t
Vektor< double, 3 > Vector_t
Definition: Vektor.h:6
SecondaryEmissionPhysics::calcDeltar
double calcDeltar(const double incEnergy, const double cosTheta)
Definition: SecondaryEmissionPhysics.cpp:734
SecondaryEmissionPhysics::betacf
double betacf(const double a, const double b, const double x)
Definition: SecondaryEmissionPhysics.cpp:1110
SecondaryEmissionPhysics::gcf
double gcf(const double a, const double x)
Definition: SecondaryEmissionPhysics.cpp:946
PartBunchBase::getLocalNum
size_t getLocalNum() const
Definition: PartBunchBase.hpp:2230
SecondaryEmissionPhysics::calcEmiNum
void calcEmiNum(const double incEnergy, const double cosTheta, int &seNum, const double &vSeyZero, const double &vEzero, const double &vSeyMax, const double &vEmax, const double &vKenergy, const double &vKtheta, double &seyNum)
Definition: SecondaryEmissionPhysics.cpp:645
sqrt
Tps< T > sqrt(const Tps< T > &x)
Square root.
Definition: TpsMath.h:91
PartBunchBase::dt
ParticleAttrib< double > dt
Definition: PartBunchBase.h:493
PAssert_GT
#define PAssert_GT(a, b)
Definition: PAssert.h:123
SecondaryEmissionPhysics::betaiapprox
double betaiapprox(double a, double b, double x)
Definition: SecondaryEmissionPhysics.cpp:1156
SecondaryEmissionPhysics::betai
double betai(const double x, const double a, const double b)
Definition: SecondaryEmissionPhysics.cpp:1067
SecondaryEmissionPhysics::gser
double gser(const double a, const double x)
Definition: SecondaryEmissionPhysics.cpp:929
cos
Tps< T > cos(const Tps< T > &x)
Cosine.
Definition: TpsMath.h:129
SecondaryEmissionPhysics::setSeMaterial
void setSeMaterial(int material_num)
Definition: SecondaryEmissionPhysics.cpp:787
PartBunchBase::Bin
ParticleAttrib< int > Bin
Definition: PartBunchBase.h:492
PartBunchBase::Bf
ParticleAttrib< Vector_t > Bf
Definition: PartBunchBase.h:491
PartBunchBase::R
ParticlePos_t & R
Definition: PartBunchBase.h:479
SecondaryEmissionPhysics::calcDeltats
double calcDeltats(const double incEnergy, const double cosTheta)
Definition: SecondaryEmissionPhysics.cpp:723
SecondaryEmissionPhysics::invgammp
double invgammp(double p, double a)
Definition: SecondaryEmissionPhysics.cpp:1026
IpplTimings::getTimer
static TimerRef getTimer(const char *nm)
Definition: IpplTimings.h:182
myeps::EPS
const double EPS
Definition: SecondaryEmissionPhysics.h:15
IpplTimings::stopTimer
static void stopTimer(TimerRef t)
Definition: IpplTimings.h:192
Options::ppdebug
bool ppdebug
ppdebug flag.
Definition: Options.cpp:13
myeps::FPMIN
const double FPMIN
Definition: SecondaryEmissionPhysics.h:16
SecondaryEmissionPhysics::nSec
void nSec(const double &incEnergy, const double &cosTheta, const int &matNumber, int &seNum, int &seType, const double &incQ, const Vector_t &TriNorm, const Vector_t &inteCoords, const Vector_t &localX, PartBunchBase< double, 3 > *itsBunch, double &seyNum, const double &ppVw, const double &vVThermal, const bool nEmissionMode)
Definition: SecondaryEmissionPhysics.cpp:24
Generated on Sat Jan 4 2020 14:48:06 for OPAL (Object Oriented Parallel Accelerator Library) by   doxygen 1.8.5