cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
helike_cs.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 /*HeCollid evaluate collisional rates */
4 /*HeCSInterp interpolate on He1 collision strengths */
5 /*AtomCSInterp do the atom */
6 /*IonCSInterp do the ions */
7 /*CS_l_mixing_PS64 - find rate for l-mixing collisions by protons, for neutrals */
8 #include "cddefines.h"
9 #include "atmdat.h"
10 #include "conv.h"
11 #include "dense.h"
12 #include "helike.h"
13 #include "helike_cs.h"
14 #include "hydro_vs_rates.h"
15 #include "iso.h"
16 #include "opacity.h"
17 #include "phycon.h"
18 #include "physconst.h"
19 #include "rfield.h"
20 #include "taulines.h"
21 #include "thirdparty.h"
22 #include "trace.h"
23 
24 /* returns thermally-averaged Seaton 62 collision strength. */
25 STATIC double S62_Therm_ave_coll_str( double EProjectile_eV );
26 
27 /* all of these are used in the calculation of Stark collision strengths
28  * following the algorithms in Vrinceanu & Flannery (2001). */
29 STATIC double Therm_ave_coll_str_int_VF01( double EProjectileRyd);
30 STATIC double collision_strength_VF01( double velOrEner, bool lgParamIsRedVel );
31 STATIC double L_mix_integrand_VF01( double alpha );
32 STATIC double StarkCollTransProb_VF01( long int n, long int l, long int lp, double alpha, double deltaPhi);
33 
37 static double kTRyd;
38 
39 /* These are masses relative to the proton mass of the electron, proton, he+, and alpha particle. */
40 static double ColliderMass[4] = {ELECTRON_MASS/PROTON_MASS, 1.0, 4.0, 4.0};
41 static double ColliderCharge[4] = {1.0, 1.0, 1.0, 2.0};
42 
43 void HeCollidSetup( void )
44 {
45  /* this must be longer than data path string, set in path.h/cpu.cpp */
46  long i, i1, j, nelem, ipHi, ipLo;
47  bool lgEOL, lgHIT;
48  FILE *ioDATA;
49  long ipISO = ipHE_LIKE;
50 
51 # define chLine_LENGTH 1000
52  char chLine[chLine_LENGTH];
53 
54  DEBUG_ENTRY( "HeCollidSetup()" );
55 
56  /* get the collision strength data for the He 1 lines */
57  if( trace.lgTrace )
58  fprintf( ioQQQ," HeCollidSetup opening he1_cs.dat:");
59 
60  ioDATA = open_data( "he1_cs.dat", "r" );
61 
62  /* check that magic number is ok */
63  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
64  {
65  fprintf( ioQQQ, " HeCollidSetup could not read first line of he1_cs.dat.\n");
66  cdEXIT(EXIT_FAILURE);
67  }
68  i = 1;
69  i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
70  /*i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
71  i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);*/
72 
73  /* the following is to check the numbers that appear at the start of he1_cs.dat */
74  if( i1 !=COLLISMAGIC )
75  {
76  fprintf( ioQQQ,
77  " HeCollidSetup: the version of he1_cs.dat is not the current version.\n" );
78  fprintf( ioQQQ,
79  " HeCollidSetup: I expected to find the number %i and got %li instead.\n" ,
80  COLLISMAGIC, i1 );
81  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
82  cdEXIT(EXIT_FAILURE);
83  }
84 
85  /* get the array of temperatures */
86  lgHIT = false;
87  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
88  {
89  /* only look at lines without '#' in first col */
90  if( chLine[0] != '#')
91  {
92  lgHIT = true;
93  iso.nCS[ipISO] = 0;
94  lgEOL = false;
95  i = 1;
96  while( !lgEOL && iso.nCS[ipISO] < HE1CSARRAY)
97  {
98  iso.CSTemp[ipISO][iso.nCS[ipISO]] = FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
99  ++iso.nCS[ipISO];
100  }
101  break;
102  }
103  }
104  --iso.nCS[ipISO];
105  if( !lgHIT )
106  {
107  fprintf( ioQQQ, " HeCollidSetup could not find line in CS temperatures in he1_cs.dat.\n");
108  cdEXIT(EXIT_FAILURE);
109  }
110 
111  /* create space for array of CS values, if not already done */
112  iso.HeCS.reserve( NISO );
113  iso.HeCS.reserve( ipISO, LIMELM );
114 
115  for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
116  {
120  if( nelem != ipHELIUM )
121  continue;
122 
123  /* only grab core for elements that are turned on */
124  if( nelem == ipHELIUM || dense.lgElmtOn[nelem])
125  {
126  iso.HeCS.reserve( ipISO, nelem, iso.numLevels_max[ipHE_LIKE][nelem]-iso.nCollapsed_max[ipHE_LIKE][nelem] );
127 
128  for( ipHi=ipHe2s3S; ipHi < iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem];++ipHi )
129  {
130  iso.HeCS.reserve( ipISO, nelem, ipHi, ipHi );
131 
132  for( ipLo=ipHe1s1S; ipLo<ipHi; ++ipLo )
133  iso.HeCS.reserve( ipISO, nelem, ipHi, ipLo, iso.nCS[ipISO] );
134  }
135  }
136  }
137 
138  iso.HeCS.alloc();
139  iso.HeCS.zero();
140 
141  /* now read in the CS data */
142  lgHIT = false;
143  nelem = ipHELIUM;
144  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
145  {
146  char *chTemp;
147  /* only look at lines without '#' in first col */
148  if( (chLine[0] == ' ') || (chLine[0]=='\n') )
149  break;
150  if( chLine[0] != '#')
151  {
152  lgHIT = true;
153 
154  /* get lower and upper level index */
155  j = 1;
156  ipLo = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
157  ipHi = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
158  ASSERT( ipLo < ipHi );
159  if( ipHi >= iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] )
160  continue;
161  else
162  {
163  chTemp = chLine;
164  /* skip over 4 tabs to start of cs data */
165  for( i=0; i<3; ++i )
166  {
167  if( (chTemp = strchr( chTemp, '\t' )) == NULL )
168  {
169  fprintf( ioQQQ, " HeCollidSetup no init cs\n" );
170  cdEXIT(EXIT_FAILURE);
171  }
172  ++chTemp;
173  }
174 
175  for( i=0; i<iso.nCS[ipISO]; ++i )
176  {
177  double a;
178  if( (chTemp = strchr( chTemp, '\t' )) == NULL )
179  {
180  fprintf( ioQQQ, " HeCollidSetup not scan cs, current indices: %li %li\n", ipHi, ipLo );
181  cdEXIT(EXIT_FAILURE);
182  }
183  ++chTemp;
184  sscanf( chTemp , "%le" , &a );
185  iso.HeCS[ipISO][nelem][ipHi][ipLo][i] = (realnum)a;
186  }
187  }
188  }
189  }
190 
191  /* close the data file */
192  fclose( ioDATA );
193 
194  return;
195 }
196 
197 /* Choose either AtomCSInterp or IonCSInterp */
198 realnum HeCSInterp(long int nelem,
199  long int ipHi,
200  long int ipLo,
201  long int Collider )
202 {
203  realnum cs, factor1;
204 
205  /* This variable is for diagnostic purposes:
206  * a string used in the output to designate where each cs comes from. */
207  const char *where = " ";
208 
209  DEBUG_ENTRY( "HeCSInterp()" );
210 
211  if( !iso.lgColl_excite[ipHE_LIKE] )
212  {
213  return (realnum)1E-10;
214  }
215 
216  if( nelem == ipHELIUM )
217  {
218  /* do for helium */
219  cs = AtomCSInterp( nelem, ipHi , ipLo, &factor1, &where, Collider );
220  }
221  else
222  {
223  /* get collision strengths for an ion */
224  cs = IonCSInterp( nelem, ipHi , ipLo, &factor1, &where, Collider );
225  }
226 
227  /* set 15% errors on all collision rates. */
228  iso_put_error(ipHE_LIKE, nelem, ipHi, ipLo, IPCOLLIS, 0.15f );
229 
230  ASSERT( cs >= 0.f );
231 
232  /* in many cases the correction factor for split states has already been made,
233  * if not then factor is still negative */
234  /* Remove the second test here when IonCSInterp is up to par with AtomCSInterp */
235  ASSERT( factor1 >= 0.f || nelem!=ipHELIUM );
236  if( factor1 < 0.f )
237  {
239 
240  factor1 = 1.f;
241  }
242 
243  /* take factor into account, usually 1, ratio of stat weights if within 2 3P
244  * and with collisions from collapsed to resolved levels */
245  cs *= factor1;
246 
247  {
248  /*@-redef@*/
249  enum {DEBUG_LOC=false};
250  /*@+redef@*/
251 
252  if( DEBUG_LOC && ( nelem==ipOXYGEN ) && (cs > 0.f) && (StatesElem[ipHE_LIKE][nelem][ipHi].n == 2)
253  && ( StatesElem[ipHE_LIKE][nelem][ipLo].n <= 2 ) )
254  fprintf(ioQQQ,"%li\t%li\t%li\t-\t%li\t%li\t%li\t%.2e\t%s\n",
255  StatesElem[ipHE_LIKE][nelem][ipLo].n , StatesElem[ipHE_LIKE][nelem][ipLo].S ,
256  StatesElem[ipHE_LIKE][nelem][ipLo].l , StatesElem[ipHE_LIKE][nelem][ipHi].n ,
257  StatesElem[ipHE_LIKE][nelem][ipHi].S , StatesElem[ipHE_LIKE][nelem][ipHi].l , cs,where);
258  }
259 
260  return MAX2(cs, 1.e-10f);
261 }
262 
263 realnum AtomCSInterp(long int nelem,
264  long int ipHi,
265  long int ipLo,
266  realnum *factor1,
267  const char **where,
268  long int Collider )
269 {
270  long ipArray;
271  realnum cs;
272 
273  DEBUG_ENTRY( "AtomCSInterp()" );
274 
275  ASSERT( nelem == ipHELIUM );
276 
277  /* init values, better be positive when we exit */
278  cs = -1.f;
279 
280  /* this may be used for splitting up the collision strength within 2 3P when
281  * the lower level is withint 2 3P, and for collisions between resolved and collapsed levels.
282  * It may be set somewhere in routine, so set to negative value here as flag saying not set */
283  *factor1 = -1.f;
284 
285  /* for most of the helium iso sequence, the order of the J levels within 2 3P
286  * in increasing energy, is 0 1 2 - the exception is atomic helium itself,
287  * which is swapped, 2 1 0 */
288 
289  /* this branch is for upper and lower levels within 2p3P */
290  if( ipLo >= ipHe2p3P0 && ipHi <= ipHe2p3P2 && Collider==ipELECTRON )
291  {
292  *factor1 = 1.f;
293  /* atomic helium, use Berrington private comm cs */
294 
295  /* >>refer he1 cs Berrington, Keith, 2001, private communication - email follows
296  > Dear Gary,
297  > I could not find any literature on the He fine-structure
298  > problem (but I didn't look very hard, so there may be
299  > something somewhere). However, I did a quick R-matrix run to
300  > see what the magnitude of the collision strengths are... At
301  > 1000K, I get the effective collision strength for 2^3P J=0-1,
302  > 0-2, 1-2 as 0.8,0.7,2.7; for 10000K I get 1.2, 2.1, 6.0
303  */
304  /* indices are the same and correct, only thing special is that energies are in inverted order...was right first time. */
305  if( ipLo == ipHe2p3P0 && ipHi == ipHe2p3P1 )
306  {
307  cs = 1.2f;
308  }
309  else if( ipLo == ipHe2p3P0 && ipHi == ipHe2p3P2 )
310  {
311  cs = 2.1f;
312  }
313  else if( ipLo == ipHe2p3P1 && ipHi == ipHe2p3P2 )
314  {
315  cs = 6.0f;
316  }
317  else
318  {
319  cs = 1.0f;
320  TotalInsanity();
321  }
322 
323  *where = "Berr ";
324  /* statistical weights included */
325  }
326  /* >>chng 02 feb 25, Bray data should come first since it is the best we have. */
327  /* this branch is the Bray et al data, for n <= 5, where quantal calcs exist
328  * must exclude ipLo >= ipHe2p1P because they give no numbers for those */
329  else if( StatesElem[ipHE_LIKE][nelem][ipHi].n <= 5 &&
330  ( ipHi < iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] ) &&
331  iso.HeCS[ipHE_LIKE][nelem][ipHi][ipLo][0] >= 0.f && Collider== ipELECTRON )
332  {
333  ASSERT( *factor1 == -1.f );
334  ASSERT( ipLo < ipHi );
335  ASSERT( ipHe2p3P0 == 3 );
336 
337  /* ipLo is within 2^3P */
338  if( ipLo >= ipHe2p3P0 && ipLo <= ipHe2p3P2 )
339  {
340  /* *factor1 is ratio of statistical weights of level to term */
341 
342  /* ipHe2p3P0, ipHe2p3P1, ipHe2p3P2 have indices 3,4,and 5, but j=0,1,and 2. */
343  *factor1 = (2.f*((realnum)ipLo-3.f)+1.f) / 9.f;
344 
345  /* ipHi must be above ipHe2p3P2 since 2p3Pj->2p3Pk taken care of above */
346  ASSERT( ipHi > ipHe2p3P2 );
347  }
348  /* ipHi is within 2^3P */
349  else if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
350  {
351  ASSERT( ipLo < ipHe2p3P0 );
352 
353  *factor1 = (2.f*((realnum)ipHi-3.f)+1.f) / 9.f;
354  }
355  /* neither are within 2^3P...no splitting necessary */
356  else
357  {
358  *factor1 = 1.f;
359  }
360 
361  /* SOME OF THESE ARE NOT N-CHANGING! */
362  /* Must be careful about turning each one on or off. */
363 
364  /* this is the case where we have quantal calculations */
365  /* >>refer He1 cs Bray, I., Burgess, A., Fursa, D.V., & Tully, J.A., 2000, A&AS, 146, 481-498 */
366  /* check whether we are outside temperature array bounds,
367  * and use extreme value if we are */
368  if( phycon.alogte <= iso.CSTemp[ipHE_LIKE][0] )
369  {
370  cs = iso.HeCS[ipHE_LIKE][nelem][ipHi][ipLo][0];
371  }
372  else if( phycon.alogte >= iso.CSTemp[ipHE_LIKE][iso.nCS[ipHE_LIKE]-1] )
373  {
374  cs = iso.HeCS[ipHE_LIKE][nelem][ipHi][ipLo][iso.nCS[ipHE_LIKE]-1];
375  }
376  else
377  {
378  realnum flow;
379 
380  /* find which array element within the cs vs temp array */
381  ipArray = (long)((phycon.alogte - iso.CSTemp[ipHE_LIKE][0])/(iso.CSTemp[ipHE_LIKE][1]-iso.CSTemp[ipHE_LIKE][0]));
382  ASSERT( ipArray < iso.nCS[ipHE_LIKE] );
383  ASSERT( ipArray >= 0 );
384  /* when taking the average, this is the fraction from the lower temperature value */
385  flow = (realnum)( (phycon.alogte - iso.CSTemp[ipHE_LIKE][ipArray])/
386  (iso.CSTemp[ipHE_LIKE][ipArray+1]-iso.CSTemp[ipHE_LIKE][ipArray]));
387  ASSERT( (flow >= 0.f) && (flow <= 1.f) );
388 
389  cs = iso.HeCS[ipHE_LIKE][nelem][ipHi][ipLo][ipArray] * (1.f-flow) +
390  iso.HeCS[ipHE_LIKE][nelem][ipHi][ipLo][ipArray+1] * flow;
391  }
392 
393  *where = "Bray ";
394 
395  /* options to kill collisional excitation and/or l-mixing */
396  if( StatesElem[ipHE_LIKE][nelem][ipHi].n == StatesElem[ipHE_LIKE][nelem][ipLo].n )
397  /* iso.lgColl_l_mixing turned off with atom he-like l-mixing collisions off command */
399  else
400  {
401  /* iso.lgColl_excite turned off with atom he-like collisional excitation off command */
403  }
404 
405  ASSERT( cs >= 0.f );
406  /* statistical weights included */
407  }
408  /* this branch, n-same, l-changing collision, and not case of He with quantal data */
409  else if( (StatesElem[ipHE_LIKE][nelem][ipHi].n == StatesElem[ipHE_LIKE][nelem][ipLo].n ) &&
410  (StatesElem[ipHE_LIKE][nelem][ipHi].S == StatesElem[ipHE_LIKE][nelem][ipLo].S ) )
411  {
412  ASSERT( *factor1 == -1.f );
413  *factor1 = 1.f;
414 
415  /* ASSERT( StatesElem[ipHE_LIKE][nelem][ipHi].n >= 3 ); */
416  ASSERT( StatesElem[ipHE_LIKE][nelem][ipHi].n <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] );
417 
418  if( (StatesElem[ipHE_LIKE][nelem][ipLo].l <=2) &&
419  abs(StatesElem[ipHE_LIKE][nelem][ipHi].l - StatesElem[ipHE_LIKE][nelem][ipLo].l)== 1 )
420  {
421  /* Use the method given in
422  * >>refer He CS Seaton, M. J. 1962, Proc. Phys. Soc. 79, 1105
423  * statistical weights included */
424  cs = (realnum)CS_l_mixing_S62(ipHE_LIKE, nelem, ipLo, ipHi, (double)phycon.te, Collider);
425  }
426  else if( iso.lgCS_Vrinceanu[ipHE_LIKE] )
427  {
428  if( StatesElem[ipHE_LIKE][nelem][ipLo].l >=3 &&
429  StatesElem[ipHE_LIKE][nelem][ipHi].l >=3 )
430  {
431  /* Use the method given in
432  * >>refer He CS Vrinceanu, D. \& Flannery, M. R. 2001, PhysRevA 63, 032701
433  * statistical weights included */
435  nelem,
436  StatesElem[ipHE_LIKE][nelem][ipLo].n,
437  StatesElem[ipHE_LIKE][nelem][ipLo].l,
438  StatesElem[ipHE_LIKE][nelem][ipHi].l,
439  StatesElem[ipHE_LIKE][nelem][ipLo].S,
440  (double)phycon.te,
441  Collider );
442  }
443  else
444  {
445  cs = 0.f;
446  }
447  }
448  /* this branch, l changing by one */
449  else if( abs(StatesElem[ipHE_LIKE][nelem][ipHi].l - StatesElem[ipHE_LIKE][nelem][ipLo].l)== 1)
450  {
451  /* >>refer He cs Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */
452  /* statistical weights included */
453  cs = (realnum)CS_l_mixing_PS64( ipHE_LIKE, nelem, ipLo, ipHi, Collider);
454  }
455  else
456  {
457  /* l changes by more than 1, but same-n collision */
458  cs = 0.f;
459  }
460 
461  /* ipLo is within 2^3P */
462  if( ipLo >= ipHe2p3P0 && ipLo <= ipHe2p3P2 )
463  {
464  *factor1 = (2.f*((realnum)ipLo-3.f)+1.f) / 9.f;
465  }
466 
467  /* ipHi is within 2^3P */
468  if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
469  {
470  *factor1 = (2.f*((realnum)ipHi-3.f)+1.f) / 9.f;
471  }
472 
473  *where = "lmix ";
475  }
476 
477  /* this is an atomic n-changing collision with no quantal calculations */
478  else if( StatesElem[ipHE_LIKE][nelem][ipHi].n != StatesElem[ipHE_LIKE][nelem][ipLo].n )
479  {
480  ASSERT( *factor1 == -1.f );
481  /* this is an atomic n-changing collision with no quantal calculations */
482  /* gbar g-bar goes here */
483 
484  /* >>chng 02 jul 25, add option for fits to quantal cs data */
485  if( iso.lgCS_Vriens[ipHE_LIKE] )
486  {
487  /* >>refer He CS Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940
488  * statistical weight IS included in the routine */
489  cs = (realnum)CS_VS80(
490  ipHE_LIKE,
491  nelem,
492  ipHi,
493  ipLo,
494  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Aul,
495  phycon.te,
496  Collider );
497 
498  *factor1 = 1.f;
499  *where = "Vriens";
500  }
501  else if( iso.lgCS_None[ipHE_LIKE] )
502  {
503  cs = 0.f;
504  *factor1 = 1.f;
505  *where = "no gb";
506  }
507  else if( iso.nCS_new[ipHE_LIKE] )
508  {
509  *factor1 = 1.f;
510  /* Don't know if stat weights are included in this, but they're probably
511  * wrong anyway since they are based in part on the former (incorrect)
512  * implementation of Vriens and Smeets rates */
513 
514  /* two different fits, allowed and forbidden */
515  if( Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Aul > 1. )
516  {
517  /* permitted lines - large A */
518  double x =
519  log10(MAX2(34.7,Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN));
520 
521  if( iso.nCS_new[ipHE_LIKE] == 1 )
522  {
523  /* this is the broken power law fit, passing through both quantal
524  * calcs at high energy and asymptotically goes to VS at low energies */
525  if( x < 4.5 )
526  {
527  /* low energy fit for permitted transitions */
528  cs = (realnum)pow( 10. , -1.45*x + 6.75);
529  }
530  else
531  {
532  /* higher energy fit for permitted transitions */
533  cs = (realnum)pow( 10. , -3.33*x+15.15);
534  }
535  }
536  else if( iso.nCS_new[ipHE_LIKE] == 2 )
537  {
538  /* single parallel fit for permitted transitions, runs parallel to VS */
539  cs = (realnum)pow( 10. , -2.3*x+10.3);
540  }
541  else
542  TotalInsanity();
543  }
544  else
545  {
546  /* fit for forbidden transitions */
547  if( Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN < 25119.f )
548  {
549  cs = 0.631f;
550  }
551  else
552  {
553  cs = (realnum)pow(10.,
554  -3.*log10(Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN)+12.8);
555  }
556  }
557 
558  *where = "newgb";
559  }
560  else
561  TotalInsanity();
562 
563  /* ipLo is within 2^3P */
564  if( ipLo >= ipHe2p3P0 && ipLo <= ipHe2p3P2 )
565  {
566  *factor1 = (2.f*((realnum)ipLo-3.f)+1.f) / 9.f;
567  }
568 
569  /* ipHi is within 2^3P */
570  if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
571  {
572  *factor1 = (2.f*((realnum)ipHi-3.f)+1.f) / 9.f;
573  }
574 
575  /* options to turn off collisions */
577 
578  }
579  else
580  {
581  /* If spin changing collisions are prohibited in the l-mixing routine,
582  * they will fall here, and will have been assigned no collision strength.
583  * Assign zero for now. */
584  ASSERT( StatesElem[ipHE_LIKE][nelem][ipHi].S != StatesElem[ipHE_LIKE][nelem][ipLo].S );
585  cs = 0.f;
586  *factor1 = 1.f;
587  }
588 
589  ASSERT( cs >= 0.f );
590 
591  /*iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPCOLLIS,-1);*/
592 
593  return(cs);
594 
595 }
596 
597 /* IonCSInterp interpolate on collision strengths for element nelem */
598 realnum IonCSInterp( long nelem , long ipHi , long ipLo, realnum *factor1, const char **where, long Collider )
599 {
600  realnum cs;
601 
602  DEBUG_ENTRY( "IonCSInterp()" );
603 
604  ASSERT( nelem > ipHELIUM );
605  ASSERT( nelem < LIMELM );
606 
607  /* init values, better be positive when we exit */
608  cs = -1.f;
609 
610  /* this may be used for splitting up the collision strength for collisions between
611  * resolved and collapsed levels. It may be set somewhere in routine, so set to
612  * negative value here as flag saying not set */
613  *factor1 = -1.f;
614 
615 
616  /* >>chng 02 mar 04, the approximation here for transitions within 2p3P was not needed,
617  * because the Zhang data give these transitions. They are of the same order, but are
618  * specific to the three transitions */
619 
620  /* this branch is ground to n=2 or from n=2 to n=2, for ions only */
621  /*>>refer Helike CS Zhang, Honglin, & Sampson, Douglas H. 1987, ApJS 63, 487 */
622  if( StatesElem[ipHE_LIKE][nelem][ipHi].n==2
623  && StatesElem[ipHE_LIKE][nelem][ipLo].n<=2 && Collider==ipELECTRON)
624  {
625  *where = "Zhang";
626  *factor1 = 1.;
627 
628  /* Collisions from gound */
629  if( ipLo == ipHe1s1S )
630  {
631  switch( ipHi )
632  {
633  case 1: /* to 2tripS */
634  cs = 0.25f/POW2(nelem+1.f);
635  break;
636  case 2: /* to 2singS */
637  cs = 0.4f/POW2(nelem+1.f);
638  break;
639  case 3: /* to 2tripP0 */
640  cs = 0.15f/POW2(nelem+1.f);
641  break;
642  case 4: /* to 2tripP1 */
643  cs = 0.45f/POW2(nelem+1.f);
644  break;
645  case 5: /* to 2tripP2 */
646  cs = 0.75f/POW2(nelem+1.f);
647  break;
648  case 6: /* to 2singP */
649  cs = 1.3f/POW2(nelem+1.f);
650  break;
651  default:
652  TotalInsanity();
653  break;
654  }
656  }
657  /* collisions from 2tripS to n=2 */
658  else if( ipLo == ipHe2s3S )
659  {
660  switch( ipHi )
661  {
662  case 2: /* to 2singS */
663  cs = 2.75f/POW2(nelem+1.f);
664  break;
665  case 3: /* to 2tripP0 */
666  cs = 60.f/POW2(nelem+1.f);
667  break;
668  case 4: /* to 2tripP1 */
669  cs = 180.f/POW2(nelem+1.f);
670  break;
671  case 5: /* to 2tripP2 */
672  cs = 300.f/POW2(nelem+1.f);
673  break;
674  case 6: /* to 2singP */
675  cs = 5.8f/POW2(nelem+1.f);
676  break;
677  default:
678  TotalInsanity();
679  break;
680  }
682  }
683  /* collisions from 2singS to n=2 */
684  else if( ipLo == ipHe2s1S )
685  {
686  switch( ipHi )
687  {
688  case 3: /* to 2tripP0 */
689  cs = 0.56f/POW2(nelem+1.f);
690  break;
691  case 4: /* to 2tripP1 */
692  cs = 1.74f/POW2(nelem+1.f);
693  break;
694  case 5: /* to 2tripP2 */
695  cs = 2.81f/POW2(nelem+1.f);
696  break;
697  case 6: /* to 2singP */
698  cs = 190.f/POW2(nelem+1.f);
699  break;
700  default:
701  TotalInsanity();
702  break;
703  }
705  }
706  /* collisions from 2tripP0 to n=2 */
707  else if( ipLo == ipHe2p3P0 )
708  {
709  switch( ipHi )
710  {
711  case 4: /* to 2tripP1 */
712  cs = 8.1f/POW2(nelem+1.f);
713  break;
714  case 5: /* to 2tripP2 */
715  cs = 8.2f/POW2(nelem+1.f);
716  break;
717  case 6: /* to 2singP */
718  cs = 3.9f/POW2(nelem+1.f);
719  break;
720  default:
721  TotalInsanity();
722  break;
723  }
725  }
726  /* collisions from 2tripP1 to n=2 */
727  else if( ipLo == ipHe2p3P1 )
728  {
729  switch( ipHi )
730  {
731  case 5: /* to 2tripP2 */
732  cs = 30.f/POW2(nelem+1.f);
733  break;
734  case 6: /* to 2singP */
735  cs = 11.7f/POW2(nelem+1.f);
736  break;
737  default:
738  TotalInsanity();
739  break;
740  }
742  }
743  /* collisions from 2tripP2 to n=2 */
744  else if( ipLo == ipHe2p3P2 )
745  {
746  /* to 2singP */
747  cs = 19.5f/POW2(nelem+1.f) * (realnum)iso.lgColl_l_mixing[ipHE_LIKE];
748  }
749  else
750  TotalInsanity();
751 
752  /* statistical weights included */
753  }
754 
755  /* this branch, n-same, l-changing collisions */
756  else if( (StatesElem[ipHE_LIKE][nelem][ipHi].n == StatesElem[ipHE_LIKE][nelem][ipLo].n ) &&
757  (StatesElem[ipHE_LIKE][nelem][ipHi].S == StatesElem[ipHE_LIKE][nelem][ipLo].S ) )
758  {
759  ASSERT( *factor1 == -1.f );
760  *factor1 = 1.f;
761 
762  ASSERT( StatesElem[ipHE_LIKE][nelem][ipHi].n <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] );
763 
764  if( (StatesElem[ipHE_LIKE][nelem][ipLo].l <=2) &&
765  abs(StatesElem[ipHE_LIKE][nelem][ipHi].l - StatesElem[ipHE_LIKE][nelem][ipLo].l)== 1 )
766  {
767  /* Use the method given in
768  * >>refer He CS Seaton, M. J. 1962, Proc. Phys. Soc. 79, 1105
769  * statistical weights included */
770  cs = (realnum)CS_l_mixing_S62(ipHE_LIKE, nelem, ipLo, ipHi, (double)phycon.te, Collider);
771  }
772  else if( iso.lgCS_Vrinceanu[ipHE_LIKE] )
773  {
774  if( StatesElem[ipHE_LIKE][nelem][ipLo].l >=3 &&
775  StatesElem[ipHE_LIKE][nelem][ipHi].l >=3 )
776  {
777  /* Use the method given in
778  * >>refer He CS Vrinceanu, D. \& Flannery, M. R. 2001, PhysRevA 63, 032701
779  * statistical weights included */
781  nelem,
782  StatesElem[ipHE_LIKE][nelem][ipLo].n,
783  StatesElem[ipHE_LIKE][nelem][ipLo].l,
784  StatesElem[ipHE_LIKE][nelem][ipHi].l,
785  StatesElem[ipHE_LIKE][nelem][ipLo].S,
786  (double)phycon.te,
787  Collider );
788  }
789  else
790  {
791  cs = 0.f;
792  }
793  }
794  /* this branch, l changing by one */
795  else if( abs(StatesElem[ipHE_LIKE][nelem][ipHi].l - StatesElem[ipHE_LIKE][nelem][ipLo].l)== 1)
796  {
797  /* >>refer He cs Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */
798  /* statistical weights included */
799  cs = (realnum)CS_l_mixing_PS64( ipHE_LIKE, nelem, ipLo, ipHi, Collider);
800  }
801  else
802  {
803  /* l changes by more than 1, but same-n collision */
804  cs = 0.f;
805  }
806 
807  /* ipHi is within 2^3P, do not need to split on ipLo. */
808  if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 )
809  {
810  *factor1 = (2.f*((realnum)ipHi-3.f)+1.f) / 9.f;
811  }
812 
813  *where = "lmix ";
815  }
816 
817  /* this branch, n changing collisions for ions */
818  else if( StatesElem[ipHE_LIKE][nelem][ipHi].n != StatesElem[ipHE_LIKE][nelem][ipLo].n )
819  {
820  if( iso.lgCS_Vriens[ipHE_LIKE] )
821  {
822  /* this is the default branch */
823  /* >>refer He CS Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940
824  * statistical weight is NOT included in the routine */
825  cs = (realnum)CS_VS80(
826  ipHE_LIKE,
827  nelem,
828  ipHi,
829  ipLo,
830  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Aul,
831  phycon.te,
832  Collider );
833 
834  *factor1 = 1.f;
835  *where = "Vriens";
836  }
837  else
838  {
839  /* \todo 2 this branch and above now do the same thing. Change something. */
840  /* this branch is for testing and reached with command ATOM HE-LIKE COLLISIONS VRIENS OFF */
841  fixit(); /* use Percival and Richards here */
842 
843  cs = 0.f;
844  *where = "hydro";
845  }
846  }
847  else
848  {
849  /* what's left are deltaN=0, spin changing collisions.
850  * These have not been accounted for. */
851  /* Make sure what comes here is what we think it is. */
852  ASSERT( StatesElem[ipHE_LIKE][nelem][ipHi].n == StatesElem[ipHE_LIKE][nelem][ipLo].n );
853  ASSERT( StatesElem[ipHE_LIKE][nelem][ipHi].S != StatesElem[ipHE_LIKE][nelem][ipLo].S );
854  cs = 0.f;
855  *where = "spin ";
856  }
857 
858  ASSERT( cs >= 0.f );
859 
860  /*iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPCOLLIS,-1);*/
861 
862  return(cs);
863 }
864 
865 
866 /*CS_l_mixing_S62 - find rate for l-mixing collisions by protons, for neutrals */
867 /* The S62 stands for Seaton 1962 */
869  long ipISO,
870  long nelem /* the chemical element, 1 for He */,
871  long ipLo /* lower level, 0 for ground */,
872  long ipHi /* upper level, 0 for ground */,
873  double temp,
874  long Collider)
875 {
876  /* >>refer He l-mixing Seaton, M.J., 1962, Proc. Phys. Soc. */
877  double coll_str, deltaE;
878 
879  DEBUG_ENTRY( "CS_l_mixing_S62()" );
880 
881  if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 )
882  {
883  return 0.;
884  }
885 
886  global_temp = temp;
887  global_stat_weight = StatesElem[ipISO][nelem][ipLo].g;
888  /* >>chng 05 sep 06, RP - update energies of excited states */
889  /* global_deltaE = EVRYD*(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] -
890  iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]); */
891  global_deltaE = Transitions[ipISO][nelem][ipHi][ipLo].EnergyErg/EN1EV;
892  deltaE = global_deltaE;
893  global_I_energy_eV = EVRYD*iso.xIsoLevNIonRyd[ipISO][nelem][0];
894  global_Collider = Collider;
895 
896  ASSERT( TRANS_PROB_CONST*POW2(deltaE/WAVNRYD/EVRYD) > 0. );
897 
898  global_osc_str = Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul/
900 
901  /* This returns a thermally averaged collision strength */
902  coll_str = qg32( 0.0, 1.0, S62_Therm_ave_coll_str);
903  coll_str += qg32( 1.0, 10.0, S62_Therm_ave_coll_str);
904  ASSERT( coll_str > 0. );
905 
906  /*iso_put_error(ipISO,nelem,ipHi,ipLo,ipCollis,-1);*/
907 
908  return coll_str;
909 }
910 
911 /* The integrand for calculating the thermal average of collision strengths */
912 STATIC double S62_Therm_ave_coll_str( double proj_energy_overKT )
913 {
914 
915  double integrand, cross_section, /*Rnot,*/ osc_strength, coll_str, zOverB2;
916  double deltaE, /* betanot, */ betaone, zeta_of_betaone, /* cs1, */ cs2, temp, stat_weight;
917  double I_energy_eV, Dubya, proj_energy;
918  long i, Collider;
919 
920  DEBUG_ENTRY( "S62_Therm_ave_coll_str()" );
921 
922  /* projectile energy in eV */
923  proj_energy = proj_energy_overKT * phycon.te / EVDEGK;
924 
925  deltaE = global_deltaE;
926  osc_strength = global_osc_str;
927  temp = (double)global_temp;
928  stat_weight = global_stat_weight;
929  I_energy_eV = global_I_energy_eV;
930  Collider = global_Collider;
931 
932  /* Rnot = 1.1229*H_BAR/sqrt(ELECTRON_MASS*deltaE*EN1EV)/Bohr_radius; in units of Bohr_radius */
933 
934  proj_energy *= ColliderMass[ipELECTRON]/ColliderMass[Collider];
935  /* The deltaE here is to make sure that the collider has no less
936  * than the energy difference between the initial and final levels. */
937  proj_energy += deltaE;
938  Dubya = 0.5*(2.*proj_energy-deltaE);
939  ASSERT( Dubya > 0. );
940 
941  /* betanot = sqrt(proj_energy/I_energy_eV)*(deltaE/2./Dubya)*Rnot; */
942 
943  ASSERT( I_energy_eV > 0. );
944  ASSERT( osc_strength > 0. );
945 
946  /* from equation 33 */
947  zOverB2 = 0.5*POW2(Dubya/deltaE)*deltaE/I_energy_eV/osc_strength;
948 
949  ASSERT( zOverB2 > 0. );
950 
951  if( zOverB2 > 100. )
952  {
953  betaone = sqrt( 1./zOverB2 );
954  }
955  else if( zOverB2 < 0.54 )
956  {
957  /* Low betaone approximation */
958  betaone = (1./3.)*(log(PI)-log(zOverB2)+1.28);
959  if( betaone > 2.38 )
960  {
961  /* average with this over approximation */
962  betaone += 0.5*(log(PI)-log(zOverB2));
963  betaone *= 0.5;
964  }
965  }
966  else
967  {
968  long ip_zOverB2 = 0;
969  double zetaOVERbeta2[101] = {
970  1.030E+02,9.840E+01,9.402E+01,8.983E+01,8.583E+01,8.200E+01,7.835E+01,7.485E+01,
971  7.151E+01,6.832E+01,6.527E+01,6.236E+01,5.957E+01,5.691E+01,5.436E+01,5.193E+01,
972  4.961E+01,4.738E+01,4.526E+01,4.323E+01,4.129E+01,3.943E+01,3.766E+01,3.596E+01,
973  3.434E+01,3.279E+01,3.131E+01,2.989E+01,2.854E+01,2.724E+01,2.601E+01,2.482E+01,
974  2.369E+01,2.261E+01,2.158E+01,2.059E+01,1.964E+01,1.874E+01,1.787E+01,1.705E+01,
975  1.626E+01,1.550E+01,1.478E+01,1.409E+01,1.343E+01,1.280E+01,1.219E+01,1.162E+01,
976  1.107E+01,1.054E+01,1.004E+01,9.554E+00,9.094E+00,8.655E+00,8.234E+00,7.833E+00,
977  7.449E+00,7.082E+00,6.732E+00,6.397E+00,6.078E+00,5.772E+00,5.481E+00,5.202E+00,
978  4.937E+00,4.683E+00,4.441E+00,4.210E+00,3.989E+00,3.779E+00,3.578E+00,3.387E+00,
979  3.204E+00,3.031E+00,2.865E+00,2.707E+00,2.557E+00,2.414E+00,2.277E+00,2.148E+00,
980  2.024E+00,1.907E+00,1.795E+00,1.689E+00,1.589E+00,1.493E+00,1.402E+00,1.316E+00,
981  1.235E+00,1.157E+00,1.084E+00,1.015E+00,9.491E-01,8.870E-01,8.283E-01,7.729E-01,
982  7.206E-01,6.712E-01,6.247E-01,5.808E-01,5.396E-01};
983 
984  ASSERT( zOverB2 >= zetaOVERbeta2[100] );
985 
986  /* find beta in the table */
987  for( i=0; i< 100; ++i )
988  {
989  if( ( zOverB2 < zetaOVERbeta2[i] ) && ( zOverB2 >= zetaOVERbeta2[i+1] ) )
990  {
991  /* found the temperature - use it */
992  ip_zOverB2 = i;
993  break;
994  }
995  }
996 
997  ASSERT( (ip_zOverB2 >=0) && (ip_zOverB2 < 100) );
998 
999  betaone = (zOverB2 - zetaOVERbeta2[ip_zOverB2]) /
1000  (zetaOVERbeta2[ip_zOverB2+1] - zetaOVERbeta2[ip_zOverB2]) *
1001  (pow(10., ((double)ip_zOverB2+1.)/100. - 1.) - pow(10., ((double)ip_zOverB2)/100. - 1.)) +
1002  pow(10., (double)ip_zOverB2/100. - 1.);
1003 
1004  }
1005 
1006  zeta_of_betaone = zOverB2 * POW2(betaone);
1007 
1008  /* cs1 = betanot * bessel_k0(betanot) * bessel_k1(betanot); */
1009  cs2 = 0.5*zeta_of_betaone + betaone * bessel_k0(betaone) * bessel_k1(betaone);
1010 
1011  /* cross_section = MIN2(cs1, cs2); */
1012  cross_section = cs2;
1013 
1014  /* cross section in units of PI * a_o^2 */
1015  cross_section *= 8. * (I_energy_eV/deltaE) * osc_strength * (I_energy_eV/proj_energy);
1016 
1017  /* convert to collision strength */
1018  coll_str = cross_section * stat_weight * ( proj_energy/EVRYD );
1019 
1020  integrand = exp( -1.*(proj_energy-deltaE)*EVDEGK/temp ) * coll_str;
1021  return integrand;
1022 }
1023 
1024 /*CS_l_mixing_PS64 - find rate for l-mixing collisions by protons, for neutrals */
1026  long ipISO,
1027  long nelem /* the chemical element, 1 for He */,
1028  long ipLo /* lower level, 0 for ground */,
1029  long ipHi /* upper level, 0 for ground */,
1030  long Collider)
1031 {
1032  /* >>refer H-like l-mixing Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */
1033  /* >>refer He-like l-mixing Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */
1034  double cs, Dnl,
1035  TwoLogDebye, TwoLogRc1,
1036  factor1, factor2,
1037  bestfactor, factorpart,
1038  reduced_mass, reduced_mass_2_emass,
1039  rate, tau;
1040  const double ChargIncoming = ColliderCharge[Collider];
1041 
1042  DEBUG_ENTRY( "CS_l_mixing_PS64()" );
1043 
1044  ASSERT( ipHi > ipLo );
1045  /* In this routine, two different cutoff radii are calculated, and as per PS64,
1046  * the least of these is selected. We take the least positive result. */
1047 
1048  /* This reduced mass is in grams. */
1049  reduced_mass = dense.AtomicWeight[nelem]*ColliderMass[Collider]/
1050  (dense.AtomicWeight[nelem]+ColliderMass[Collider])*ATOMIC_MASS_UNIT;
1051  /* this mass always appears relative to the electron mass, so define it that way */
1052  reduced_mass_2_emass = reduced_mass / ELECTRON_MASS;
1053 
1054  /* This is the lifetime of ipLo. */
1055  tau = StatesElem[ipISO][nelem][ipLo].lifetime;
1056 
1057  if( ipISO == ipH_LIKE && Transitions[ipISO][nelem][ipLo][ipH1s].ipCont > 0 )
1058  {
1059  tau = (1./tau) + Transitions[ipISO][nelem][ipLo][ipH1s].Emis->Aul;
1060  tau = 1./tau;
1061  }
1062 
1063  /* equation 46 of PS64 */
1064  /* min on density added to prevent this from becoming large and negative
1065  * at very high n_e - Pengelly & Seaton did not apply this above 1e12 cm^-3 */
1066  /* This is 2 times the log of the Debye radius. */
1067  //TwoLogDebye = 1.68 + log10( phycon.te / MIN2(1e11 , dense.eden ) );
1068  /* Brocklehurst (1971, equation 3.40) has 1.181 instead of 1.68. This appears to be due to
1069  * Pengelly and Seaton neglecting one of the two factors of PI in their Equation 33 */
1070  TwoLogDebye = 1.181 + log10( phycon.te / MIN2(1e10 , dense.eden ) );
1071 
1072  /* This is w times the log of cutoff = 0.72v(tau), where tau is the lifetime of the level nl. */
1073  TwoLogRc1 = 10.95 + log10( phycon.te * tau * tau / reduced_mass_2_emass );
1074 
1075  /* this is equation 44 of PS64 */
1076  Dnl = POW2( ChargIncoming / (double)(nelem + 1. - ipISO)) * 6. * POW2( (double)StatesElem[ipISO][nelem][ipLo].n ) *
1077  ( POW2((double)StatesElem[ipISO][nelem][ipLo].n) -
1078  POW2((double)StatesElem[ipISO][nelem][ipLo].l) - StatesElem[ipISO][nelem][ipLo].l - 1);
1079 
1080  ASSERT( Dnl > 0. );
1081  ASSERT( phycon.te / Dnl / reduced_mass_2_emass > 0. );
1082 
1083  factorpart = (11.54 + log10( phycon.te / Dnl / reduced_mass_2_emass ) );
1084 
1085  if( (factor1 = factorpart + TwoLogDebye) <= 0.)
1086  factor1 = BIGDOUBLE;
1087  if( (factor2 = factorpart + TwoLogRc1) <= 0.)
1088  factor2 = BIGDOUBLE;
1089 
1090  /* Now we find the least positive result. */
1091  bestfactor = MIN2(factor1,factor2);
1092 
1093  ASSERT( bestfactor > 0. );
1094 
1095  /* If both factors are bigger than 100, toss out the result, and return SMALLFLOAT. */
1096  if( bestfactor > 100. )
1097  {
1098  return SMALLFLOAT;
1099  }
1100 
1101  /* This is the rate coefficient. Units: cm^3 s-1 */
1102  rate = 9.93e-6 * sqrt( reduced_mass_2_emass ) * Dnl / phycon.sqrte * bestfactor;
1103 
1104  /***** NB NB NB NB
1105  * Brocklehurst (1971) has a factor of electron density in the rate coefficient (equation 3.38).
1106  * This is NOT a proper rate, as can be seen in his equations 2.2 and 2.4. This differs
1107  * from the formulism given by PS64. */
1108  //rate *= dense.eden;
1109 
1110  /* this is the TOTAL rate from nl to nl+/-1. So assume we can
1111  * divide by two to get the rate nl to either nl+1 or nl-1. */
1112  if( StatesElem[ipISO][nelem][ipLo].l > 0 )
1113  rate /=2;
1114 
1115  /* convert rate to collision strength */
1116  /* NB - the term in parentheses corrects for the fact that COLL_CONST is only appropriate
1117  * for electron colliders and is off by reduced_mass_2_emass^-1.5 */
1118  cs = rate / ( COLL_CONST * pow(reduced_mass_2_emass, -1.5) ) *
1119  phycon.sqrte * StatesElem[ipISO][nelem][ipHi].g;
1120 
1121  ASSERT( cs > 0. );
1122 
1123  /*iso_put_error(ipISO,nelem,ipHi,ipLo,ipCollis,-1);*/
1124 
1125  return cs;
1126 }
1127 
1128 /*CS_l_mixing - find collision strength for l-mixing collisions for neutrals */
1129 /* The VF stands for Vrinceanu & Flannery 2001 */
1130 /* >>refer He l-mixing Vrinceanu, D. & Flannery, M. R. 2001, PhysRevA 63, 032701 */
1131 /* >>refer He l-mixing Hezel, T. P., Burkhardt, C. E., Ciocca, M., He, L-W., */
1132 /* >>refercon Leventhal, J. J. 1992, Am. J. Phys. 60, 329 */
1133 double CS_l_mixing_VF01(long int ipISO,
1134  long int nelem,
1135  long int n,
1136  long int l,
1137  long int lp,
1138  long int s,
1139  double temp,
1140  long int Collider )
1141 {
1142 
1143  double coll_str;
1144 
1145  DEBUG_ENTRY( "CS_l_mixing_VF01()" );
1146 
1147  global_ipISO = ipISO;
1148  global_z = nelem;
1149  global_n = n;
1150  global_l = l;
1151  global_l_prime = lp;
1152  global_s = s;
1153  global_temp = temp;
1154  global_Collider = Collider;
1157 
1158  /* no need to do this for h-like */
1159  if( ipISO > ipH_LIKE )
1160  {
1161  ASSERT( l != 0 );
1162  ASSERT( lp != 0 );
1163  }
1164 
1165  kTRyd = temp / TE1RYD;
1166  if( !iso.lgCS_therm_ave[ipISO] )
1167  {
1168  /* Must do some thermal averaging for densities greater
1169  * than about 10000 and less than about 1e10,
1170  * because kT gives significantly different results.
1171  * Still, do sparser integration than is done below */
1172  if( (dense.eden > 10000.) && (dense.eden < 1E10 ) )
1173  {
1174  coll_str = qg32( 0.0, 6.0, Therm_ave_coll_str_int_VF01);
1175  }
1176  else
1177  {
1178  /* Do NOT average over Maxwellian */
1179  coll_str = collision_strength_VF01( kTRyd, false );
1180  }
1181  }
1182  else
1183  {
1184  /* DO average over Maxwellian */
1185  coll_str = qg32( 0.0, 1.0, Therm_ave_coll_str_int_VF01);
1186  coll_str += qg32( 1.0, 10.0, Therm_ave_coll_str_int_VF01);
1187  }
1188 
1189  return coll_str;
1190 }
1191 
1192 /* The integrand for calculating the thermal average of collision strengths */
1193 STATIC double Therm_ave_coll_str_int_VF01( double EOverKT )
1194 {
1195  double integrand;
1196 
1197  DEBUG_ENTRY( "Therm_ave_coll_str_int_VF01()" );
1198 
1199  integrand = exp( -1.*EOverKT ) * collision_strength_VF01( EOverKT * kTRyd, false );
1200 
1201  return integrand;
1202 }
1203 
1204 STATIC double collision_strength_VF01( double velOrEner, bool lgParamIsRedVel )
1205 {
1206 
1207  double cross_section, coll_str, RMSv, aveRadius, reduced_vel, E_Proj_Ryd;
1208  double ConstantFactors, reduced_mass, CSIntegral, stat_weight;
1210  double quantum_defect1, quantum_defect2, omega_qd1, omega_qd2, omega_qd;
1211  double reduced_b_max, reduced_b_min, alphamax, alphamin, step, alpha1 /*, alpha2*/;
1212 
1213  long ipISO, nelem, n, l, lp, s, ipLo, ipHi, Collider = global_Collider;
1214 
1215  DEBUG_ENTRY( "collision_strength_VF01()" );
1216 
1217  nelem = global_z;
1218  n = global_n;
1219  ASSERT( n > 0 );
1220  l = global_l;
1221  lp = global_l_prime;
1222  s = global_s;
1223  ipISO = global_ipISO;
1224 
1225  /* >>chng 06 may 30, move these up from below. Also ipHi needs lp not l. */
1226  ipLo = iso.QuantumNumbers2Index[ipISO][nelem][n][l][s];
1227  ipHi = iso.QuantumNumbers2Index[ipISO][nelem][n][lp][s];
1228 
1229  stat_weight = StatesElem[ipISO][nelem][ipLo].g;
1230 
1231  /* no need to do this for h-like */
1232  if( ipISO > ipH_LIKE )
1233  {
1234  /* these shut up the lint, already done above */
1235  ASSERT( l > 0 );
1236  ASSERT( lp > 0 );
1237  }
1238 
1239  /* This reduced mass is in grams. */
1240  reduced_mass = dense.AtomicWeight[nelem]*ColliderMass[Collider]/
1241  (dense.AtomicWeight[nelem]+ColliderMass[Collider])*ATOMIC_MASS_UNIT;
1242  ASSERT( reduced_mass > 0. );
1243 
1244  /* this is root mean squared velocity */
1245  /* use this as projectile velocity for thermally-averaged cross-section */
1246  aveRadius = (BOHR_RADIUS_CM/((double)nelem+1.-(double)ipISO))*POW2((double)n);
1247  ASSERT( aveRadius < 1.e-4 );
1248  /* >>chng 05 jul 14, as per exchange with Ryan Porter & Peter van Hoof, avoid
1249  * roundoff error and give ability to go beyond zinc */
1250  /*ASSERT( aveRadius >= BOHR_RADIUS_CM );*/
1251  ASSERT( aveRadius > 3.9/LIMELM * BOHR_RADIUS_CM );
1252  global_an = aveRadius;
1253 
1254  /* vn = n * H_BAR/ m / r = Z * e^2 / n / H_BAR
1255  * where Z is the effective charge. */
1256  RMSv = ((double)nelem+1.-(double)ipISO)*POW2(ELEM_CHARGE_ESU)/(double)n/H_BAR;
1257  ASSERT( RMSv > 0. );
1258 
1259  ASSERT( ColliderMass[Collider] > 0. );
1260 
1261  if( lgParamIsRedVel )
1262  {
1263  /* velOrEner is a reduced velocity */
1264  reduced_vel = velOrEner;
1265  /* The proton mass is needed here because the ColliderMass array is
1266  * expressed in units of the proton mass, but here we need absolute mass. */
1267  E_Proj_Ryd = 0.5 * POW2( reduced_vel * RMSv ) * ColliderMass[Collider] *
1268  PROTON_MASS / EN1RYD;
1269  }
1270  else
1271  {
1272  /* velOrEner is a projectile energy in Rydbergs */
1273  E_Proj_Ryd = velOrEner;
1274  reduced_vel = sqrt( 2.*E_Proj_Ryd*EN1RYD/ColliderMass[Collider]/PROTON_MASS )/RMSv;
1275  }
1276 
1277  /* put limits on the reduced velocity. These limits should be more than fair. */
1278  ASSERT( reduced_vel > 1.e-10 );
1279  ASSERT( reduced_vel < 1.e10 );
1280 
1281  global_red_vel = reduced_vel;
1282 
1283  /* Factors outside integral */
1284  ConstantFactors = 4.5*PI*POW2(ColliderCharge*aveRadius/reduced_vel);
1285 
1286  /* Reduced here means in units of aveRadius: */
1287  reduced_b_min = 1.5 * ColliderCharge / reduced_vel;
1288  ASSERT( reduced_b_min > 1.e-10 );
1289  ASSERT( reduced_b_min < 1.e10 );
1290 
1291  if( ipISO == ipH_LIKE )
1292  {
1293  /* Debye radius: appears to be too large, results in 1/v^2 variation. */
1294  reduced_b_max = sqrt( BOLTZMANN*global_temp/ColliderCharge/dense.eden )
1295  / (PI2*ELEM_CHARGE_ESU)/aveRadius;
1296  }
1297  else if( ipISO == ipHE_LIKE )
1298  {
1299  quantum_defect1 = (double)n- (double)nelem/sqrt(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo]);
1300  quantum_defect2 = (double)n- (double)nelem/sqrt(iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]);
1301 
1302  /* The magnitude of each quantum defect must be between zero and one. */
1303  ASSERT( fabs(quantum_defect1) < 1.0 );
1304  ASSERT( fabs(quantum_defect1) > 0.0 );
1305  ASSERT( fabs(quantum_defect2) < 1.0 );
1306  ASSERT( fabs(quantum_defect2) > 0.0 );
1307 
1308  /* The quantum defect precession frequencies */
1309  omega_qd1 = fabs( 5.* quantum_defect1 * (1.-0.6*POW2((double)l/(double)n)) / POW3( (double)n ) / (double)l );
1310  /* >>chng 06 may 30, this needs lp not l. */
1311  omega_qd2 = fabs( 5.* quantum_defect2 * (1.-0.6*POW2((double)lp/(double)n)) / POW3( (double)n ) / (double)lp );
1312  /* Take the average for the two levels, for reciprocity. */
1313  omega_qd = 0.5*( omega_qd1 + omega_qd2 );
1314 
1315  ASSERT( omega_qd > 0. );
1316 
1317  reduced_b_max = sqrt( 1.5 * ColliderCharge * n / omega_qd )/aveRadius;
1318  }
1319  else
1320  /* rethink this before using on other iso sequences. */
1321  TotalInsanity();
1322 
1323  reduced_b_max = MAX2( reduced_b_max, reduced_b_min );
1324  ASSERT( reduced_b_max > 0. );
1325  global_bmax = reduced_b_max*aveRadius;
1326  alphamin = 1.5*ColliderCharge/(reduced_vel * reduced_b_max);
1327  alphamax = 1.5*ColliderCharge/(reduced_vel * reduced_b_min);
1328 
1329  ASSERT( alphamin > 0. );
1330  ASSERT( alphamax > 0. );
1331 
1332  alphamin = MAX2( alphamin, 1.e-30 );
1333  alphamax = MAX2( alphamax, 1.e-20 );
1334 
1335  CSIntegral = 0.;
1336 
1337  if( alphamax > alphamin )
1338  {
1339 
1340  step = (alphamax - alphamin)/5.;
1341  alpha1 = alphamin;
1342  CSIntegral += qg32( alpha1, alpha1+step, L_mix_integrand_VF01);
1343  CSIntegral += qg32( alpha1+step, alpha1+4.*step, L_mix_integrand_VF01);
1344  }
1345 
1346  /* Calculate cross section */
1347  cross_section = ConstantFactors * CSIntegral;
1348 
1349  /* convert to collision strength */
1350  coll_str = cross_section * stat_weight / PI / BOHR_RADIUS_CM / BOHR_RADIUS_CM;
1351  coll_str *= E_Proj_Ryd;
1352 
1353  coll_str = MAX2( SMALLFLOAT, coll_str);
1354 
1355  return coll_str;
1356 }
1357 
1358 STATIC double L_mix_integrand_VF01( double alpha )
1359 {
1360  double integrand, deltaPhi, b;
1361 
1362  DEBUG_ENTRY( "L_mix_integrand_VF01()" );
1363 
1364  ASSERT( alpha >= 1.e-30 );
1365  ASSERT( global_bmax > 0. );
1366  ASSERT( global_red_vel > 0. );
1367 
1368  /* >>refer He l-mixing Kazansky, A. K. & Ostrovsky, V. N. 1996, JPhysB: At. Mol. Opt. Phys. 29, 3651 */
1370  /* deltaPhi is the effective angle swept out by the projector as viewed by the target. */
1371  if( b < global_bmax )
1372  {
1373  deltaPhi = -1.*PI + 2.*asin(b/global_bmax);
1374  }
1375  else
1376  {
1377  deltaPhi = 0.;
1378  }
1379  integrand = 1./alpha/alpha/alpha;
1380  integrand *= StarkCollTransProb_VF01( global_n, global_l, global_l_prime, alpha, deltaPhi );
1381 
1382  return integrand;
1383 }
1384 
1385 STATIC double StarkCollTransProb_VF01( long n, long l, long lp, double alpha, double deltaPhi)
1386 {
1387  double probability;
1388  double cosU1, cosU2, sinU1, sinU2, cosChiOver2, sinChiOver2, cosChi, A, B;
1389 
1390  DEBUG_ENTRY( "StarkCollTransProb_VF01()" );
1391 
1392  ASSERT( alpha > 0. );
1393 
1394  /* These are defined on page 11 of VF01 */
1395  cosU1 = 2.*POW2((double)l/(double)n) - 1.;
1396  cosU2 = 2.*POW2((double)lp/(double)n) - 1.;
1397 
1398  sinU1 = sqrt( 1. - cosU1*cosU1 );
1399  sinU2 = sqrt( 1. - cosU2*cosU2 );
1400 
1401 
1402  cosChiOver2 = (1. + alpha*alpha*cos( sqrt(1.+alpha*alpha) * deltaPhi ) )/(1.+alpha*alpha);
1403  sinChiOver2 = sqrt( 1. - cosChiOver2*cosChiOver2 );
1404  cosChi = 2. * POW2( cosChiOver2 ) - 1.;
1405 
1406  if( l == 0 )
1407  {
1408  if( -1.*cosU2 - cosChi < 0. )
1409  {
1410  probability = 0.;
1411  }
1412  else
1413  {
1414  /* Here is the initial state S case */
1415  ASSERT( sinChiOver2 > 0. );
1416  ASSERT( sinChiOver2*sinChiOver2 > POW2((double)lp/(double)n) );
1417  /* This is equation 35 of VF01. There is a factor of hbar missing in the denominator of the
1418  * paper, but it's okay if you use atomic units (hbar = 1). */
1419  probability = (double)lp/(POW2((double)n)*sinChiOver2*sqrt( POW2(sinChiOver2) - POW2((double)lp/(double)n) ) );
1420  }
1421  }
1422  else
1423  {
1424  double OneMinusCosChi = 1. - cosChi;
1425 
1426  if( OneMinusCosChi == 0. )
1427  {
1428  double hold = sin( deltaPhi / 2. );
1429  /* From approximation at bottom of page 10 of VF01. */
1430  OneMinusCosChi = 8. * alpha * alpha * POW2( hold );
1431  }
1432 
1433  if( OneMinusCosChi == 0. )
1434  {
1435  probability = 0.;
1436  }
1437  else
1438  {
1439  /* Here is the general case */
1440  A = (cosU1*cosU2 - sinU1*sinU2 - cosChi)/OneMinusCosChi;
1441  B = (cosU1*cosU2 + sinU1*sinU2 - cosChi)/OneMinusCosChi;
1442 
1443  ASSERT( B > A );
1444 
1445  /* These are the three cases of equation 34. */
1446  if( B <= 0. )
1447  {
1448  probability = 0.;
1449  }
1450  else
1451  {
1452  ASSERT( POW2( sinChiOver2 ) > 0. );
1453 
1454  probability = 2.*lp/(PI* /*H_BAR* */ n*n*POW2( sinChiOver2 ));
1455 
1456  if( A < 0. )
1457  {
1458  probability *= ellpk( -A/(B-A) );
1459  probability /= sqrt( B - A );
1460  }
1461  else
1462  {
1463  probability *= ellpk( A/B );
1464  probability /= sqrt( B );
1465  }
1466  }
1467  }
1468 
1469  }
1470 
1471  return probability;
1472 
1473 }
1474 
1475 #if 0
1476 STATIC void updateHeColl(long int nelem)
1477 {
1478  long int ipLo , ipHi, i;
1479 
1480  /* collision strengths are assumed to be roughly constant for small changes in temperature
1481  * and are not recalculated as often as other data. This array stores the last temperature
1482  * at which collision strengths were evaluated for each species of the isoelectronic sequence. */
1483  static double TeUsedForCS[LIMELM]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
1484  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0};
1485 
1486  DEBUG_ENTRY("updateHeColl()");
1487  /* Calculate initial value of collision strengths, OR
1488  * Reevaluate collision strengths if the temperature has changed by more than 15%. */
1489 
1490  /* This was the case in the initial coding -- if the number/order of
1491  * colliders changes, must make sure all the auxiliary information
1492  * is correct. In the long term the code should be robust to
1493  * removing these conditions */
1494  ASSERT(ipELECTRON == 0 && ipPROTON == 1 && ipHE_PLUS == 2);
1495 
1496  if( (TeUsedForCS[nelem] == 0.) ||
1497  ( TeUsedForCS[nelem]/phycon.te > 1.15 ) ||
1498  ( TeUsedForCS[nelem]/phycon.te < 0.85 ) )
1499  {
1500  for( ipHi=ipHe2s3S; ipHi <iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
1501  {
1502  for( ipLo=ipHe1s1S; ipLo < ipHi; ipLo++ )
1503  {
1504  for (i=0;i<=ipNCOLLIDER;i++)
1505  {
1506  /* Should remove this limit when collider data is complete */
1507  if (i > ipHE_PLUS)
1508  continue;
1509 
1510  /* collsion strengths due to impact */
1511  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Coll.csi[i] =
1512  HeCSInterp( nelem , ipHi , ipLo, i );
1513 
1514  /* check for sanity */
1515  ASSERT( Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Coll.csi[i] >= 0. );
1516  }
1517  }
1518  }
1519  /* Update temperature at which collision strengths were evaluated. */
1520  TeUsedForCS[nelem] = phycon.te;
1521  }
1522 }
1523 #endif
1524 

Generated for cloudy by doxygen 1.8.4