cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
hydrocollid.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*HCSAR_interp interpolate on collision strengths */
4 /*C6cs123 line collision rates for lower levels of hydrogenic carbon, n=1,2,3 */
5 /*Ca20cs123 line collision rates for lower levels of hydrogenic calcium, n=1,2,3 */
6 /*Hydcs123 Hydrogenic de-excitation collision rates n=1,2,3 */
7 /*H1cs123 hydrogen collision data levels involving 1s,2s,2p,3. */
8 /*Ne10cs123 line collision rates for lower levels of hydrogenic neon, n=1,2,3 */
9 /*He2cs123 line collision strengths for lower levels of helium ion, n=1,2,3, by K Korista */
10 /*Fe26cs123 line collision rates for lower levels of hydrogenic iron, n=1,2,3 */
11 #include "cddefines.h"
12 #include "atmdat.h"
13 #include "atmdat_adfa.h"
14 #include "helike_cs.h"
15 #include "hydrogenic.h"
16 #include "hydro_vs_rates.h"
17 #include "iso.h"
18 #include "opacity.h"
19 #include "phycon.h"
20 #include "thirdparty.h"
21 #include "integrate.h"
22 #include "freebound.h"
23 #include "hydroeinsta.h"
24 
25 STATIC double Fe26cs123(long int i, long int j);
26 STATIC double He2cs123(long int i, long int j);
27 STATIC double Hydcs123(long int ilow, long int ihigh, long int iz, long int chType);
28 STATIC double C6cs123(long int i, long int j);
29 STATIC double Ca20cs123(long int i, long int j);
30 STATIC double Ne10cs123(long int i, long int j);
31 
32 STATIC realnum HCSAR_interp( int ipLo , int ipHi );
33 STATIC realnum HlikeCSInterp( long nelem, long Collider, long nHi, long lHi, long sHi, long nLo, long lLo, long sLo );
34 STATIC double CS_ThermAve_PR78(long ipISO, long nelem, long nHi, long nLo, double deltaE, double temp );
35 STATIC double Therm_ave_coll_str_int_PR78( double EOverKT );
36 inline double C2_PR78(double x, double y);
37 STATIC double CS_PercivalRichards78( double Ebar );
38 
40 static double kTRyd, global_deltaE;
41 
42 static const realnum HCSTE[NHCSTE] = {5802.f,11604.f,34812.f,58020.f,116040.f,174060.f,232080.f,290100.f};
43 
44 /*HCSAR_interp interpolate on collision strengths */
45 STATIC realnum HCSAR_interp( int ipLo , int ipHi )
46 {
47 
48  static int ip=1;
49  realnum cs;
50 
51  DEBUG_ENTRY( "HCSAR_interp()" );
52 
53  if( ipLo==1 && ipHi==2 )
54  {
55  fprintf(ioQQQ,"HCSAR_interp was called for the 2s-2p transition, which it cannot do\n");
57  }
58  if( phycon.te <= HCSTE[0] )
59  {
60  cs = t_ADfA::Inst().h_coll_str( ipLo , ipHi , 0 );
61  }
62  else if( phycon.te >= HCSTE[NHCSTE-1] )
63  {
64  cs = t_ADfA::Inst().h_coll_str( ipLo , ipHi , NHCSTE-1 );
65  }
66  else
67  {
68  /* the ip index is most likely correct since it points to the last temperature */
69  if( (HCSTE[ip-1] >= phycon.te ) || ( phycon.te > HCSTE[ip]) )
70  {
71  /* we must find the temperature in the array */
72  for( ip=1; ip<NHCSTE; ++ip )
73  {
74  if( (HCSTE[ip-1] < phycon.te ) && ( phycon.te <= HCSTE[ip]) )
75  break;
76  }
77  }
78  /* we now have the index */
79  cs = t_ADfA::Inst().h_coll_str( ipLo , ipHi , ip-1 ) +
80  (t_ADfA::Inst().h_coll_str( ipLo , ipHi , ip ) - t_ADfA::Inst().h_coll_str( ipLo , ipHi , ip-1 ) ) / (HCSTE[ip]-HCSTE[ip-1] ) *
81  ((realnum)phycon.te - HCSTE[ip-1] );
82  if( cs <= 0.)
83  {
84  fprintf(ioQQQ," insane cs returned by HCSAR_interp, values are\n");
85  fprintf(ioQQQ,"%.3f %.3f \n", t_ADfA::Inst().h_coll_str( ipLo , ipHi , ip-1 ),t_ADfA::Inst().h_coll_str( ipLo , ipHi , ip ) );
86  }
87  }
88  return(cs);
89 }
90 
91 /*Hydcs123 Hydrogenic de-excitation collision strengths between levels n=1,2,3,
92  * for any charge. routine only called by HydroCSInterp to fill in hydroline arrays
93  * with collision strengths */
95  /* lower principal quantum number, 1, 2, or 3, in this routine
96  * 1 is 1s, 2 is 2s, 3 is 2p, and 4 is 3s, 5 is 3p, and 6 is 3d */
97  long int ipLow,
98  /* upper principal quantum nubmer, 2, 3, or 4 */
99  long int ipHi,
100  /* charge, 0 for hydrogen, 1 for helium, etc */
101  long int nelem,
102  /* = 'e' for electron collisions, ='p' for proton */
103  long int chType)
104 {
105  long int i,
106  j,
107  k;
108  double C,
109  D,
110  EE,
111  expq ,
112  Hydcs123_v,
113  logx,
114  Ratehigh,
115  Ratelow,
116  TeUse,
117  gLo,
118  gHi,
119  q,
120  rate,
121  slope,
122  temp,
123  temphigh,
124  templow,
125  tev,
126  x,
127  QuanNLo,
128  QuanNUp,
129  Charge,
130  ChargeSquared,
131  zhigh,
132  zlow;
133  static const double ap[5] = {-2113.113,729.0084,1055.397,854.632,938.9912};
134  static const double bp[5] = {-6783.515,-377.7190,724.1936,493.1107,735.7466};
135  static const double cp[5] = {-3049.719,226.2320,637.8630,388.5465,554.6369};
136  static const double dp[5] = {3514.5153,88.60169,-470.4055,-329.4914,-450.8459};
137  static const double ep[5] = {0.005251557,0.009059154,0.008725781,0.009952418,0.01098687};
138  static const double ae[5] = {-767.5859,-643.1189,-461.6836,-429.0543,-406.5285};
139  static const double be[5] = {-1731.9178,-1442.548,-1055.364,-980.3079,-930.9266};
140  static const double ce[5] = {-939.1834,-789.9569,-569.1451,-530.1974,-502.0939};
141  static const double de[5] = {927.4773,773.2008,564.3272,524.2944,497.7763};
142  static const double ee[5] = {-0.002528027,-0.003793665,-0.002122103,-0.002234207,-0.002317720};
143  static const double A[2] = {4.4394,0.0};
144  static const double B[2] = {0.8949,0.8879};
145  static const double C0[2] = {-0.6012,-0.2474};
146  static const double C1[2] = {-3.9710,-3.7562};
147  static const double C2[2] = {-4.2176,2.0491};
148  static const double D0[2] = {2.930,0.0539};
149  static const double D1[2] = {1.7990,3.4009};
150  static const double D2[2] = {4.9347,-1.7770};
151 
152  DEBUG_ENTRY( "Hydcs123()" );
153  /* Hydrogenic de-excitation collision rates n=1,2,3
154  * >>refer h1 cs Callaway, J. 1983, Phys Let A, 96, 83
155  * >>refer h1 cs Zygelman, B., & Dalgarno, A. 1987, Phys Rev A, 35, 4085
156  * for 2p-2s only.
157  * The fit from Callaway is in nuclear charge for 1s - 2s,2p only.
158  * For transtions involving level 3, interpolation in Z involving
159  * the functions He2cs123,C6cs123,Ne10cs123,Ca20cs123, Fe26cs123.
160  *
161  * The fits from ZD are for 2p-2s for Z=2,6,12,16,18 only other charges are
162  * interpolated, both electron and proton rates are included,
163  * the variable chType is either 'e' or 'p'.
164  *
165  * ipLow is the lower level and runs from 1 to 3 (1s, 2s, 2p)
166  * ipHi is the upper level and runs from 2 to 6 (2s, 2p, 3s, 3p, 3d) */
167 
168  /* for Callaway fit: */
169  /* for Zygelman and Dalgarno: */
170 
171  /* first entry is 2p, then 2s */
172 
173  /* fit in nuclear charge Z for 2p-2s collisions in hydrogenic species
174  * equation is a+bx+cx^2ln(x)+dexp(x)+eln(x)/x^2, where x=te/Z^2 in a.u.
175  * first are the proton rates: */
176  /* these are electron rates: */
177 
178  /* following is for charged species */
179  /* charge is on scale with He+=1, Li++=2, etc */
180  ASSERT( nelem > ipHYDROGEN );
181  ASSERT( nelem < LIMELM );
182 
183  /* these are the pointers to upper and lower levels. 1=1s, 2=2s, 3=2p, 4=3 */
184  ASSERT( ipLow > 0);
185  ASSERT( ipLow <= 3);
186  ASSERT( ipHi > 1 );
187  ASSERT( ipHi <=6 );
188 
189  /* set quantum numbers and stat. weights of the transitions: */
190  if( ipHi == 6 )
191  {
192  /* upper is n=3 then set level, stat. weight */
193  QuanNUp = 3.;
194  gHi = 10.;
195  /* following will be set here even though it is not used in this case,
196  * to prevent good compilers from falsing on i not set,
197  * there is assert when used to make sure it is ok */
198  i = -1;
199  }
200  else if( ipHi == 5 )
201  {
202  /* upper is n=3 then set level, stat. weight */
203  QuanNUp = 3.;
204  gHi = 6.;
205  /* following will be set here even though it is not used in this case,
206  * to prevent good compilers from falsing on i not set,
207  * there is assert when used to make sure it is ok */
208  i = -1;
209  }
210  else if( ipHi == 4 )
211  {
212  /* upper is n=3 then set level, stat. weight */
213  QuanNUp = 3.;
214  gHi = 2.;
215  /* following will be set here even though it is not used in this case,
216  * to prevent good compilers from falsing on i not set,
217  * there is assert when used to make sure it is ok */
218  i = -1;
219  }
220  else if( ipHi == 3 )
221  {
222  /* upper is nl=2p then set level, stat. weight */
223  QuanNUp = 2.;
224  gHi = 6.;
225  /* used to point within vectors defined above */
226  i = 0;
227  }
228  else if( ipHi == 2 )
229  {
230  /* upper is nl=2s then set level, stat. weight */
231  QuanNUp = 2.;
232  gHi = 2.;
233  /* used to point within vectors defined above */
234  i = 1;
235  }
236  else
237  {
238  fprintf( ioQQQ, " Insane levels in Hydcs123\n" );
240  }
241 
242  /* which lower level? */
243  if( ipLow == 1 )
244  {
245  /* lower is n=1 then set level, stat. weight */
246  QuanNLo = 1.;
247  gLo = 2.;
248  }
249  else if( ipLow == 2 )
250  {
251  /* lower is nl=2s then set level, stat. weight */
252  QuanNLo = 2.;
253  gLo = 2.;
254  }
255  else if( ipLow == 3 )
256  {
257  QuanNLo = 2.;
258  gLo = 6.;
259  }
260  else
261  {
262  fprintf( ioQQQ, " Insane levels in Hydcs123\n" );
264  }
265 
266  /* following is the physical charge */
267  Charge = (double)(nelem + 1);
268  /* square of charge */
269  ChargeSquared = Charge*Charge;
270 
271  if( ipLow == 2 && ipHi == 3 )
272  {
273  /*************************************** this is 2p-2s:
274  * series of if statements determines which entries Charge is between. */
275  if( nelem < 1 )
276  {
277  /* this can't happen since routine returned above when ip=1 and
278  * special atomic hydrogen routine called */
279  fprintf( ioQQQ, " insane charge given to Hydcs123\n" );
281  }
282  else if( nelem == 1 )
283  {
284  zlow = 2.;
285  j = 1;
286  zhigh = 2.;
287  k = 1;
288  }
289  /* Li through C */
290  else if( nelem <= 5 )
291  {
292  zlow = 2.;
293  j = 1;
294  zhigh = 6.;
295  k = 2;
296  }
297  else if( nelem <= 11 )
298  {
299  zlow = 6.;
300  j = 2;
301  zhigh = 12.;
302  k = 3;
303  }
304  else if( nelem <= 15 )
305  {
306  zlow = 12.;
307  j = 3;
308  zhigh = 16.;
309  k = 4;
310  }
311  else if( nelem <= 17 )
312  {
313  zlow = 16.;
314  j = 4;
315  zhigh = 18.;
316  k = 5;
317  }
318  /* following changed to else from else if,
319  * to prevent false comment in good compilers */
320  /*else if( nelem > 18 )*/
321  else
322  {
323  zlow = 18.;
324  j = 5;
325  zhigh = 18.;
326  k = 5;
327  }
328 
329  /* convert Te to a.u./Z^2
330  * determine rate at the low Charge */
331  x = EVRYD/TE1RYD*phycon.te/(27.211396*pow2(zlow));
332  TeUse = MIN2(x,0.80);
333  x = MAX2(0.025,TeUse);
334  logx = log(x);
335 
336  /* what type of collision are we dealing with? */
337  if( chType == 'e' )
338  {
339  /* electron collisions */
340  Ratelow = ae[j-1] + be[j-1]*x + ce[j-1]*pow2(x)*logx + de[j-1]*
341  exp(x) + ee[j-1]*logx/pow2(x);
342  }
343  else if( chType == 'p' )
344  {
345  Ratelow = ap[j-1] + bp[j-1]*x + cp[j-1]*pow2(x)*logx + dp[j-1]*
346  exp(x) + ep[j-1]*logx/pow2(x);
347  }
348  else
349  {
350  /* this can't happen */
351  fprintf( ioQQQ, " insane collision species given to Hydcs123\n" );
353  }
354 
355  if( fp_equal( zlow, zhigh ) )
356  {
357  rate = Ratelow;
358  }
359  else
360  {
361  /* determine rate at the high Charge */
362  x = EVRYD/TE1RYD*phycon.te/(27.211396*pow2(zhigh));
363  TeUse = MIN2(x,0.80);
364  x = MAX2(0.025,TeUse);
365  logx = log(x);
366  if( chType == 'e' )
367  {
368  Ratehigh = ae[k-1] + be[k-1]*x + ce[k-1]*pow2(x)*logx +
369  de[k-1]*exp(x) + ee[k-1]*logx/pow2(x);
370  }
371  else
372  {
373  Ratehigh = ap[k-1] + bp[k-1]*x + cp[k-1]*pow2(x)*logx +
374  dp[k-1]*exp(x) + ep[k-1]*logx/pow2(x);
375  }
376  /* linearly interpolate in charge */
377  slope = (Ratehigh - Ratelow)/(zhigh - zlow);
378  rate = slope*(Charge - zlow) + Ratelow;
379  }
380  rate = rate/ChargeSquared/Charge*1.0e-7;
381  /* must convert to cs and need to know the valid temp range */
382  templow = 0.025*27.211396*TE1RYD/EVRYD*ChargeSquared;
383  temphigh = 0.80*27.211396*TE1RYD/EVRYD*ChargeSquared;
384  TeUse = MIN2((double)phycon.te,temphigh);
385  temp = MAX2(TeUse,templow);
386  Hydcs123_v = rate*gHi*sqrt(temp)/COLL_CONST;
387 
388  if( chType == 'p' )
389  {
390  /* COLL_CONST is incorrect for protons, correct here */
391  Hydcs123_v *= powpq( PROTON_MASS/ELECTRON_MASS, 3, 2 );
392  }
393  }
394  else if( ipHi == 4 || ipHi == 5 || ipHi == 6 )
395  {
396  /* n = 3
397  * for the rates involving n = 3, must do something different. */
398  if( nelem < 1 )
399  {
400  fprintf( ioQQQ, " insane charge given to Hydcs123\n" );
402  }
403  else if( nelem == 1 )
404  {
405  zlow = 2.;
406  Ratelow = He2cs123(ipLow,ipHi);
407  zhigh = 2.;
408  Ratehigh = Ratelow;
409  }
410  else if( nelem > 1 && nelem <= 5 )
411  {
412  zlow = 2.;
413  Ratelow = He2cs123(ipLow,ipHi);
414  zhigh = 6.;
415  Ratehigh = C6cs123(ipLow,ipHi);
416  }
417  else if( nelem > 5 && nelem <= 9 )
418  {
419  zlow = 6.;
420  Ratelow = C6cs123(ipLow,ipHi);
421  zhigh = 10.;
422  Ratehigh = Ne10cs123(ipLow,ipHi);
423  }
424  else if( nelem > 9 && nelem <= 19 )
425  {
426  zlow = 10.;
427  Ratelow = Ne10cs123(ipLow,ipHi);
428  zhigh = 20.;
429  Ratehigh = Ca20cs123(ipLow,ipHi);
430  }
431  else if( nelem > 19 && nelem <= 25 )
432  {
433  zlow = 20.;
434  Ratelow = Ca20cs123(ipLow,ipHi);
435  zhigh = 26.;
436  Ratehigh = Fe26cs123(ipLow,ipHi);
437  }
438  /*>>>chng 98 dec 17, to else to stop comment from good compilers*/
439  /*else if( nelem > 26 )*/
440  else
441  {
442  Charge = 26.;
443  zlow = 26.;
444  Ratelow = Fe26cs123(ipLow,ipHi);
445  zhigh = 26.;
446  Ratehigh = Ratelow;
447  }
448 
449  /* linearly interpolate */
450  if( fp_equal( zlow, zhigh ) )
451  {
452  rate = Ratelow;
453  }
454  else
455  {
456  slope = (Ratehigh - Ratelow)/(zhigh - zlow);
457  rate = slope*(Charge - zlow) + Ratelow;
458  }
459 
461  Hydcs123_v = rate;
462 
463  /* the routines C6cs123, Ne10cs123, etc... do not resolve L for n>2 */
464  /* dividing by N should roughly recover the original l-resolved data */
465  Hydcs123_v /= 3.;
466  }
467  else
468  {
469  /* this branch 1-2s, 1-2p */
470  if( nelem == 1 )
471  {
472  /* this brance for helium, then return */
473  Hydcs123_v = He2cs123(ipLow,ipHi);
474  return( Hydcs123_v );
475  }
476 
477  /* electron temperature in eV */
478  tev = phycon.te / EVDEGK;
479  /* energy in eV for hydrogenic species and these quantum numbers */
480  EE = ChargeSquared*EVRYD*(1./QuanNLo/QuanNLo - 1./QuanNUp/QuanNUp);
481  /* EE/kT for this transion */
482  q = EE/tev;
483  TeUse = MIN2(q,10.);
484  /* q is now EE/kT but between 1 and 10 */
485  q = MAX2(1.,TeUse);
486  expq = exp(q);
487 
488  /* i must be 0 or 1 */
489  ASSERT( i==0 || i==1 );
490  C = C0[i] + C1[i]/Charge + C2[i]/ChargeSquared;
491  D = D0[i] + D1[i]/Charge + D2[i]/ChargeSquared;
492 
493  /* following code changed so that e1 always returns e1,
494  * orifinal version only returned e1 for x < 1 */
495  /* use disabled e1: */
496  /*if( q < 1. )*/
497  /*{*/
498  /*rate = (B[i-1] + D*q)*exp(-q) + (A[i-1] + C*q - D*q*q)**/
499  /*e1(q);*/
500  /*}*/
501  /*else*/
502  /*{*/
503  /*rate = (B[i-1] + D*q) + (A[i-1] + C*q - D*q*q)*e1(q);*/
504  /*}*/
505  /*rate *= 8.010e-8/2./ChargeSquared/tev*sqrt(tev);*/
506  /* convert to de-excitation */
507  /*if( q < 1. )*/
508  /*{*/
509  /*rate = rate*exp(q)*gLo/gHi;*/
510  /*}*/
511  /*else*/
512  /*{*/
513  /*rate = rate*gLo/gHi;*/
514  /*}*/
515 
516  /*>>>chng 98 dec 17, e1 always returns e1 */
517  rate = (B[i] + D*q)/expq + (A[i] + C*q - D*q*q)*
518  e1(q);
519  // Equation (2) of Callaway 1983 -- can simplify by cancelling tev factors
520  // with R_0(T)
521  //rate *= 8.010e-8/2./ChargeSquared/tev*sqrt(tev);
522  rate *= 8.010e-8/(2. * ChargeSquared * sqrt(tev));
523  /* convert to de-excitation */
524  rate *= expq*gLo/gHi;
525 
526  /* convert to cs */
527  Hydcs123_v = rate*gHi*phycon.sqrte/COLL_CONST;
528  }
529  return( Hydcs123_v );
530 }
531 
532 /*C6cs123 line collision rates for lower levels of hydrogenic carbon, n=1,2,3 */
533 STATIC double C6cs123(long int i,
534  long int j)
535 {
536  double C6cs123_v,
537  TeUse,
538  t,
539  logx,
540  x;
541  static const double a[3] = {-92.23774,-1631.3878,-6326.4947};
542  static const double b[3] = {-11.93818,-218.3341,-849.8927};
543  static const double c[3] = {0.07762914,1.50127,5.847452};
544  static const double d[3] = {78.401154,1404.8475,5457.9291};
545  static const double e[3] = {332.9531,5887.4263,22815.211};
546 
547  DEBUG_ENTRY( "C6cs123()" );
548 
549  /* These are fits to Table 5 of
550  * >>refer c6 cs Aggarwal, K.M., & Kingston, A.E. 1991, J Phys B, 24, 4583
551  * C VI collision rates for 1s-3l, 2s-3l, and 2p-3l,
552  * principal quantum numbers n and l.
553  *
554  * i is the lower level and runs from 1 to 3 (1s, 2s, 2p)
555  * j is the upper level and runs from 2 to 6 (2s, 2p, 3s, 3p, 3d)
556  * 1s-2s,2p is not done here.
557  * check temperature: fits only good between 3.8 < log Te < 6.2
558  */
559  /* arrays for fits of 3 transitions see the code below for key: */
560 
561  TeUse = MAX2(phycon.te,6310.);
562  t = MIN2(TeUse,1.6e6);
563  x = log10(t);
564  logx = log(x);
565 
566  if( i == 1 && j == 2 )
567  {
568  /* 1s - 2s (first entry) */
569  fprintf( ioQQQ, " Carbon VI 2s-1s not done in C6cs123\n" );
571  }
572 
573  else if( i == 1 && j == 3 )
574  {
575  /* 1s - 2p (second entry) */
576  fprintf( ioQQQ, " Carbon VI 2p-1s not done in C6cs123\n" );
578  }
579 
580  else if( i == 1 && ( j == 4 || j == 5 || j == 6 ) )
581  {
582  /* 1s - 3 (first entry) */
583  C6cs123_v = a[0] + b[0]*x + c[0]*pow2(x)*sqrt(x) + d[0]*logx +
584  e[0]*logx/pow2(x);
585  }
586  else if( i == 2 && ( j == 4 || j == 5 || j == 6 ) )
587  {
588  /* 2s - 3 (second entry) */
589  C6cs123_v = a[1] + b[1]*x + c[1]*pow2(x)*sqrt(x) + d[1]*logx +
590  e[1]*logx/pow2(x);
591  }
592  else if( i == 3 && ( j == 4 || j == 5 || j == 6 ) )
593  {
594  /* 2p - 3s (third entry) */
595  C6cs123_v = a[2] + b[2]*x + c[2]*pow2(x)*sqrt(x) + d[2]*logx +
596  e[2]*logx/pow2(x);
597  }
598  else
599  {
600  fprintf( ioQQQ, " insane levels for C VI n=1,2,3 !!!\n" );
602  }
603  return( C6cs123_v );
604 }
605 
606 /*Ca20cs123 line collision rates for lower levels of hydrogenic calcium, n=1,2,3 */
607 STATIC double Ca20cs123(long int i,
608  long int j)
609 {
610  double Ca20cs123_v,
611  TeUse,
612  t,
613  logx,
614  x;
615  static const double a[3] = {-12.5007,-187.2303,-880.18896};
616  static const double b[3] = {-1.438749,-22.17467,-103.1259};
617  static const double c[3] = {0.008219688,0.1318711,0.6043752};
618  static const double d[3] = {10.116516,153.2650,717.4036};
619  static const double e[3] = {45.905343,685.7049,3227.2836};
620 
621  DEBUG_ENTRY( "Ca20cs123()" );
622 
623  /*
624  * These are fits to Table 5 of
625  * >>refer ca20 cs Aggarwal, K.M., & Kingston, A.E. 1992, J Phys B, 25, 751
626  * Ca XX collision rates for 1s-3l, 2s-3l, and 2p-3l,
627  * principal quantum numbers n and l.
628  *
629  * i is the lower level and runs from 1 to 3 (1s, 2s, 2p)
630  * j is the upper level and runs from 2 to 6 (2s, 2p, 3s, 3p, 3d)
631  * 1s-2s,2p is not done here.
632  * check temperature: fits only good between 5.0 < log Te < 7.2
633  */
634 
635  /* arrays for fits of 3 transitions see the code below for key: */
636 
637  TeUse = MAX2(phycon.te,1.0e5);
638  t = MIN2(TeUse,1.585e7);
639  x = log10(t);
640  logx = log(x);
641 
642  if( i == 1 && j == 2 )
643  {
644  /* 1s - 2s (first entry) */
645  fprintf( ioQQQ, " Ca XX 2s-1s not done in Ca20cs123\n" );
647  }
648 
649  else if( i == 1 && j == 3 )
650  {
651  /* 1s - 2p (second entry) */
652  fprintf( ioQQQ, " Ca XX 2p-1s not done in Ca20cs123\n" );
654  }
655 
656  else if( i == 1 && ( j == 4 || j == 5 || j == 6 ))
657  {
658  /* 1s - 3 (first entry) */
659  Ca20cs123_v = a[0] + b[0]*x + c[0]*pow2(x)*sqrt(x) + d[0]*logx +
660  e[0]*logx/pow2(x);
661  }
662  else if( i == 2 && ( j == 4 || j == 5 || j == 6 ))
663  {
664  /* 2s - 3 (second entry) */
665  Ca20cs123_v = a[1] + b[1]*x + c[1]*pow2(x)*sqrt(x) + d[1]*logx +
666  e[1]*logx/pow2(x);
667  }
668  else if( i == 3 && ( j == 4 || j == 5 || j == 6 ))
669  {
670  /* 2p - 3s (third entry) */
671  Ca20cs123_v = a[2] + b[2]*x + c[2]*pow2(x)*sqrt(x) + d[2]*logx +
672  e[2]*logx/pow2(x);
673  }
674  else
675  {
676  fprintf( ioQQQ, " insane levels for Ca XX n=1,2,3 !!!\n" );
678  }
679  return( Ca20cs123_v );
680 }
681 
682 /*Ne10cs123 line collision rates for lower levels of hydrogenic neon, n=1,2,3 */
683 STATIC double Ne10cs123(long int i,
684  long int j)
685 {
686  double Ne10cs123_v,
687  TeUse,
688  logx,
689  t,
690  x;
691  static const double a[3] = {3.346644,151.2435,71.7095};
692  static const double b[3] = {0.5176036,20.05133,13.1543};
693  static const double c[3] = {-0.00408072,-0.1311591,-0.1099238};
694  static const double d[3] = {-3.064742,-129.8303,-71.0617};
695  static const double e[3] = {-11.87587,-541.8599,-241.2520};
696 
697  DEBUG_ENTRY( "Ne10cs123()" );
698 
699  /* These are fits to Table 5 of
700  * >>refer ne10 cs Aggarwal, K.M., & Kingston, A.E. 1991, PhyS, 44, 517
701  * Ne X collision rates for 1s-3, 2s-3l, and 2p-3l,
702  * principal quantum numbers n and l.
703  *
704  * i is the lower level and runs from 1 to 3 (1s, 2s, 2p)
705  * j is the upper level and runs from 2 to 6 (2s, 2p, 3s, 3p, 3d)
706  * 1s-2s,2p is not done here.
707  * check temperature: fits only good between 3.8 < log Te < 6.2
708  * */
709  /* arrays for fits of 3 transitions see the code below for key: */
710 
711  TeUse = MAX2(phycon.te,6310.);
712  t = MIN2(TeUse,1.6e6);
713  x = log10(t);
714  logx = log(x);
715 
716  if( i == 1 && j == 2 )
717  {
718  /* 1s - 2s (first entry) */
719  fprintf( ioQQQ, " Neon X 2s-1s not done in Ne10cs123\n" );
721  }
722 
723  else if( i == 1 && j == 3 )
724  {
725  /* 1s - 2p (second entry) */
726  fprintf( ioQQQ, " Neon X 2p-1s not done in Ne10cs123\n" );
728  }
729 
730  else if( i == 1 && ( j == 4 || j == 5 || j == 6 ) )
731  {
732  /* 1s - 3 (first entry) */
733  Ne10cs123_v = a[0] + b[0]*x + c[0]*pow2(x)*sqrt(x) + d[0]*logx +
734  e[0]*logx/pow2(x);
735  }
736  else if( i == 2 && ( j == 4 || j == 5 || j == 6 ) )
737  {
738  /* 2s - 3 (second entry) */
739  Ne10cs123_v = a[1] + b[1]*x + c[1]*pow2(x)*sqrt(x) + d[1]*logx +
740  e[1]*logx/pow2(x);
741  }
742  else if( i == 3 && ( j == 4 || j == 5 || j == 6 ) )
743  {
744  /* 2p - 3s (third entry) */
745  Ne10cs123_v = a[2] + b[2]*x + c[2]*pow2(x)*sqrt(x) + d[2]*logx +
746  e[2]*logx/pow2(x);
747  }
748  else
749  {
750  fprintf( ioQQQ, " insane levels for Ne X n=1,2,3 !!!\n" );
752  }
753  return( Ne10cs123_v );
754 }
755 
756 /*He2cs123 line collision strengths for lower levels of helium ion, n=1,2,3, by K Korista */
757 STATIC double He2cs123(long int i,
758  long int j)
759 {
760  double He2cs123_v,
761  t;
762  static const double a[11]={0.12176209,0.32916723,0.46546497,0.044501688,
763  0.040523277,0.5234889,1.4903214,1.4215094,1.0295881,4.769306,9.7226127};
764  static const double b[11]={0.039936166,2.9711166e-05,-0.020835863,3.0508137e-04,
765  -2.004485e-15,4.41475e-06,1.0622666e-05,2.0538877e-06,0.80638448,2.0967075e-06,
766  7.6089851e-05};
767  static const double c[11]={143284.77,0.73158545,-2.159172,0.43254802,2.1338557,
768  8.9899702e-06,-2.9001451e-12,1.762076e-05,52741.735,-2153.1219,-3.3996921e-11};
769 
770  DEBUG_ENTRY( "He2cs123()" );
771 
772  /* These are fits to Table 2
773  * >>refer he2 cs Aggarwal, K.M., Callaway, J., Kingston, A.E., Unnikrishnan, K.
774  * >>refercon 1992, ApJS, 80, 473
775  * He II collision rates for 1s-2s, 1s-2p, 1s-3s, 1s-3p, 1s-3d, 2s-3s, 2s-3p, 2s-3d,
776  * 2p-3s, 2p-3p, and 2p-3d.
777  * principal quantum numbers n and l.
778  *
779  * i is the lower level and runs from 1 to 3 (1s, 2s, 2p)
780  * j is the upper level and runs from 2 to 6 (2s, 2p, 3s, 3p, 3d)
781  * check temperature: fits only good between 5,000K and 500,000K
782  * */
783  /* array for fits of 11 transitions see the code below for key: */
784 
785  t = phycon.te;
786  if( t < 5000. )
787  {
788  t = 5000.;
789  }
790  else if( t > 5.0e05 )
791  {
792  t = 5.0e05;
793  }
794 
795  /**************fits begin here**************
796  * */
797  if( i == 1 && j == 2 )
798  {
799  /* 1s - 2s (first entry) */
800  He2cs123_v = a[0] + b[0]*exp(-t/c[0]);
801  }
802  else if( i == 1 && j == 3 )
803  {
804  /* 1s - 2p (second entry) */
805  He2cs123_v = a[1] + b[1]*pow(t,c[1]);
806  }
807  else if( i == 1 && j == 4 )
808  {
809  /* 1s - 3s (third entry) */
810  double logt = log(t);
811  He2cs123_v = a[2] + b[2]*logt + c[2]/logt;
812  }
813  else if( i == 1 && j == 5 )
814  {
815  /* 1s - 3p (fourth entry) */
816  He2cs123_v = a[3] + b[3]*pow(t,c[3]);
817  }
818  else if( i == 1 && j == 6 )
819  {
820  /* 1s - 3d (fifth entry) */
821  He2cs123_v = a[4] + b[4]*pow(t,c[4]);
822  }
823  else if( i == 2 && j == 4 )
824  {
825  /* 2s - 3s (sixth entry) */
826  He2cs123_v = (a[5] + c[5]*t)/(1 + b[5]*t);
827  }
828  else if( i == 2 && j == 5 )
829  {
830  /* 2s - 3p (seventh entry) */
831  He2cs123_v = a[6] + b[6]*t + c[6]*t*t;
832  }
833  else if( i == 2 && j == 6 )
834  {
835  /* 2s - 3d (eighth entry) */
836  He2cs123_v = (a[7] + c[7]*t)/(1 + b[7]*t);
837  }
838  else if( i == 3 && j == 4 )
839  {
840  /* 2p - 3s (ninth entry) */
841  He2cs123_v = a[8] + b[8]*exp(-t/c[8]);
842  }
843  else if( i == 3 && j == 5 )
844  {
845  /* 2p - 3p (tenth entry) */
846  He2cs123_v = a[9] + b[9]*t + c[9]/t;
847  }
848  else if( i == 3 && j == 6 )
849  {
850  /* 2p - 3d (eleventh entry) */
851  He2cs123_v = a[10] + b[10]*t + c[10]*t*t;
852  }
853  else
854  {
855  /**************fits end here************** */
856  fprintf( ioQQQ, " insane levels for He II n=1,2,3 !!!\n" );
858  }
859  return( He2cs123_v );
860 }
861 
862 /*Fe26cs123 line collision rates for lower levels of hydrogenic iron, n=1,2,3 */
863 STATIC double Fe26cs123(long int i,
864  long int j)
865 {
866  double Fe26cs123_v,
867  TeUse,
868  logx,
869  t,
870  x;
871  static const double a[3] = {-4.238398,-238.2599,-1211.5237};
872  static const double b[3] = {-0.4448177,-27.06869,-136.7659};
873  static const double c[3] = {0.0022861,0.153273,0.7677703};
874  static const double d[3] = {3.303775,191.7165,972.3731};
875  static const double e[3] = {15.82689,878.1333,4468.696};
876 
877  DEBUG_ENTRY( "Fe26cs123()" );
878 
879  /* These are fits to Table 5 of
880  * >>refer fe26 cs Aggarwal, K.M., & Kingston, A.E. 1993, ApJS, 85, 187
881  * Fe XXVI collision rates for 1s-3, 2s-3, and 2p-3,
882  * principal quantum numbers n and l.
883  *
884  * i is the lower level and runs from 1 to 3 (1s, 2s, 2p)
885  * j is the upper level and runs from 2 to 6 (2s, 2p, 3s, 3p, 3d)
886  * 1s-2s,2p is not done here.
887  * check temperature: fits only good between 5.2 < log Te < 7.2
888  * */
889  /* arrays for fits of 3 transitions see the code below for key: */
890 
891  TeUse = MAX2(phycon.te,1.585e5);
892  t = MIN2(TeUse,1.585e7);
893  x = log10(t);
894  logx = log(x);
895 
896  if( i == 1 && j == 2 )
897  {
898  /* 1s - 2s (first entry) */
899  fprintf( ioQQQ, " Fe XXVI 2s-1s not done in Fe26cs123\n" );
901  }
902 
903  else if( i == 1 && j == 3 )
904  {
905  /* 1s - 2p (second entry) */
906  fprintf( ioQQQ, " Fe XXVI 2p-1s not done in Fe26cs123\n" );
908  }
909 
910  else if( i == 1 && ( j == 4 || j == 5 || j == 6 ) )
911  {
912  /* 1s - 3 (first entry) */
913  Fe26cs123_v = a[0] + b[0]*x + c[0]*pow2(x)*sqrt(x) + d[0]*logx +
914  e[0]*logx/pow2(x);
915  }
916  else if( i == 2 && ( j == 4 || j == 5 || j == 6 ) )
917  {
918  /* 2s - 3 (second entry) */
919  Fe26cs123_v = a[1] + b[1]*x + c[1]*pow2(x)*sqrt(x) + d[1]*logx +
920  e[1]*logx/pow2(x);
921  }
922  else if( i == 3 && ( j == 4 || j == 5 || j == 6 ) )
923  {
924  /* 2p - 3s (third entry) */
925  Fe26cs123_v = a[2] + b[2]*x + c[2]*pow2(x)*sqrt(x) + d[2]*logx +
926  e[2]*logx/pow2(x);
927  }
928  else
929  {
930  fprintf( ioQQQ, " insane levels for Ca XX n=1,2,3 !!!\n" );
932  }
933  return( Fe26cs123_v );
934 }
935 
936 
937 STATIC double CS_ThermAve_PR78(long ipISO, long nelem, long nHi, long nLo, double deltaE, double temp )
938 {
939 
940  double coll_str;
941 
942  DEBUG_ENTRY( "CS_ThermAve_PR78()" );
943 
944  global_ipISO = ipISO;
945  global_nelem = nelem;
946  global_nHi = nHi;
947  global_nLo = nLo;
948  global_deltaE = deltaE;
949 
950  kTRyd = temp / TE1RYD;
951 
952  if( !iso_ctrl.lgCS_therm_ave[ipISO] )
953  {
954  /* Do NOT average over Maxwellian */
955  coll_str = CS_PercivalRichards78( kTRyd );
956  }
957  else
958  {
959  /* DO average over Maxwellian */
960  coll_str = qg32( 0.0, 1.0, Therm_ave_coll_str_int_PR78);
961  coll_str += qg32( 1.0, 10.0, Therm_ave_coll_str_int_PR78);
962  }
963 
964  return coll_str;
965 }
966 
967 /* The integrand for calculating the thermal average of collision strengths */
968 STATIC double Therm_ave_coll_str_int_PR78( double EOverKT )
969 {
970  double integrand;
971 
972  DEBUG_ENTRY( "Therm_ave_coll_str_int_PR78()" );
973 
974  integrand = exp( -1.*EOverKT ) * CS_PercivalRichards78( EOverKT * kTRyd );
975 
976  return integrand;
977 }
978 
979 inline double C2_PR78(double x, double y)
980 {
981  return pow2(x)*log1p(2./3.*x)/(2.*y + 1.5*x);
982 }
983 
984 STATIC double CS_PercivalRichards78( double Ebar )
985 {
986  DEBUG_ENTRY( "CS_PercivalRichards78()" );
987 
988  if( Ebar < global_deltaE )
989  return 0.;
990 
991  long ipISO = global_ipISO;
992  long nelem = global_nelem;
993  double np = (double)global_nHi;
994  double n = (double)global_nLo;
995  double n2 = pow2(n);
996  double nnp = n*np;
997  double Ebar2 = pow2(Ebar);
998  double s = np - n;
999  long is = global_nHi - global_nLo ;
1000 
1001  ASSERT( s > 0. );
1002 
1003  double A = 8./(3.*s) * pow3(np/(s*n)) * (0.184 - 0.04/pow2(cbrt(s))) * powi(1.-0.2*s/nnp, 1+2*is);
1004 
1005  double Z3 = (double)(nelem + 1 - ipISO);
1006  double Z32 = pow2(Z3);
1007 
1008  double D = exp(-Z32/(nnp*Ebar2));
1009 
1010  double L = log((1.+0.53*Ebar2*nnp/Z32) / (1.+0.4*Ebar));
1011 
1012  double F = powi(1.-0.3*s*D/nnp, 1+2*is);
1013 
1014  double G = 0.5*pow3(Ebar*n2/(Z3*np));
1015 
1016  double h1 = sqrt(2. - pow2(n/np));
1017  double h2 = 2.*Z3/(n2*Ebar);
1018  double xPlus = h2/(h1+1.);
1019  double xMinus = h2/(h1-1.);
1020 
1021  double y = 1./(1.-D*log(18.*s)/(4.*s));
1022 
1023  double H = C2_PR78(xMinus,y) - C2_PR78(xPlus,y);
1024 
1025  /* this is the LHS of equation 1 of PR78 */
1026  double cross_section = A*D*L + F*G*H;
1027  /* this is the result after solving equation 1 for the cross section */
1028  cross_section *= PI*pow2(n2*BOHR_RADIUS_CM/Z3) / Ebar;
1029 
1030  double stat_weight;
1031  if( ipISO == ipH_LIKE )
1032  stat_weight = 2.*n2;
1033  else if( ipISO == ipHE_LIKE )
1034  stat_weight = 4.*n2;
1035  else
1036  TotalInsanity();
1037 
1038  /* convert to collision strength */
1039  return cross_section*stat_weight*Ebar / (PI*pow2(BOHR_RADIUS_CM));
1040 }
1041 
1042 #if 0
1043 STATIC void TestPercivalRichards( void )
1044 {
1045  double CStemp;
1046 
1047  /* this reproduces Table 1 of PR78 */
1048  for( long i=0; i<5; i++ )
1049  {
1050  double Ebar[5] = {0.1, 0.4, 0.8, 1.0, 10.};
1051 
1052  CStemp = CS_PercivalRichards78( 0, 2, 12, 10, Ebar[i] );
1053  }
1054 
1055  /* this reproduces Table 2 of PR78 */
1056  for( long i=0; i<5; i++ )
1057  {
1058  double Ebar[5] = {0.1, 0.4, 0.8, 1.0, 10.};
1059 
1060  CStemp = CS_ThermAve_PR78( ipISO, 0, N_(ipHi), N_(ipLo), phycon.te );
1061  }
1062 
1063  return;
1064 }
1065 #endif
1066 
1068  long ipHi,
1069  long ipLo,
1070  long ipCollider )
1071 {
1072  DEBUG_ENTRY( "HydroCSInterp()" );
1073 
1074  long nHi = iso_sp[ipH_LIKE][nelem].st[ipHi].n();
1075  long lHi = iso_sp[ipH_LIKE][nelem].st[ipHi].l();
1076  long sHi = iso_sp[ipH_LIKE][nelem].st[ipHi].S();
1077  long gHi = iso_sp[ipH_LIKE][nelem].st[ipHi].g();
1078  double IP_Ryd_Hi = iso_sp[ipH_LIKE][nelem].fb[ipHi].xIsoLevNIonRyd;
1079  long nLo = iso_sp[ipH_LIKE][nelem].st[ipLo].n();
1080  long lLo = iso_sp[ipH_LIKE][nelem].st[ipLo].l();
1081  long sLo = iso_sp[ipH_LIKE][nelem].st[ipLo].S();
1082  long gLo = iso_sp[ipH_LIKE][nelem].st[ipLo].g();
1083  double IP_Ryd_Lo = iso_sp[ipH_LIKE][nelem].fb[ipLo].xIsoLevNIonRyd;
1084  // collisions are from high to low level, then initial level lifetime is from higher level
1085  double tauLo = iso_sp[ipH_LIKE][nelem].st[ipHi].lifetime();
1086  double EnerErg = iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).EnergyErg();
1087  const char *where=" ";
1088 
1089  double CStemp = GetHlikeCollisionStrength( nelem, ipCollider,
1090  nHi, lHi, sHi, gHi, IP_Ryd_Hi,
1091  nLo, lLo, sLo, gLo, IP_Ryd_Lo,
1092  tauLo, EnerErg, &where );
1093  return (realnum)CStemp;
1094 }
1095 
1096 realnum GetHlikeCollisionStrength( long nelem, long ipCollider,
1097  long nHi, long lHi, long sHi, long gHi, double IP_Ryd_Hi,
1098  long nLo, long lLo, long sLo, long gLo, double IP_Ryd_Lo,
1099  double tauLo, double EnerErg, const char **where )
1100 {
1101  DEBUG_ENTRY( "GetHlikeCollisionStrength()" );
1102 
1104  return 0.;
1105 
1106  double CStemp;
1107 
1108  // Energy gap between levels in eV
1109  double deltaE_eV = EnerErg/EN1EV;
1110 
1111  *where = " ";
1112  /* This set of collision strengths should only be used
1113  * if the Storey and Hummer flag is set */
1115  {
1116  if( nLo == nHi )
1117  {
1118  if( nHi <= iso_sp[ipH_LIKE][nelem].n_HighestResolved_max &&
1119  abs( lLo - lHi ) != 1 )
1120  {
1121  /* if delta L is not +/- 1, set collision strength to zero. */
1122  CStemp = 0.;
1123  }
1124  else if (ipCollider == ipELECTRON)
1125  {
1126  /* l-changing collisions are calculated for heavy slow projectile */
1127  CStemp = 0.;
1128  }
1129  else
1130  {
1131  //classical PS64 for Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165
1132  CStemp = CS_l_mixing_PS64(
1133  nelem,
1134  tauLo,
1135  nelem+1.-ipH_LIKE,
1136  nHi,
1137  lHi,
1138  gHi,
1139  lLo,
1140  deltaE_eV,
1141  ipCollider);
1142  if (0)
1143  {
1144  // New PS-M Refer F. Guzman MNRAS (2016)
1145  CStemp = CS_l_mixing_PS64_expI(
1146  nelem,
1147  tauLo,
1148  nelem+1.-ipH_LIKE,
1149  nHi,
1150  lHi,
1151  gHi,
1152  lLo,
1153  deltaE_eV,
1154  ipCollider);
1155  }
1156  }
1157  }
1158 
1159  else
1160  {
1161  if( nHi <= iso_sp[ipH_LIKE][nelem].n_HighestResolved_max &&
1162  abs( lLo - lHi ) != 1 )
1163  {
1164  /* if delta L is not +/- 1, set collision strength to zero. */
1165  CStemp = 0.;
1166  }
1167  else if( ipCollider == ipELECTRON )
1168  {
1169  CStemp = CS_ThermAve_PR78( ipH_LIKE, nelem, nHi, nLo,
1170  EnerErg / EN1RYD , phycon.te );
1171  if ( lHi >= 0 )
1172  CStemp /= 2.*pow2(nHi);
1173  }
1174  else
1175  CStemp = 0.;
1176  }
1177  }
1178  else
1179  {
1180  /* HCSAR_interp interpolates on a table to return R-matrix collision strengths
1181  * for hydrogen only */
1182  if( (CStemp = HlikeCSInterp( nelem, ipCollider, nHi, lHi, sHi, nLo, lLo, sLo )) >= 0.f )
1183  {
1184  fixit( "do something here and throughout this routine with where as in helike_cs.cpp." );
1185  *where = "HlikeI";
1186  }
1187  else if( nelem==ipHYDROGEN && ipCollider==ipELECTRON && ( nHi != nLo ) )
1188  {
1189  double Aul_j = HydroEinstA(nHi,nLo);
1190  CStemp = hydro_vs_deexcit( nHi, gHi, IP_Ryd_Hi, nLo, gLo, IP_Ryd_Lo, Aul_j );
1191  *where = "Vriens";
1192  /* There are NO resolved data for levels higher than n=4 */
1193  if ( lHi >= 0 )
1194  CStemp /= 2.*pow2(nHi);
1195  }
1196  else if( nLo == nHi )
1197  {
1198  if ( ipCollider == ipELECTRON)
1199  /* l-changing collisions are calculated for heavy slow projectile */
1200  CStemp = 0.;
1201  else if( iso_ctrl.lgCS_Vrinceanu[ipH_LIKE] )
1202  {
1203  CStemp = CS_l_mixing_VF01(ipH_LIKE, nelem,
1204  nHi,
1205  lHi,
1206  lLo,
1207  sLo,
1208  gHi,
1209  tauLo,
1210  IP_Ryd_Hi,
1211  IP_Ryd_Lo,
1212  phycon.te,
1213  ipCollider );
1214  *where ="VF01 ";
1215  }
1216  else if( iso_ctrl.lgCS_VOS12[ipH_LIKE] )
1217  {
1218  /* Use the method given in
1219  * >>refer He CS Vrinceau etal ApJ 747, 56 (2012); semiclasical treatment: eq. (7)
1220  * statistical weights included */
1221  CStemp = CS_l_mixing_VOS12(nHi,lHi,lLo,nelem,gHi,nelem+1.-ipH_LIKE,ipCollider,phycon.sqrte);
1222  *where = "VOS12 ";
1223  }
1224  else if( iso_ctrl.lgCS_VOS12QM[ipH_LIKE] )// && abs( lLo - lHi ) == 1)
1225  {
1226  /* Use the method given in
1227  * >>refer He CS Vrinceau etal ApJ 747, 56 (2012); quantal treatment: eq. (2)
1228  * statistical weights included */
1229  CStemp = CS_l_mixing_VOS12QM(ipH_LIKE, nelem,
1230  nLo,
1231  lHi,
1232  lLo,
1233  sLo,
1234  gHi,
1235  tauLo,
1236  IP_Ryd_Hi,
1237  IP_Ryd_Lo,
1238  phycon.te,
1239  ipCollider );
1240  *where = "VOS12Q";
1241  }
1242  else if ( iso_ctrl.lgCS_PS64[ipH_LIKE] && abs( lLo - lHi ) == 1 )
1243  {
1245  {
1246  //classical PS64 for Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165
1247  CStemp = CS_l_mixing_PS64(
1248  nelem,
1249  tauLo,
1250  nelem+1.-ipH_LIKE,
1251  nHi,
1252  lHi,
1253  gHi,
1254  lLo,
1255  deltaE_eV,
1256  ipCollider);
1257  *where = "PS64 ";
1258  }
1259  else
1260  // PS-M Refer F. Guzman MNRAS (2016)
1261  CStemp = CS_l_mixing_PS64_expI(
1262  nelem,
1263  tauLo,
1264  nelem+1.-ipH_LIKE,
1265  nHi,
1266  lHi,
1267  gHi,
1268  lLo,
1269  //sHi,
1270  deltaE_eV,
1271  ipCollider);
1272  *where = "PSM ";
1273  }
1274  else
1275  CStemp = 0.;
1276  }
1277  else
1278  {
1279  ASSERT( nHi != nLo );
1280  /* highly excited levels */
1281  /* Vriens&Smeets (1980) give no Z dependence on their cross sections,
1282  * Percival and Richards (1978) have got a Z dependence so their rates are preferred */
1283  CStemp = CS_ThermAve_PR78( ipH_LIKE, nelem, nHi, nLo,
1284  EnerErg / EN1RYD , phycon.te );
1285  *where = "PR78 ";
1286  /* Thse are unresolved rates so they should e averaged for resolved transitions */
1287  if ( lHi >= 0 )
1288  CStemp /= 2.*pow2(nHi);
1289 
1290  }
1291  }
1292 
1293  return (realnum)CStemp;
1294 }
1295 
1296 STATIC realnum HlikeCSInterp( long nelem, long Collider, long nHi, long lHi, long sHi, long nLo, long lLo, long sLo )
1297 {
1298  DEBUG_ENTRY( "HlikeCSInterp()" );
1299 
1300  // Cannot do collapsed levels with unspecified l here.
1301  if( lLo < 0 || lHi < 0 )
1302  return -1.f;
1303 
1304  ASSERT( sHi==2 && sLo==2 );
1305  long ipLo = iso_sp[ipH_LIKE][nelem].QuantumNumbers2Index[nLo][lLo][sLo];
1306  long ipHi = iso_sp[ipH_LIKE][nelem].QuantumNumbers2Index[nHi][lHi][sHi];
1307  ASSERT( ipLo < ipHi );
1308 
1309  realnum cs = -1.f;
1310 
1311  if( nelem==ipHYDROGEN && Collider==ipELECTRON && nHi <= 5 && ( nHi != nLo ) )
1312  {
1313  cs = HCSAR_interp(ipLo,ipHi);
1314  }
1315  else if( nelem>ipHYDROGEN && Collider==ipELECTRON && nHi <= 3 && nLo < 3 )
1316  {
1317  cs = Hydcs123(ipLo + 1,ipHi + 1,nelem,'e');
1318  }
1319  else if( nelem>ipHYDROGEN && Collider==ipPROTON && nHi==2 && lHi==1 && nLo==2 && lLo==0 )
1320  {
1321  cs = Hydcs123(ipLo + 1,ipHi + 1,nelem,'p');
1322  }
1323 
1324  return cs;
1325 }
1326 
static long global_nelem
Definition: hydrocollid.cpp:39
realnum EnergyErg() const
Definition: transition.h:90
qList st
Definition: iso.h:482
STATIC double Ne10cs123(long int i, long int j)
bool lgCS_PSClassic[NISO]
Definition: iso.h:408
const int ipHE_LIKE
Definition: iso.h:65
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
t_opac opac
Definition: opacity.cpp:5
realnum h_coll_str(long ipLo, long ipHi, long ipTe)
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
double e1(double x)
bool lgCS_therm_ave[NISO]
Definition: iso.h:408
double CS_l_mixing_VF01(long ipISO, long nelem, long n, long l, long lp, long s, long gLo, double tauLo, double IP_Ryd_Hi, double IP_Ryd_Lo, double temp, long Collider)
Definition: helike_cs.cpp:1650
STATIC double Fe26cs123(long int i, long int j)
double C2_PR78(double x, double y)
t_phycon phycon
Definition: phycon.cpp:6
double CS_l_mixing_PS64(long nelem, double tau, double target_charge, long n, long l, double gLo, long lp, double deltaE_eV, long Collider)
Definition: helike_cs.cpp:1376
T pow3(T a)
Definition: cddefines.h:988
static long global_nHi
Definition: hydrocollid.cpp:39
static const double C1
FILE * ioQQQ
Definition: cddefines.cpp:7
vector< freeBound > fb
Definition: iso.h:481
#define MIN2(a, b)
Definition: cddefines.h:803
bool lgCS_Vrinceanu[NISO]
Definition: iso.h:408
STATIC double Hydcs123(long int ilow, long int ihigh, long int iz, long int chType)
Definition: hydrocollid.cpp:94
static t_ADfA & Inst()
Definition: cddefines.h:209
realnum GetHlikeCollisionStrength(long nelem, long ipCollider, long nHi, long lHi, long sHi, long gHi, double IP_Ryd_Hi, long nLo, long lLo, long sLo, long gLo, double IP_Ryd_Lo, double tauLo, double EnerErg, const char **where)
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
STATIC double CS_ThermAve_PR78(long ipISO, long nelem, long nHi, long nLo, double deltaE, double temp)
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:854
bool lgColl_excite[NISO]
Definition: iso.h:362
realnum HydroCSInterp(long nelem, long ipHi, long ipLo, long ipCollider)
double CS_l_mixing_VOS12(long n, long l, long lp, long nelem, double gLo, long Ztarget, long Collider, double sqrte)
Definition: helike_cs.cpp:2078
STATIC realnum HCSAR_interp(int ipLo, int ipHi)
Definition: hydrocollid.cpp:45
static const realnum HCSTE[NHCSTE]
Definition: hydrocollid.cpp:42
STATIC double He2cs123(long int i, long int j)
#define STATIC
Definition: cddefines.h:118
#define N_(A_)
Definition: iso.h:22
static long global_ipISO
Definition: hydrocollid.cpp:39
bool lgCS_PS64[NISO]
Definition: iso.h:408
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
STATIC double CS_PercivalRichards78(double Ebar)
#define cdEXIT(FAIL)
Definition: cddefines.h:482
static double kTRyd
Definition: hydrocollid.cpp:40
double powi(double, long int)
Definition: service.cpp:594
diatomics h2("h2", 4100.,&hmi.H2_total, Yan_H2_CS)
multi_arr< long, 3 > QuantumNumbers2Index
Definition: iso.h:490
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:473
STATIC double Ca20cs123(long int i, long int j)
double CS_l_mixing_VOS12QM(long ipISO, long nelem, long n, long l, long lp, long s, long gLo, double tauLo, double IP_Ryd_Hi, double IP_Ryd_Lo, double temp, long Collider)
Definition: helike_cs.cpp:1669
#define ASSERT(exp)
Definition: cddefines.h:613
double HydroEinstA(long int n1, long int n2)
Definition: hydroeinsta.cpp:13
bool lgCS_VOS12[NISO]
Definition: iso.h:408
static double global_deltaE
Definition: hydrocollid.cpp:40
double qg32(double, double, double(*)(double))
Definition: service.cpp:1175
const int ipH_LIKE
Definition: iso.h:64
const int LIMELM
Definition: cddefines.h:308
T pow2(T a)
Definition: cddefines.h:981
STATIC double C6cs123(long int i, long int j)
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
double powpq(double x, int p, int q)
Definition: service.cpp:630
#define NHCSTE
Definition: atmdat.h:268
bool lgCS_VOS12QM[NISO]
Definition: iso.h:408
STATIC double cross_section(double EgammaRyd, double EthRyd, long nelem, long n, long l, long s)
static long global_nLo
Definition: hydrocollid.cpp:39
#define MAX2(a, b)
Definition: cddefines.h:824
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
STATIC realnum HlikeCSInterp(long nelem, long Collider, long nHi, long lHi, long sHi, long nLo, long lLo, long sLo)
double hydro_vs_deexcit(long nHi, long gHi, double IP_Ryd_Hi, long nLo, long gLo, double IP_Ryd_Lo, double Aul)
double sqrte
Definition: phycon.h:58
#define fixit(a)
Definition: cddefines.h:417
double te
Definition: phycon.h:21
double CS_l_mixing_PS64_expI(long nelem, double tau, double target_charge, long n, long l, double g, long lp, double deltaE_eV, long Collider)
Definition: helike_cs.cpp:1183
STATIC double Therm_ave_coll_str_int_PR78(double EOverKT)
const int ipHYDROGEN
Definition: cddefines.h:349
bool lgCaseB_HummerStorey
Definition: opacity.h:178