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

Generated for cloudy by doxygen 1.8.3.1