cloudy  trunk
 All Data Structures Namespaces 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-2022 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 "dense.h"
10 #include "helike.h"
11 #include "helike_cs.h"
12 #include "hydro_vs_rates.h"
13 #include "iso.h"
14 #include "phycon.h"
15 #include "thirdparty.h"
16 #include "thirdparty_quadpack.h"
17 #include "trace.h"
18 #include "freebound.h"
19 #include "lines_service.h"
20 #include "integrate.h"
21 #include "vectorize.h"
22 #include "hydroeinsta.h"
23 
25 vector<double> CSTemp;
27 multi_arr<realnum,3> HeCS;
28 
29 STATIC realnum HeCSTableInterp( long nelem, long Collider,
30  long nHi, long lHi, long sHi, long jHi,
31  long nLo, long lLo, long sLo, long jLo );
32 
33 STATIC double CS_l_mixing_S62( double deltaE_eV, double IP_Ryd_ground, long gLo, double Aul, long nelem, long Collider, double temp );
34 
35 /* all of these are used in the calculation of Stark collision strengths
36  * following the algorithms in Vrinceanu & Flannery (2001). */
37 template<class P>
38 class my_Integrand_VF01_E;
39 
40 template<class P>
41 STATIC double collision_strength_VF01( long ipISO, double velOrEner,
42  const my_Integrand_VF01_E<P>& vf );
43 
44 /* These are masses relative to the proton mass of the electron, proton, he+, and alpha particle. */
45 static const double ColliderCharge[4] = {1.0, 1.0, 1.0, 2.0};
46 
47 inline double reduced_amu( long nelem, long Collider )
48 {
49  return dense.AtomicWeight[nelem]*colliders.list[Collider].mass_amu/
50  (dense.AtomicWeight[nelem]+colliders.list[Collider].mass_amu)*ATOMIC_MASS_UNIT;
51 }
52 
53 namespace
54 {
55  void compare_VOS12();
56  class StarkCollTransProb_VF01;
57  class StarkCollTransProb_VOS12QM;
58 }
59 
60 template<class P>
61 class my_Integrand_VF01_E
62 {
63 public:
64  long ipISO, nelem, n, l, lp, gLo;
65  double tauLo, IP_Ryd_Hi, IP_Ryd_Lo;
66  long Collider;
67  double ColliderCharge, temp;
68  double reduced_mass, aveRadius, RMSv, reduced_b_max, ColliderMass;
69 
70  my_Integrand_VF01_E( long ipISO_1, long nelem_1, long n_1, long l_1, long lp_1, long s, long gLo_1,
71  double tauLo_1, double IP_Ryd_Hi_1, double IP_Ryd_Lo_1, long Collider_1, double temp_1) :
72  ipISO(ipISO_1), nelem(nelem_1), n(n_1), l(l_1), lp(lp_1), gLo(gLo_1),
73  tauLo(tauLo_1), IP_Ryd_Hi(IP_Ryd_Hi_1), IP_Ryd_Lo(IP_Ryd_Lo_1), Collider(Collider_1), temp(temp_1)
74  {
75  ColliderCharge = ::ColliderCharge[Collider];
76  ASSERT( ColliderCharge > 0. );
77  reduced_mass = reduced_amu(nelem, Collider);
78  ASSERT( reduced_mass > 0. );
79  // Don't police max value of n here: this isn't a validity
80  // constraint *for this function*, and it may be useful to
81  // use higher values for testing
82  if( ! ( s > 0 && n > 0 ) ) // && n <= iso_sp[ipISO][nelem].n_HighestResolved_max ) )
83  {
84  fprintf( ioQQQ, "invalid parameter for my_Integrand_VF01_E\n" );
85  cdEXIT(EXIT_FAILURE);
86  }
87  /* this is root mean squared velocity */
88  /* use this as projectile velocity for thermally-averaged cross-section */
89  double vred = sqrt(2.*BOLTZMANN*temp/reduced_mass);
90  // This is the n-scaled Bohr radius
91  aveRadius = (BOHR_RADIUS_CM/((double)nelem+1.-(double)ipISO))*POW2((double)n);
92  ASSERT( aveRadius < 1.e-4 );
93  /* >>chng 05 jul 14, as per exchange with Ryan Porter & Peter van Hoof, avoid
94  * roundoff error and give ability to go beyond zinc */
95  /*ASSERT( aveRadius >= BOHR_RADIUS_CM );*/
96  ASSERT( aveRadius > 3.9/LIMELM * BOHR_RADIUS_CM );
97 
98  /* vn = n * H_BAR/ m / r = Z * e^2 / n / H_BAR
99  * where Z is the effective charge. */
100  RMSv = ((double)nelem+1.-(double)ipISO)*POW2(ELEM_CHARGE_ESU)/(double)n/H_BAR;
101  ASSERT( RMSv > 0. );
102 
103  /* From here Pengelly and Seaton cut-offs MNRAS (1964),127-165 are applied
104  * to H and He
105  */
106 
107  /* Those are variables used for cut-off calculation */
108  double meanb;
109  double deltaE_eV= abs(IP_Ryd_Hi - IP_Ryd_Lo)*EVRYD;
110  double reduced_b_maxa, reduced_b_maxb;
111  if( ipISO == ipH_LIKE )
112  {
113  /* Debye radius: appears to be too large, results in 1/v^2 variation. */
114  /* Reduced here means in units of aveRadius: */
115 
116  reduced_b_max = sqrt( BOLTZMANN*temp/ColliderCharge/dense.eden/PI)
117  / (2.*ELEM_CHARGE_ESU)/aveRadius;
118 
119  /* lifetime Radius in units of aveRadius */
120  reduced_b_maxb = 0.72*vred*tauLo/aveRadius;
121 
122  reduced_b_max = MIN2( reduced_b_max, reduced_b_maxb );
123 
124  }
125  else if( ipISO == ipHE_LIKE )
126  {
127  double omega_qd;
128 
129  /* Keep the precession criterion cut-off for Vrinceanu calculation
130  * as it seems to be large compared with PS cut-offs for all
131  * calculations. If other cut-offs applied it could not be enough for
132  * low temperatures and high densities
133  */
134  if(iso_ctrl.lgCS_Vrinceanu[ipISO] && (l != 0 && lp != 0 ))
135  {
136  double quantum_defect1 = (double)n- (double)nelem/sqrt(IP_Ryd_Lo);
137  double quantum_defect2 = (double)n- (double)nelem/sqrt(IP_Ryd_Hi);
138 
139  /* The magnitude of each quantum defect must be between zero and one. */
140  ASSERT( fabs(quantum_defect1) < 1.0 );
141  ASSERT( fabs(quantum_defect1) > 0.0 );
142  ASSERT( fabs(quantum_defect2) < 1.0 );
143  ASSERT( fabs(quantum_defect2) > 0.0 );
144 
145  /* The quantum defect precession frequencies */
146  double omega_qd1 = fabs( 5.* quantum_defect1 * (1.-0.6*POW2((double)l/(double)n)) / POW3( (double)n ) / (double)l );
147  /* >>chng 06 may 30, this needs lp not l. */
148  double omega_qd2 = fabs( 5.* quantum_defect2 * (1.-0.6*POW2((double)lp/(double)n)) / POW3( (double)n ) / (double)lp );
149  /* Take the average for the two levels, for reciprocity. */
150  omega_qd = 0.5*( omega_qd1 + omega_qd2 );
151 
152  ASSERT( omega_qd > 0. );
153 
154  reduced_b_max = sqrt( 1.5 * ColliderCharge * n / omega_qd )/aveRadius;
155  }
156  else
157  {
158  // Debye Radius, reduced in units of aveRadius
159  reduced_b_maxa = sqrt( BOLTZMANN*temp/ColliderCharge/dense.eden/PI )
160  / (2.*ELEM_CHARGE_ESU)/aveRadius;
161  // Lifetime cut-off
162  reduced_b_maxb = 0.72*vred*tauLo/aveRadius;
163 
164  // Criteria for degeneration of levels
165  if (iso_ctrl.lgCS_PSdeg[ipISO] && iso_ctrl.lgCS_VOS12QM[ipISO])
166  {
167  /*PS64 criterion for degeneration
168  * if PS_crit = 1.55 Rc = 0.72*v*tau
169  */
170  double PS_crit = deltaE_eV*tauLo/HBAReV;
171  if (PS_crit > 1.55)
172  reduced_b_maxb = 1.12*HBAReV*vred/deltaE_eV/aveRadius;
173  }
174  else if (iso_ctrl.lgCS_B72[ipISO] && iso_ctrl.lgCS_VOS12QM[ipISO])
175  {
176  /*B72 criterion cut off, Brocklehurst M., 1972, MNRAS, 157, 211
177  * beta = <L>DE/hbar/2/<W> < 4
178  */
179  if (dense.xIonDense[0][0] > 0 )
180  {
181  // mean impact parameter is the minimum of 1/(ion density) or Debye radius
182  meanb= powpq(1/dense.xIonDense[0][0],1,3);
183  meanb = MIN2(meanb,reduced_b_maxa);
184  }
185  else
186  meanb = reduced_b_maxa;
187  double beta_brock= meanb*abs(deltaE_eV)/(HBAReV*4.*PI*vred);
188  if (beta_brock >= 0.4)
189  reduced_b_maxb = 1.12*HBAReV*vred/deltaE_eV/aveRadius;
190  }
191 
192  reduced_b_max = min( reduced_b_maxb ,reduced_b_maxa);
193 
194  //we need to ensure a large enough integral limit for Vrinceanu
195  if (iso_ctrl.lgCS_Vrinceanu[ipISO])
196  reduced_b_max = max( reduced_b_maxb ,reduced_b_maxa);
197 
198  }
199  }
200  else
201  /* rethink this before using on other iso sequences. */
202  TotalInsanity();
203 
204  ColliderMass = colliders.list[Collider].mass_amu;
205  ASSERT( ColliderMass > 0. );
206  }
207 
208  double operator() (double EOverKT) const
209  {
210  /* This reduced mass is in grams. */
211  double col_str = collision_strength_VF01<P>( ipISO, EOverKT * temp / TE1RYD,
212  *this );
213  return exp( -1.*EOverKT ) * col_str;
214  }
215 };
216 
217 // Version of my_Integrand_VF01_E designed to be integrated over y = exp(-EOverkT) in the range (0,1]
218 // Changing the integrand allows the E integral to be formally extended to E=infty, and means the
219 // integrand should be more uniform
220 template<class P>
221 class my_Integrand_VF01_E_log
222 {
223 public:
224  my_Integrand_VF01_E<P> data;
225 
226  my_Integrand_VF01_E_log( long ipISO_1, long nelem_1, long n_1, long l_1, long lp_1, long s, long gLo_1,
227  double tauLo_1, double IP_Ryd_Hi_1, double IP_Ryd_Lo_1, long Collider_1, double temp_1) :
228  data(ipISO_1, nelem_1, n_1, l_1, lp_1, s, gLo_1,
229  tauLo_1, IP_Ryd_Hi_1, IP_Ryd_Lo_1, Collider_1, temp_1) {}
230 
231  double operator() (double y) const
232  {
233  double EOverKT = -log(y);
234  /* This reduced mass is in grams. */
235  double col_str = collision_strength_VF01<P>( data.ipISO, EOverKT * data.temp / TE1RYD,
236  data );
237  return col_str;
238  }
239 };
240 
241 void HeCollidSetup( void )
242 {
243  long i, i1, j, nelem, ipHi, ipLo;
244  bool lgEOL, lgHIT;
245  FILE *ioDATA;
246 
247  static const int chLine_LENGTH = 1000;
248  char chLine[chLine_LENGTH];
249 
250  DEBUG_ENTRY( "HeCollidSetup()" );
251 
252  if (0)
253  {
254  compare_VOS12();
255  cdEXIT(EXIT_FAILURE);
256  }
257 
258  /* get the collision strength data for the He 1 lines */
259  if( trace.lgTrace )
260  fprintf( ioQQQ," HeCollidSetup opening he1_cs.dat:");
261 
262  ioDATA = open_data( "he1_cs.dat", "r" );
263 
264  /* check that magic number is ok */
265  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
266  {
267  fprintf( ioQQQ, " HeCollidSetup could not read first line of he1_cs.dat.\n");
268  cdEXIT(EXIT_FAILURE);
269  }
270  i = 1;
271  i1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
272  /*i2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
273  i3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);*/
274 
275  /* the following is to check the numbers that appear at the start of he1_cs.dat */
276  if( i1 !=COLLISMAGIC )
277  {
278  fprintf( ioQQQ,
279  " HeCollidSetup: the version of he1_cs.dat is not the current version.\n" );
280  fprintf( ioQQQ,
281  " HeCollidSetup: I expected to find the number %i and got %li instead.\n" ,
282  COLLISMAGIC, i1 );
283  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
284  cdEXIT(EXIT_FAILURE);
285  }
286 
287  /* get the array of temperatures */
288  lgHIT = false;
289  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
290  {
291  /* only look at lines without '#' in first col */
292  if( chLine[0] == '#')
293  continue;
294 
295  lgHIT = true;
296  lgEOL = false;
297  char *chTemp = strtok(chLine," \t\n");
298  while( chTemp != NULL )
299  {
300  CSTemp.push_back( atof(chTemp) );
301  chTemp = strtok(NULL," \t\n");
302  }
303  break;
304  }
305  if( !lgHIT )
306  {
307  fprintf( ioQQQ, " HeCollidSetup could not find line in CS temperatures in he1_cs.dat.\n");
308  cdEXIT(EXIT_FAILURE);
309  }
310  ASSERT( CSTemp.size() == 9U );
311 
312  /* create space for array of CS values, if not already done */
313  {
314  long nelem = ipHELIUM;
315  long numLevs = iso_sp[ipHE_LIKE][nelem].numLevels_max - iso_sp[ipHE_LIKE][nelem].nCollapsed_max;
316  HeCS.reserve( numLevs );
317  for( long ipHi=1; ipHi < numLevs; ++ipHi )
318  {
319  HeCS.reserve( ipHi, ipHi );
320  for( long ipLo=0; ipLo<ipHi; ++ipLo )
321  HeCS.reserve( ipHi, ipLo, CSTemp.size() );
322  }
323  HeCS.alloc();
324  for( long ipHi=1; ipHi < numLevs; ++ipHi )
325  {
326  for( long ipLo=0; ipLo<ipHi; ++ipLo )
327  {
328  for( unsigned i=0; i< CSTemp.size(); ++i )
329  {
330  HeCS[ipHi][ipLo][i] = -1.f;
331  }
332  }
333  }
334  }
335 
336  /* now read in the CS data */
337  nelem = ipHELIUM;
338  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
339  {
340  char *chTemp;
341  /* only look at lines without '#' in first col */
342  if( (chLine[0] == ' ') || (chLine[0]=='\n') )
343  break;
344  if( chLine[0] != '#')
345  {
346  /* get lower and upper level index */
347  j = 1;
348  ipLo = (long)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
349  ipHi = (long)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
350  ASSERT( ipLo < ipHi );
351  if( ipHi >= iso_sp[ipHE_LIKE][nelem].numLevels_max - iso_sp[ipHE_LIKE][nelem].nCollapsed_max )
352  continue;
353  else
354  {
355  chTemp = chLine;
356  /* skip over 4 tabs to start of cs data */
357  for( long i=0; i<3; ++i )
358  {
359  if( (chTemp = strchr_s( chTemp, '\t' )) == NULL )
360  {
361  fprintf( ioQQQ, " HeCollidSetup could not init cs\n" );
362  cdEXIT(EXIT_FAILURE);
363  }
364  ++chTemp;
365  }
366 
367  for( unsigned i=0; i< CSTemp.size(); ++i )
368  {
369  double a;
370  if( (chTemp = strchr_s( chTemp, '\t' )) == NULL )
371  {
372  fprintf( ioQQQ, " HeCollidSetup could not scan cs, current indices: %li %li\n", ipHi, ipLo );
373  cdEXIT(EXIT_FAILURE);
374  }
375  ++chTemp;
376  sscanf( chTemp , "%le" , &a );
377  HeCS[ipHi][ipLo][i] = (realnum)a;
378  }
379  }
380  }
381  }
382 
383  /* close the data file */
384  fclose( ioDATA );
385 
386  return;
387 }
388 
389 /* Choose either AtomCSInterp or IonCSInterp */
390 realnum HeCSInterp(long nelem,
391  long ipHi,
392  long ipLo,
393  long Collider )
394 {
395  DEBUG_ENTRY( "HeCSInterp()" );
396 
397  ASSERT( nelem >= ipHELIUM );
398  ASSERT( nelem < LIMELM );
399 
400  long nHi = iso_sp[ipHE_LIKE][nelem].st[ipHi].n();
401  long lHi = iso_sp[ipHE_LIKE][nelem].st[ipHi].l();
402  long sHi = iso_sp[ipHE_LIKE][nelem].st[ipHi].S();
403  long jHi = iso_sp[ipHE_LIKE][nelem].st[ipHi].j();
404  long gHi = iso_sp[ipHE_LIKE][nelem].st[ipHi].g();
405  double IP_Ryd_Hi = iso_sp[ipHE_LIKE][nelem].fb[ipHi].xIsoLevNIonRyd;
406  long nLo = iso_sp[ipHE_LIKE][nelem].st[ipLo].n();
407  long lLo = iso_sp[ipHE_LIKE][nelem].st[ipLo].l();
408  long sLo = iso_sp[ipHE_LIKE][nelem].st[ipLo].S();
409  long jLo = iso_sp[ipHE_LIKE][nelem].st[ipLo].j();
410  long gLo = iso_sp[ipHE_LIKE][nelem].st[ipLo].g();
411  double IP_Ryd_Lo = iso_sp[ipHE_LIKE][nelem].fb[ipLo].xIsoLevNIonRyd;
412  double Aul = iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul();
413  // collisions are from high to low level, then initial level lifetime is from higher level
414  double tauLo = iso_sp[ipHE_LIKE][nelem].st[ipHi].lifetime();
415  double EnerWN = iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).EnergyWN();
416  double EnerErg = iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).EnergyErg();
417 
418  if( !iso_ctrl.lgColl_excite[ipHE_LIKE] ||
419  ( nHi==nLo && !iso_ctrl.lgColl_l_mixing[ipHE_LIKE] ) )
420  {
421  return 0.f;
422  }
423  const char *where=" ";
424  realnum cs = GetHelikeCollisionStrength( nelem, Collider,
425  nHi, lHi, sHi, jHi, gHi, IP_Ryd_Hi,
426  nLo, lLo, sLo, jLo, gLo, IP_Ryd_Lo,
427  Aul, tauLo, EnerWN, EnerErg, &where );
428 
429  ASSERT( cs >= 0.f );
430 
431  return cs;
432 }
433 
434 realnum GetHelikeCollisionStrength( long nelem, long Collider,
435  long nHi, long lHi, long sHi, long jHi, long gHi, double IP_Ryd_Hi,
436  long nLo, long lLo, long sLo, long jLo, long gLo, double IP_Ryd_Lo,
437  double Aul, double tauLo, double EnerWN, double EnerErg, const char **where )
438 {
439  DEBUG_ENTRY( "GetHelikeCollisionStrength()" );
440 
441  /* This variable is for diagnostic purposes:
442  * a string used in the output to designate where each cs comes from. */
443  *where = " ";
444 
445  /* init values, better be positive when we exit */
446  realnum cs = -1.f;
447 
448  /* this may be used for splitting up the collision strength within 2^3P */
449  realnum factor1 = 1.f;
450 
451  /* Energy difference in eV */
452  double deltaE_eV = EnerErg/EN1EV;
453 
454  /* lower level is within 2^3P */
455  if( nLo==2 && lLo==1 && sLo==3 )
456  {
457  factor1 *= (2.f*jLo+1.f) / 9.f;
458  }
459 
460  /* upper level is within 2^3P */
461  if( nHi==2 && lHi==1 && sHi==3 )
462  {
463  factor1 *= (2.f*jHi+1.f) / 9.f;
464  }
465 
466  /* for most of the helium iso sequence, the order of the J levels within 2 3P
467  * in increasing energy, is 0 1 2 - the exception is atomic helium itself,
468  * which is swapped, 2 1 0 */
469 
470  /* this branch is for tabulated data, mostly from Bray et al 2000. */
471  if( (cs = HeCSTableInterp( nelem, Collider, nHi, lHi, sHi, jHi, nLo, lLo, sLo, jLo )) >= 0.f )
472  {
473  // These are already j-resolved, so return this to unity.
474  if( nLo==2 && lLo==1 && sLo==3 && nHi==2 && lHi==1 && sHi==3 )
475  {
476  factor1 = 1.f;
477  }
478 
479  *where = "Table ";
480 
481  ASSERT( cs >= 0.f );
482  /* statistical weights included */
483  }
484  /* this branch is ground to n=2 or from n=2 to n=2, for ions only */
485  /*>>refer Helike CS Zhang, Honglin, & Sampson, Douglas H. 1987, ApJS 63, 487 */
486  else if( nelem!=ipHELIUM && nHi==2 && nLo<=2 && Collider==ipELECTRON)
487  {
488  *where = "Zhang ";
489  factor1 = 1.;
490 
491  /* Collisions from gound */
492  if( nLo == 1 )
493  {
494  if( lHi==0 && sHi==3 ) // to 2tripS
495  cs = 0.25f/POW2(nelem+1.f);
496  else if( lHi==0 && sHi==1 ) // to 2singS
497  cs = 0.4f/POW2(nelem+1.f);
498  else if( lHi==1 && sHi==3 && jHi==0 ) // to 2tripP0
499  cs = 0.15f/POW2(nelem+1.f);
500  else if( lHi==1 && sHi==3 && jHi==1 ) // to 2tripP1
501  cs = 0.45f/POW2(nelem+1.f);
502  else if( lHi==1 && sHi==3 && jHi==2 ) // to 2tripP2
503  cs = 0.75f/POW2(nelem+1.f);
504  else if( lHi==1 && sHi==1 ) // to 2singP
505  cs = 1.3f/POW2(nelem+1.f);
506  else
507  TotalInsanity();
508  }
509  /* collisions from 2tripS to n=2 */
510  else if( nLo==2 && lLo==0 && sLo==3 )
511  {
512  if( lHi==0 && sHi==1 ) // to 2singS
513  cs = 2.75f/POW2(nelem+1.f);
514  else if( lHi==1 && sHi==3 && jHi==0 ) // to 2tripP0
515  cs = 60.f/POW2(nelem+1.f);
516  else if( lHi==1 && sHi==3 && jHi==1 ) // to 2tripP1
517  cs = 180.f/POW2(nelem+1.f);
518  else if( lHi==1 && sHi==3 && jHi==2 ) // to 2tripP2
519  cs = 300.f/POW2(nelem+1.f);
520  else if( lHi==1 && sHi==1 ) // to 2singP
521  cs = 5.8f/POW2(nelem+1.f);
522  else
523  TotalInsanity();
524  }
525  /* collisions from 2singS to n=2 */
526  else if( nLo==2 && lLo==0 && sLo==1 )
527  {
528  if( lHi==1 && sHi==3 && jHi==0 ) // to 2tripP0
529  cs = 0.56f/POW2(nelem+1.f);
530  else if( lHi==1 && sHi==3 && jHi==1 ) // to 2tripP1
531  cs = 1.74f/POW2(nelem+1.f);
532  else if( lHi==1 && sHi==3 && jHi==2 ) // to 2tripP2
533  cs = 2.81f/POW2(nelem+1.f);
534  else if( lHi==1 && sHi==1 ) // to 2singP
535  cs = 190.f/POW2(nelem+1.f);
536  else
537  TotalInsanity();
538  }
539  /* collisions from 2tripP0 to n=2 */
540  else if( nLo==2 && lLo==1 && sLo==3 && jLo==0 )
541  {
542  if( lHi==1 && sHi==3 && jHi==1 ) // to 2tripP1
543  cs = 8.1f/POW2(nelem+1.f);
544  else if( lHi==1 && sHi==3 && jHi==2 ) // to 2tripP2
545  cs = 8.2f/POW2(nelem+1.f);
546  else if( lHi==1 && sHi==1 ) // to 2singP
547  cs = 3.9f/POW2(nelem+1.f);
548  else
549  TotalInsanity();
550  }
551  /* collisions from 2tripP1 to n=2 */
552  else if( nLo==2 && lLo==1 && sLo==3 && jLo==1 )
553  {
554  if( lHi==1 && sHi==3 && jHi==2 ) // to 2tripP2
555  cs = 30.f/POW2(nelem+1.f);
556  else if( lHi==1 && sHi==1 ) // to 2singP
557  cs = 11.7f/POW2(nelem+1.f);
558  else
559  TotalInsanity();
560  }
561  /* collisions from 2tripP2 to n=2 */
562  else if( nLo==2 && lLo==1 && sLo==3 && jLo==2 )
563  {
564  ASSERT( lHi==1 && sHi==1 );
565  /* to 2singP */
566  cs = 19.5f/POW2(nelem+1.f) * (realnum)iso_ctrl.lgColl_l_mixing[ipHE_LIKE];
567  }
568  else
569  TotalInsanity();
570 
571  /* statistical weights included */
572  }
573 
574  /* this branch, n-same, l-changing collisions, with no data */
575  else if( (nHi==nLo) && (sHi==sLo) )
576  {
577  ASSERT( nHi <= iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max );
578  if( iso_ctrl.lgCS_Seaton[ipHE_LIKE] && Collider == ipELECTRON )
579  {
580  /* Use the method given in
581  * >>refer He CS Seaton, M. J. 1962, Proc. Phys. Soc. 79, 1105
582  * statistical weights included. This must be default for non-degenerate case and l<3 */
583  /* FG 01/10/2020: Seaton does not deal well with protons, and so it is degraded to deal
584  * with only electrons. The default for l<3 is now PSM.
585  */
586  double IP_Ryd_ground = iso_sp[ipHE_LIKE][nelem].fb[0].xIsoLevNIonRyd;
587  cs = (realnum)CS_l_mixing_S62( deltaE_eV, IP_Ryd_ground, gHi, Aul, nelem, Collider, (double)phycon.te );
588  *where = "S62 ";
589  }
590  else if( iso_ctrl.lgCS_Vrinceanu[ipHE_LIKE] )
591  {
592 
593  if( (lLo>=3 && lHi>=3))
594  {
595  if (lHi != lLo )
596  {
597  /* Use the method given in
598  * >>refer He CS Vrinceanu, D. \& Flannery, M. R. 2001, PhysRevA 63, 032701
599  * statistical weights included.
600  * All multipoles included */
601  cs = (realnum)CS_l_mixing_VF01( ipHE_LIKE,
602  nelem,
603  nLo,
604  lHi,
605  lLo,
606  sLo,
607  gHi,
608  tauLo,
609  IP_Ryd_Hi,
610  IP_Ryd_Lo,
611  (double)phycon.te,
612  Collider );
613  *where = "VF01 ";
614 
615  }
616  else
617  {
618  cs = 0.f;
619  *where = "none ";
620  }
621 
622  }
623  else
624  {
625  cs = 0.f;
626  *where = "none ";
627  }
628  }
629  else if( iso_ctrl.lgCS_VOS12[ipHE_LIKE] )
630  {
631  if(lLo>=3 && lHi>=3)
632  {
633  /* That avoids problems when dl=0 in case of j changing
634  * collisions as 2^3P_0 -> 2^3P_1 transition
635  * (non degenerate case)
636  */
637  if ( lLo != lHi )
638  {
639  /* Use the method given in
640  * >>refer He CS Vrinceau etal ApJ 747, 56 (2012) equation (7) semiclassical treatment
641  * statistical weights included */
642  cs = (realnum)CS_l_mixing_VOS12(
643  nHi,
644  lHi,
645  lLo,
646  nelem,
647  gHi,
648  nelem+1.
649  -ipHE_LIKE,
650  Collider,
651  phycon.sqrte);
652  *where = "VOS12 ";
653  }
654  else
655  {
656  cs = 0.f;
657  *where = "none ";
658  }
659  }
660  else
661  {
662  cs = 0.f;
663  *where = "none ";
664  }
665  }
666  else if( iso_ctrl.lgCS_VOS12QM[ipHE_LIKE] )
667  {
668 
669  if ( lLo != lHi )
670  {
671  /* Use the method given in
672  * >>refer He CS Vrinceau et al ApJ 747, 56 (2012) eq. (2) quantal treatment
673  * statistical weights included */
674  cs = (realnum)CS_l_mixing_VOS12QM( ipHE_LIKE,
675  nelem,
676  nLo,
677  lHi,
678  lLo,
679  sLo,
680  gHi,
681  tauLo,
682  IP_Ryd_Hi,
683  IP_Ryd_Lo,
684  (double)phycon.te,
685  Collider );
686  *where = "VOS12Q";
687  }
688  else
689  {
690  cs = 0.f;
691  *where = "none ";
692  }
693  if (0)
694  {
695  /* This is for diagnostic, in order to compare different methods */
696  long nLo1 = nLo, nHi1 = nHi;
697  long lLo1 = lLo;
698  long lHi1 = lHi;
699  long ipISO1 = ipHE_LIKE;
700  double cs1 = (realnum)CS_l_mixing_VF01( ipISO1,
701  nelem,
702  nLo1,
703  lLo1,
704  lHi1,
705  sLo,
706  gLo,
707  tauLo,
708  IP_Ryd_Hi,
709  IP_Ryd_Lo,
710  (double)phycon.te,
711  Collider );
712  double cs2 = (realnum)CS_l_mixing_VF01( ipISO1,
713  nelem,
714  nHi1,
715  lHi1,
716  lLo1,
717  sHi,
718  gHi,
719  tauLo,
720  IP_Ryd_Lo,
721  IP_Ryd_Hi,
722  (double)phycon.te,
723  Collider );
724  double cs3 = (realnum)CS_l_mixing_VOS12QM( ipISO1,
725  nelem,
726  nHi1,
727  lHi1,
728  lLo1,
729  sHi,
730  gHi,
731  tauLo,
732  IP_Ryd_Lo,
733  IP_Ryd_Hi,
734  (double)phycon.te,
735  Collider );
736  fprintf(ioQQQ,"Check VF01 %ld %ld %ld %ld %ld %ld %g: %g (%g) VOS12 %g (%g) VOS12QM %g PS64 %g\n",
737  Collider,nLo1,lLo1,lHi1,sLo,nelem,phycon.te,
738  cs1,cs2,
739  CS_l_mixing_VOS12(nLo1,lLo1,lHi1,nelem,gLo,nelem+1.-ipHE_LIKE,Collider,phycon.sqrte),
740  CS_l_mixing_VOS12(nHi1,lHi1,lLo1,nelem,gHi,nelem+1.-ipHE_LIKE,Collider,phycon.sqrte),
741  cs3,
742  CS_l_mixing_PS64(nelem,tauLo,nelem+1.-ipISO1,nLo1,lLo1,gHi,lHi1,deltaE_eV,Collider));
743  double rate_t1, rate_t2, rate_t;
744  double oHi=1./(double)gHi;
745  double ogLo=1./(double)gLo;
746  double reduced_mass_collider_system = dense.AtomicWeight[nelem]*colliders.list[Collider].mass_amu/
747  (dense.AtomicWeight[nelem]+colliders.list[Collider].mass_amu)*ATOMIC_MASS_UNIT;
748  double ratef = powpq(ELECTRON_MASS/reduced_mass_collider_system,3,2) * COLL_CONST/phycon.sqrte;
749 
750  rate_t1 = cs1*ratef*oHi;
751  rate_t2 = cs2*ratef*ogLo;
752  rate_t = rate_t1+rate_t2;
753  fprintf(ioQQQ,"Rates for H %ld %ld %ld %ld %ld %ld %ld %ld %g: %g %g \n",
754  nLo1,lLo1,lHi1,sLo,sHi,jLo,jHi,nelem,phycon.te,rate_t, dense.eden);
755  cdEXIT(EXIT_FAILURE);
756  }
757  }
758  /* this branch, l changing by one */
759  else if( abs(lHi-lLo)==1 && iso_ctrl.lgCS_PS64[ipHE_LIKE])
760  {
761  /* >>refer He cs Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */
762  /* P & S 1964 only consider dipole l-changing */
763  /* statistical weights included */
764  /*FG 01/10/2020: PS64 and PSM deal with l<3 by default */
765  if (iso_ctrl.lgCS_PSClassic[ipHE_LIKE])
766  {
767  /* Classical Pengelly and Seaton */
768  cs = (realnum)CS_l_mixing_PS64(
769  nelem,
770  tauLo,
771  nelem+1.-ipHE_LIKE,
772  nLo,
773  lHi,
774  gHi,
775  lLo,
776  deltaE_eV,
777  Collider);
778  *where = "PS64 ";
779  }
780  else
781  {
782  /* PS-M: Modified PS method
783  * Refer to F. Guzman et al. MNRAS (2016) 464, 312
784  */
785  cs = CS_l_mixing_PS64_expI(
786  nelem,
787  tauLo,
788  nelem+1.-ipHE_LIKE,
789  nHi,
790  lHi,
791  gHi,
792  lLo,
793  deltaE_eV,
794  Collider);
795  *where = "PSM ";
796  }
797  }
798  else
799  {
800  /* l changes by more than 1, but same-n collision */
801  cs = 0.f;
802  *where = "none ";
803  }
804  }
805 
806  /* n-changing collision with no quantal calculations */
807  else if( nHi != nLo )
808  {
809  if( iso_ctrl.lgCS_Vriens[ipHE_LIKE] )
810  {
811  /* >>refer He CS Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940
812  * statistical weight IS included in the routine */
813  double Aul_j = HydroEinstA(nHi,nLo);
814  cs = (realnum)CS_VS80( nHi, gHi, IP_Ryd_Hi, nLo, gLo, IP_Ryd_Lo, Aul_j, nelem, Collider, phycon.te );
815  *where = "Vriens";
816  if( nelem != ipHELIUM )
817  {
818  fixit( "This should not be set here!" );
819  factor1 = 1.f;
820  }
821  }
822  else if( iso_ctrl.lgCS_None[ipHE_LIKE] )
823  {
824  cs = 0.f;
825  *where = "no gb ";
826  }
827  else if( iso_ctrl.nCS_new[ipHE_LIKE] && nelem==ipHELIUM )
828  {
829  /* Don't know if stat weights are included in this, but they're probably
830  * wrong anyway since they are based in part on the former (incorrect)
831  * implementation of Vriens and Smeets rates */
832 
833  /* two different fits, allowed and forbidden */
834  if( Aul > 1. )
835  {
836  /* permitted lines - large A */
837  double x = log10(MAX2(34.7,EnerWN));
838 
839  if( iso_ctrl.nCS_new[ipHE_LIKE] == 1 )
840  {
841  /* this is the broken power law fit, passing through both quantal
842  * calcs at high energy and asymptotically goes to VS at low energies */
843  if( x < 4.5 )
844  {
845  /* low energy fit for permitted transitions */
846  cs = (realnum)exp10( -1.45*x + 6.75);
847  }
848  else
849  {
850  /* higher energy fit for permitted transitions */
851  cs = (realnum)exp10( -3.33*x+15.15);
852  }
853  }
854  else if( iso_ctrl.nCS_new[ipHE_LIKE] == 2 )
855  {
856  /* single parallel fit for permitted transitions, runs parallel to VS */
857  cs = (realnum)exp10( -2.3*x+10.3);
858  }
859  else
860  TotalInsanity();
861  }
862  else
863  {
864  /* fit for forbidden transitions */
865  if( EnerWN < 25119.f )
866  {
867  cs = 0.631f;
868  }
869  else
870  {
871  cs = (realnum)exp10( -3.*log10(EnerWN)+12.8);
872  }
873  }
874 
875  *where = "newgb ";
876  }
877  else if( nelem != ipHELIUM )
878  {
879  /* \todo 2 this branch and above now do the same thing. Change something. */
880  /* this branch is for testing and reached with command SPECIES HE-LIKE COLLISIONS VRIENS OFF */
881  fixit( "Use Percival and Richards here." );
882 
883  cs = 0.f;
884  *where = "hydro ";
885  }
886  else
887  TotalInsanity();
888  }
889  else
890  {
891  /* what's left are deltaN=0, spin changing collisions.
892  * These have not been accounted for. */
893  /* Make sure what comes here is what we think it is. */
894  ASSERT( nHi == nLo );
895  ASSERT( sHi != sLo );
896  cs = 0.f;
897  *where = "spin ";
898  }
899 
900  /* take factor into account, usually 1, ratio of stat weights if within 2 3P
901  * and with collisions from collapsed to resolved levels */
902  cs *= factor1;
903 
904  {
905  /*@-redef@*/
906  enum {DEBUG_LOC=false};
907  /*@+redef@*/
908 
909  if( DEBUG_LOC && /* ( nelem==ipOXYGEN ) &&*/ (cs > 0.f) ) //&& (iso_sp[ipHE_LIKE][nelem].st[ipHi].n() == 2)
910  //&& ( iso_sp[ipHE_LIKE][nelem].st[ipLo].n() <= 2 ) )
911  fprintf(ioQQQ,"DEBUGGG HeCSInterp %li\t%li\t%li\t%li\t-\t%li\t%li\t%li\t%e\t%s\n",
912  nelem,
913  nLo, sLo, lLo,
914  nHi, sHi, lHi,
915  cs, *where );
916  }
917 
918  ASSERT( cs >= 0.f );
919 
920  return cs;
921 }
922 
923 STATIC realnum HeCSTableInterp( long nelem, long Collider,
924  long nHi, long lHi, long sHi, long jHi,
925  long nLo, long lLo, long sLo, long jLo )
926 {
927  DEBUG_ENTRY( "HeCSTableInterp()" );
928 
929  if( nelem!=ipHELIUM )
930  return -1.f;
931  else if( Collider!= ipELECTRON )
932  return -1.f;
933  else if( nHi > 5 )
934  return -1.f;
935  else if( nHi > iso_sp[ipHE_LIKE][nelem].n_HighestResolved_max )
936  {
937  fixit( "HeCS allocation must be changed in order to remove this and do collapsed levels here." );
938  return -1.f;
939  }
940  // Cannot do collapsed levels with unspecified l here.
941  else if( lLo < 0 || lHi < 0 )
942  return -1.f;
943 
944  long ipLo = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[nLo][lLo][sLo];
945  long ipHi = iso_sp[ipHE_LIKE][nelem].QuantumNumbers2Index[nHi][lHi][sHi];
946 
947  ASSERT( ipHe2p3P0 == 3 );
948  // make sure 2^3Pj terms are fully specified
949  if( nLo==2 && lLo==1 && sLo==3 )
950  {
951  ASSERT( jLo>=0 && jLo<=2 );
952  ipLo -= (2 - jLo);
953  }
954  if( nHi==2 && lHi==1 && sHi==3 )
955  {
956  ASSERT( jHi>=0 && jHi<=2 );
957  ipHi -= (2 - jHi);
958  }
959 
960  // Ensure indices make sense.
961  ASSERT( ipLo < ipHi );
962 
963  if( HeCS[ipHi][ipLo][0] < 0.f )
964  return -1.f;
965 
966  realnum cs = -1.f;
967  /* this is the case where we have quantal calculations */
968  /* >>refer He1 cs Bray, I., Burgess, A., Fursa, D.V., & Tully, J.A., 2000, A&AS, 146, 481-498 */
969  /* check whether we are outside temperature array bounds,
970  * and use extreme value if we are */
971  if( phycon.alogte <= CSTemp[0] )
972  {
973  cs = HeCS[ipHi][ipLo][0];
974  }
975  else if( phycon.alogte >= CSTemp.back() )
976  {
977  cs = HeCS[ipHi][ipLo][CSTemp.size()-1];
978  }
979  else
980  {
981  /* find which array element within the cs vs temp array */
982  long ipArray = (long)((phycon.alogte - CSTemp[0])/(CSTemp[1]-CSTemp[0]));
983  ASSERT( (unsigned)ipArray < CSTemp.size() );
984  ASSERT( ipArray >= 0 );
985  /* when taking the average, this is the fraction from the lower temperature value */
986  double flow = ( (phycon.alogte - CSTemp[ipArray])/
987  (CSTemp[ipArray+1]-CSTemp[ipArray]));
988  ASSERT( (flow >= 0.) && (flow <= 1.) );
989 
990  cs = HeCS[ipHi][ipLo][ipArray] * (1.-flow) +
991  HeCS[ipHi][ipLo][ipArray+1] * flow;
992  }
993 
994  ASSERT( cs >= 0.f );
995  return cs;
996 }
997 
998 class my_Integrand_S62
999 {
1000  double deltaE, osc_strength, temp, I_energy_eV;
1001 public:
1002  my_Integrand_S62(double deltaE, double osc_strength, double temp, double I_energy_eV) :
1003  deltaE(deltaE), osc_strength(osc_strength), temp(temp), I_energy_eV(I_energy_eV)
1004  {}
1005  void operator() (const double proj_energy_overKT[], double res[], long n) const;
1006 };
1007 
1008 /*CS_l_mixing_S62 - find rate for l-mixing collisions by protons, for neutrals */
1009 /* The S62 stands for Seaton 1962 */
1010 STATIC double CS_l_mixing_S62( double deltaE_eV, double IP_Ryd_ground, long gLo, double Aul, long nelem, long Collider, double temp )
1011 {
1012  /* >>refer He l-mixing Seaton, M.J., 1962, Proc. Phys. Soc. */
1013  DEBUG_ENTRY( "CS_l_mixing_S62()" );
1014 
1015  if( Aul <= iso_ctrl.SmallA )
1016  {
1017  return 0.;
1018  }
1019 
1020  ASSERT( TRANS_PROB_CONST*POW2(deltaE_eV/WAVNRYD/EVRYD) > 0. );
1021 
1022  double osc_strength = Aul / (TRANS_PROB_CONST*POW2(deltaE_eV/WAVNRYD/EVRYD));
1023  double I_energy_eV = EVRYD * IP_Ryd_ground;
1024  double reduced_mass(reduced_amu(nelem,Collider));
1025 
1026  // Integral is over [0,infty) in principle
1027  // Dimensionless parameters for integral can then be taken as
1028  // 1. kT/I_energy_eV; 2. osc_strength; 3. deltaE/I_energy_eV
1029  // Tabulation may be possible? NB 1. is independent of atomic physics
1030  // parameters specific to the system. At present, the integral is independent of
1031  // the collider type, but seems likely that this will need to be corrected.
1032 
1033  // Major cost is very accurate Bessel function evaluations, which
1034  // doesn't seem consistent with accuracy of interpolation in
1035  // deriving zeta
1036  my_Integrand_S62 func(deltaE_eV, osc_strength, temp, I_energy_eV);
1037 
1038  double energy_factor = phycon.te/EVDEGK * colliders.list[ipELECTRON].mass_amu/colliders.list[Collider].mass_amu;
1039  double elow = 0., emid = 1.*energy_factor, ehigh = 10.*energy_factor;
1040  VecIntegrator<my_Integrand_S62,VecGaussian32> S62(func);
1041  /* This returns a thermally averaged collision strength */
1042  double coll_str = S62.sum( elow, emid );
1043  coll_str += S62.sum( emid, ehigh );
1044  ASSERT( coll_str > 0. );
1045  coll_str /= energy_factor;
1046 
1047  // Moved constant factors in conversion from cross-section to collision strength out of integral,
1048  // hence 1.0 for energy argument here
1049  coll_str = ConvCrossSect2CollStr(coll_str*PI*BOHR_RADIUS_CM*BOHR_RADIUS_CM, (double)gLo, 1.0, reduced_mass);
1050  return coll_str;
1051 }
1052 
1053 //* Calculate beta(F), where F(beta) = K0(beta)^2 + K1(beta)^2
1054 STATIC double S62BesselInvert(double zOverB2)
1055 {
1056  DEBUG_ENTRY( "S62BesselInvert()" );
1057  double betaone;
1058  if( zOverB2 > 100. )
1059  {
1060  betaone = sqrt( 1./zOverB2 );
1061  }
1062  else if( zOverB2 < 0.54 )
1063  {
1064  /* Low betaone approximation */
1065  double logz = log(zOverB2/PI);
1066  betaone = (1./3.)*(1.28-logz);
1067  if( betaone > 2.38 )
1068  {
1069  /* average with this over approximation */
1070  betaone += -0.5*logz;
1071  betaone *= 0.5;
1072  }
1073  }
1074  else
1075  {
1076  long ip_zOverB2 = 0;
1077  static const double zetaOVERbeta2[101] = {
1078  1.030E+02,9.840E+01,9.402E+01,8.983E+01,8.583E+01,8.200E+01,7.835E+01,7.485E+01,
1079  7.151E+01,6.832E+01,6.527E+01,6.236E+01,5.957E+01,5.691E+01,5.436E+01,5.193E+01,
1080  4.961E+01,4.738E+01,4.526E+01,4.323E+01,4.129E+01,3.943E+01,3.766E+01,3.596E+01,
1081  3.434E+01,3.279E+01,3.131E+01,2.989E+01,2.854E+01,2.724E+01,2.601E+01,2.482E+01,
1082  2.369E+01,2.261E+01,2.158E+01,2.059E+01,1.964E+01,1.874E+01,1.787E+01,1.705E+01,
1083  1.626E+01,1.550E+01,1.478E+01,1.409E+01,1.343E+01,1.280E+01,1.219E+01,1.162E+01,
1084  1.107E+01,1.054E+01,1.004E+01,9.554E+00,9.094E+00,8.655E+00,8.234E+00,7.833E+00,
1085  7.449E+00,7.082E+00,6.732E+00,6.397E+00,6.078E+00,5.772E+00,5.481E+00,5.202E+00,
1086  4.937E+00,4.683E+00,4.441E+00,4.210E+00,3.989E+00,3.779E+00,3.578E+00,3.387E+00,
1087  3.204E+00,3.031E+00,2.865E+00,2.707E+00,2.557E+00,2.414E+00,2.277E+00,2.148E+00,
1088  2.024E+00,1.907E+00,1.795E+00,1.689E+00,1.589E+00,1.493E+00,1.402E+00,1.316E+00,
1089  1.235E+00,1.157E+00,1.084E+00,1.015E+00,9.491E-01,8.870E-01,8.283E-01,7.729E-01,
1090  7.206E-01,6.712E-01,6.247E-01,5.808E-01,5.396E-01};
1091 
1092  ASSERT( zOverB2 >= zetaOVERbeta2[100] );
1093 
1094  /* find beta in the table */
1095  long ilo=0, ihi = 100;
1096  for(;;)
1097  {
1098  long imid = (ilo+ihi)/2.;
1099  if ( zetaOVERbeta2[imid] > zOverB2 )
1100  ilo = imid;
1101  else
1102  ihi = imid;
1103  if (ihi == ilo+1)
1104  break;
1105  }
1106  ASSERT( ( zOverB2 < zetaOVERbeta2[ilo] ) && ( zOverB2 >= zetaOVERbeta2[ilo+1] ) );
1107 
1108  ip_zOverB2 = ilo;
1109 
1110  ASSERT( (ip_zOverB2 >=0) && (ip_zOverB2 < 100) );
1111 
1112  const double fp = 1.023292992280754131; // exp10(1./100.);
1113  betaone = exp10( (double)ip_zOverB2/100. - 1.)*
1114  ( (zOverB2 - zetaOVERbeta2[ip_zOverB2]) * fp
1115  -zOverB2 + zetaOVERbeta2[ip_zOverB2+1] )/
1116  (zetaOVERbeta2[ip_zOverB2+1] - zetaOVERbeta2[ip_zOverB2]);
1117  }
1118  //fprintf(ioQQQ,"%g %g %g\n",betaone,zOverB2,POW2(bessel_k0(betaone))+POW2(bessel_k1(betaone)));
1119  return betaone;
1120 }
1121 
1122 /* The integrand for calculating the thermal average of collision strengths */
1123 void my_Integrand_S62::operator() (const double proj_energy[], double res[], long n) const
1124 {
1125  DEBUG_ENTRY( "S62_Therm_ave_coll_str()" );
1126 
1127  ASSERT( deltaE > 0. );
1128 
1129  avx_ptr<double> x(n), val(n);
1130  for( long i=0; i < n; ++i )
1131  x[i] = -proj_energy[i]*EVDEGK/temp;
1132  vexp(x.ptr0(), val.ptr0(), 0, n);
1133 
1134  for( long i=0; i < n; ++i )
1135  {
1136  /* projectile energy in eV */
1137  /* Rnot = 1.1229*H_BAR/sqrt(ELECTRON_MASS*deltaE*EN1EV)/Bohr_radius; in units of Bohr_radius */
1138  /* The deltaE here is to make sure that the collider has no less
1139  * than the energy difference between the initial and final levels. */
1140  double Dubya = proj_energy[i] + 0.5*deltaE;
1141  ASSERT( Dubya > 0. );
1142 
1143  /* betanot = sqrt((proj_energy[i]+deltaE)/I_energy_eV)*(deltaE/2./Dubya)*Rnot; */
1144 
1145  ASSERT( I_energy_eV > 0. );
1146  ASSERT( osc_strength > 0. );
1147 
1148  // z/b1**2 = 0.5*W**2/(deltaE**2*Phi_ij) -- (33)
1149  // Phi_ij = I_energy_eV*osc_strength/deltaE -- (19)
1150  /* from equation 33 */
1151  double zOverB2 = 0.5*Dubya*Dubya/(deltaE*I_energy_eV*osc_strength);
1152 
1153  // zOverB2 = K_0^2+K_1^2 -- (11)
1154  ASSERT( zOverB2 > 0. );
1155 
1156  double betaone = S62BesselInvert(zOverB2);
1157 
1158  double zeta_of_betaone = zOverB2 * POW2(betaone);
1159 
1160  // Equation (39)
1161  /* cs1 = betanot * bessel_k0(betanot) * bessel_k1(betanot); */
1162  double k0val, k1val;
1163  bessel_k0_k1(betaone, &k0val, &k1val);
1164  double cs2 = 0.5*zeta_of_betaone + betaone * k0val * k1val;
1165 
1166  /* cross_section = MIN2(cs1, cs2); */
1167  double cross_section = cs2;
1168 
1169  /* cross section in units of PI * a_o^2 */
1170  cross_section *= 8. * (I_energy_eV/deltaE) * osc_strength * (I_energy_eV/(proj_energy[i]+deltaE));
1171 
1172  // Maxwell weighting cf Burgess & Tully 1992 A&A 254, 436 eq (21), Seaton (1953)
1173  // Moved constant factors in conversion from cross-section to collision strength out of integral
1174  res[i] = val[i] * (proj_energy[i]+deltaE)/EVRYD * cross_section;
1175  }
1176 }
1177 /*CS_l_mixing_PS64 - find rate for l-mixing collisions by protons, for neutrals */
1178 /* This version does not assume that E_min/kt is small and works with the exponential integral */
1179 /* Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165
1180  * assume Emin/KT << 1 which limites the range of application of the treatment.
1181  * This routine is based on N. Badnell developments on PS formulation.
1182  */
1183 double CS_l_mixing_PS64_expI(
1184  long nelem,
1185  double tau,
1186  double target_charge,
1187  long n,
1188  long l,
1189  double g,
1190  long lp,
1191  double deltaE_eV,
1192  long Collider)
1193 {
1194  double cs;
1195  double RD,R12,Sij,Plowb,RC,RC1,/*R1,*/EC,ED,
1196  reduced_mass, reduced_mass_2_emass,bracket,contr, rate;
1197  double fb1, fb2;
1198  double eEm,eED,eEC,eEmt1Em;
1199 
1200  DEBUG_ENTRY( "CS_l_mixing_PS64_expI()" );
1201 
1202  const double ChargIncoming = ColliderCharge[Collider];
1203  // Bohr time
1204  const double tau_zero = 2.41889e-17;
1205  long lmax = MAX2(l,lp);
1206  double n2 = (double)n*(double)n;
1207  double lmax2 = (double)lmax*(double)lmax;
1208  reduced_mass = reduced_amu(nelem, Collider);
1209  reduced_mass_2_emass = reduced_mass / ELECTRON_MASS;
1210 
1211  //Projectile energy
1212  double tempryd=phycon.te/TE1RYD;
1213  // density in au
1214  double dens_au = dense.eden*pow(BOHR_RADIUS_CM,3);
1215  // projectile velocity in au
1216  double vred = sqrt(tempryd/reduced_mass_2_emass);
1217  // lifetime in au
1218  double tau_ua = tau/tau_zero;
1219 
1220  Plowb = 0.5;
1221  fb1 = 1.;
1222  fb2 = 1.;
1223  /*line strength multiplied by 2 (this factor comes from the spin of the electron) */
1224  Sij = 9.*n2*lmax*(n2-lmax2)/(2.*target_charge*target_charge);
1225 
1226  /* set uf cut offs* and Emin*/
1227  /*First, R1, low impact parameter cut-off. R12 = R1^2*/
1228  R12 = 2. * pow2(ChargIncoming)/(3.*Plowb*vred*vred*(2.*l+1));
1229 
1230  R12 *= Sij;
1231  //R1 = sqrt(R12);
1232 
1233  /* Debye cut-off in au*/
1234  RD = sqrt(tempryd/(8.*PI*dens_au));
1235 
1236  /*Lifetime cut off in au*/
1237  RC1 = 0.72*tau_ua*vred;
1238 
1239  //degeneration energy dependent cuts off
1240  if (nelem > 0)
1241  {
1242  double meanb;
1243  double PS_crit = deltaE_eV*tau/HBAReV;
1244 
1245  //double aveRadius = (BOHR_RADIUS_CM/((double)nelem+1.-(double)ipHE_LIKE))*POW2((double)n);
1246  double bmax = sqrt( BOLTZMANN*phycon.te/dense.eden )
1247  / (PI2*ELEM_CHARGE_ESU);
1248  if (dense.xIonDense[0][0] > 0 )
1249  {
1250  // mean impact parameter is the minimum of 1/(ion density) or Debye radius
1251  meanb= powpq(1/dense.xIonDense[0][0],1,3);
1252  meanb = MIN2(meanb,bmax);
1253  }
1254  else
1255  meanb = bmax;
1256  double v = sqrt(2.*BOLTZMANN*phycon.te/reduced_mass);
1257  double beta_brock= meanb*abs(deltaE_eV)/(HBAReV*4.*PI*v);
1258 
1259  /*PS64 criterion for degeneration
1260  * if PS_crit = 1.55 Rc = 0.72*v*tau.
1261  * PS64 eq. 29 */
1262  if (PS_crit > 1.55 && iso_ctrl.lgCS_PSdeg[ipHE_LIKE])
1263  RC1 = 1.12*HBAReV/*EXPEULER2*/*vred/deltaE_eV/tau_zero;
1264 
1265  // Brocklehurst criterion on eq. 3.12 MNRAS (1972) 157, 211
1266  if (beta_brock >= 0.4 && iso_ctrl.lgCS_B72[ipHE_LIKE])
1267  RC1 = 1.12*HBAReV/*EXPEULER2*/*vred/deltaE_eV/tau_zero;
1268 
1269  /* Emin for energy dependent cut off */
1270  // Emin = R1/RC1;
1271  }
1272  RC = min(RC1,RD); /* use the minimum cut-off*/
1273 
1274  /* E/Emin=RC1**2/R12 or
1275  * E/Emin = Rc1/R1
1276  */
1277  //if (RC < sqrt(R12) )
1278  //{
1279  fb1 = 2./3.;
1280  R12 = R12*2.;
1281  //R1 = sqrt(R12);
1282  //RC = R1;
1283  //}
1284  double Emin = R12/(RC*RC);
1285 
1286  if (RC == RD)
1287  fb2=1.;
1288  else
1289  {
1290  fb2=2.;
1291  Emin = sqrt(Emin);
1292  }
1293 
1294  //Emin *= Eproj;
1295 
1296  /* Resolved Dnl depending on Sij */
1297  double Dnl = 2. * pow2(ChargIncoming)*Sij/(3.*(2.*l+1.));
1298 
1299  ASSERT( Dnl > 0. );
1300  ASSERT( phycon.te / Dnl / reduced_mass_2_emass > 0. );
1301 
1302  EC = RD*RD/(RC1*RC1);
1303  ED = R12/(RD*RD);
1304  eEm = exp(-1.*Emin);
1305  eED = exp(-1.*ED);
1306  eEC = exp(-1.*EC);
1307  eEmt1Em = eEm*(1.+Emin);
1308 
1309  /* First exponential integral is used as the analitical solution of Maxwell averaged PS64 cross sections
1310  * Different cases are used depending on the cut-off
1311  */
1312  if ( fb1 == 1 && fb2 == 1)
1313  bracket = eEm + e1(Emin);
1314  else
1315  {
1316  if (fb2 ==1 )
1317  bracket = fb1*eEm+ fb2*e1(Emin);
1318  else
1319  {
1320  if (EC > Emin)
1321  bracket = fb1*eEm + 2.*e1(Emin) - e1(EC);
1322  else
1323  bracket = fb1*eED + e1(ED);
1324  }
1325  }
1326 
1327  //contribution 0<E<Emin, important for Emin/kt >>
1328  contr = 0.;
1329  if(fb1 != 1 )
1330  {
1331  if ( fb2 == 1 )
1332  contr = (1.-eEmt1Em)/Emin;
1333  else if (fb2 !=1 )
1334  {
1335  if (EC >= Emin)
1336  {
1337  contr = 2.*( 1. -eEmt1Em)/pow2(Emin);
1338  contr -= eEm;
1339  }
1340  else
1341  {
1342  contr = ( 2. -eEC*(2.+EC))/pow2(Emin);
1343  contr -= eEm*(1.+1./ED);
1344  }
1345  }
1346 
1347  contr *= 2./3.;
1348 
1349  bracket += contr;
1350  }
1351 
1352 
1353  ASSERT( bracket >= 0.);
1354 
1355  if (bracket == 0. )
1356  return SMALLFLOAT;
1357 
1358  /* This is the rate coefficient. Units: cm^3 s-1 */
1359 
1360  double units = 2.*pow(BOHR_RADIUS_CM,3)*sqrt(PI)/vred/tau_zero;
1361 
1362  rate = units * Dnl* bracket;
1363 
1364  /* convert rate to collision strength */
1365  /* NB - the term in parentheses corrects for the fact that COLL_CONST is only appropriate
1366  * for electron colliders and is off by reduced_mass_2_emass^-1.5 */
1367  cs = rate / ( COLL_CONST * powpq(reduced_mass_2_emass, -3, 2) ) * phycon.sqrte * g;
1368 
1369  ASSERT( cs > 0. );
1370 
1371  return cs;
1372 }
1373 
1374 /* Classical PS64, refer to eq. 43 Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165
1375  CS_l_mixing_PS64 - find rate for l-mixing collisions by protons, for neutrals */
1376 double CS_l_mixing_PS64(
1377  long nelem, /* the chemical element, 1 for He */
1378  double tau, /* the radiative lifetime of the initial level. */
1379  double target_charge,
1380  long n,
1381  long l,
1382  double gLo,
1383  long lp,
1384  double deltaE_eV,
1385  long Collider)
1386 {
1387  /* >>refer H-like l-mixing Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */
1388  /* >>refer He-like l-mixing Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */
1389  double cs,
1390  TwoLogDebye, TwoLogRc1,
1391  factor1, factor2,
1392  bestfactor, factorpart,
1393  reduced_mass, reduced_mass_2_emass,
1394  rate;
1395  const double ChargIncoming = ColliderCharge[Collider];
1396 
1397  DEBUG_ENTRY( "CS_l_mixing_PS64()" );
1398 
1399  /* In this routine, two different cutoff radii are calculated, and as per PS64,
1400  * the least of these is selected. We take the least positive result. */
1401 
1402  /* Only dipole l-changing is considered in this approach */
1403  ASSERT( abs(l-lp) == 1);
1404  /* This reduced mass is in grams. */
1405  reduced_mass = reduced_amu(nelem, Collider);
1406  /* this mass always appears relative to the electron mass, so define it that way */
1407  reduced_mass_2_emass = reduced_mass / ELECTRON_MASS;
1408 
1409  /* equation 46 of PS64 */
1410  /* min on density added to prevent this from becoming large and negative
1411  * at very high n_e - Pengelly & Seaton did not apply this above 1e12 cm^-3 */
1412  /* This is 2 times the log of the Debye radius. */
1413  TwoLogDebye = 1.68 + log10( phycon.te / MIN2(1e11 , dense.eden ) );
1414  /* Brocklehurst (1971, equation 3.40) has 1.181 instead of 1.68. This appears to be due to
1415  * Pengelly and Seaton neglecting one of the two factors of PI in their Equation 33 */
1416  //TwoLogDebye = 1.181 + log10( phycon.te / MIN2(1e10 , dense.eden ) );
1417 
1418  /* This is 2 times the log of cutoff = 0.72v(tau), where tau is the lifetime of the level nl.
1419  * This is PS64 equation 45 (same as Brocklehurst 1971 equation 3.41) */
1420  TwoLogRc1 = 10.95 + log10( phycon.te * tau * tau / reduced_mass_2_emass );
1421 
1422  /* non-degenerate case */
1423  if (nelem == 1)
1424  {
1425  double meanb;
1426  double PS_crit = deltaE_eV*tau/HBAReV;
1427 
1428  // This is Debye Radius
1429  double bmax = sqrt( BOLTZMANN*phycon.te/dense.eden/PI )
1430  / (2*ELEM_CHARGE_ESU);
1431  double vred = sqrt(2.*BOLTZMANN*phycon.te/reduced_mass);
1432 
1433  if (dense.xIonDense[0][0] > 0 )
1434  {
1435  // mean impact parameter is the minimum of 1/(ion density) or Debye radius
1436  meanb= powpq(1/dense.xIonDense[0][0],1,3);
1437  meanb = MIN2(meanb,bmax);
1438  }
1439  else
1440  meanb = bmax;
1441 
1442  double beta_brock= meanb*abs(deltaE_eV)/(HBAReV*4.*PI*vred);
1443  double deltaE_cm = deltaE_eV*EN1EV/ERG1CM;
1444 
1445  /*PS64 criterion for degeneration
1446  * if PS_crit = 1.55 Rc = 0.72*v*tau.
1447  * PS64 eq. 29
1448  */
1449  if (PS_crit > 1.55 && iso_ctrl.lgCS_PSdeg[ipHE_LIKE])
1450  /* Direct calculation */
1451  //TwoLogRc1 = 2.*log10(1.12*HBAReV*vred/**EXPEULER2*//deltaE_eV);
1452  /* Eq. 3.13 Brocklehurst 1971 for degenerate case */
1453  TwoLogRc1 = log10(phycon.te*ELECTRON_MASS/pow2(deltaE_cm)/reduced_mass)-11.22;
1454 
1455  // Brocklehurst criterion on eq. 3.12 MNRAS (1972) 157, 211
1456  if (beta_brock >= 0.4 && iso_ctrl.lgCS_B72[ipHE_LIKE])
1457  /* Direct calculation */
1458  //TwoLogRc1 = 2.*log10(1.12*HBAReV*vred*EXPEULER2/deltaE_eV);
1459  /* Eq. 3.13 Brocklehurst 1971 for degenerate case */
1460  TwoLogRc1 = log10(phycon.te*ELECTRON_MASS/pow2(deltaE_cm)/reduced_mass)-11.22;
1461 
1462  }
1463 
1464  /* this is equation 44 of PS64
1465  * unresolved Dnl
1466  */
1467  double Dnl = POW2( ChargIncoming / target_charge) * 6. * POW2( (double)n) *
1468  ( POW2((double)n) - POW2((double)l) - l - 1);
1469 
1470 
1471  /***** NB NB NB NB
1472  * Brocklehurst (1971) has a factor of electron density in the rate coefficient (equation 3.38).
1473  * This is NOT a proper rate, as can be seen in his equations 2.2 and 2.4. This differs
1474  * from the formulism given by PS64. */
1475  //rate *= dense.eden;
1476 
1477 
1478  // Separate up and down rates in PS64 expression using
1479  // (2l+1) Q_{l,l-1} = (2l-1) Q_{l-1,l} = n^2 l - l^3
1480  // So (2l+1) Dnl' = (n^2 l - l^3) + (n^2 [l+1] - [l+1]^3) = (2l+1) (n^2 - l^2 - l - 1)
1481 
1482  // l, gLo are value for *initial* state. It is likely (but not
1483  // certain) that l < lp.
1484 
1485  long lmax = MAX2(l,lp);
1486  /* This is resolved l->lp Dnl */
1487  // double Dnlup = POW2( ChargIncoming / target_charge) * 6. * POW2( (double)n) *
1488  // lmax * (n*n-lmax*lmax) / double( 2*l+1 );
1489  double Dnlup = POW2( ChargIncoming / target_charge) * 6. * POW2( (double)n) * lmax * (n*n-lmax*lmax) / double( 2*l+1 );
1490 
1491 
1492  ASSERT( Dnl > 0. );
1493  ASSERT( phycon.te / Dnl / reduced_mass_2_emass > 0. );
1494 
1495  /* In the current configuration is the resolved/unresolved Dnl what get's inside the logarithm.
1496  * That depends on the implementation. Resolved Dnl has been used in F. Guzman MNRAS (2016) 459, 3498
1497  * Summers MNRAS (1977) uses unresolved Dnl. That allows to use resolved results as a branching ratio
1498  * of total unresolved PS64 crates.
1499  */
1500  factorpart = (11.54 + log10( phycon.te / Dnl / reduced_mass_2_emass ) );
1501 
1502  if( (factor1 = factorpart + TwoLogDebye) <= 0.)
1503  factor1 = BIGDOUBLE;
1504 
1505  if( (factor2 = factorpart + TwoLogRc1) <= 0.)
1506  factor2 = BIGDOUBLE;
1507 
1508  /* Now we find the least positive result. */
1509  bestfactor = MIN2(factor1,factor2);
1510 
1511  ASSERT( bestfactor > 0. );
1512 
1513  /* If both factors are bigger than 100, toss out the result, and return SMALLFLOAT. */
1514  if( bestfactor > 100. )
1515  {
1516  return SMALLFLOAT;
1517  }
1518 
1519  /* This is the rate coefficient. Units: cm^3 s-1 */
1520  rate = 9.93e-6 * sqrt( reduced_mass_2_emass ) * Dnlup / phycon.sqrte * bestfactor;
1521 
1522  /* convert rate to collision strength */
1523  /* NB - the term in parentheses corrects for the fact that COLL_CONST is only appropriate
1524  * for electron colliders and is off by reduced_mass_2_emass^-1.5 */
1525  cs = rate / ( COLL_CONST * powpq(reduced_mass_2_emass, -3, 2) ) *
1526  phycon.sqrte * gLo;
1527 
1528  ASSERT( cs > 0. );
1529 
1530  return cs;
1531 }
1532 
1533 /*CS_l_mixing - find collision strength for l-mixing collisions for neutrals */
1534 /* The VF stands for Vrinceanu & Flannery 2001 */
1535 /* >>refer He l-mixing Vrinceanu, D. & Flannery, M. R. 2001, PhysRevA 63, 032701 */
1536 /* >>refer He l-mixing Hezel, T. P., Burkhardt, C. E., Ciocca, M., He, L-W., */
1537 /* >>refercon Leventhal, J. J. 1992, Am. J. Phys. 60, 329 */
1538 
1539 template<class P>
1540 inline double CS_l_mixing(long ipISO,
1541  long nelem,
1542  long n,
1543  long l,
1544  long lp,
1545  long s,
1546  long gLo,
1547  double tauLo,
1548  double IP_Ryd_Hi,
1549  double IP_Ryd_Lo,
1550  double temp,
1551  long Collider )
1552 {
1553  double coll_str;
1554 
1555  DEBUG_ENTRY( "CS_l_mixing()" );
1556 
1557  typedef my_Integrand_VF01_E<P> integType;
1558  integType func(ipISO,nelem,n,l,lp,s,gLo,tauLo,IP_Ryd_Hi,IP_Ryd_Lo,Collider,temp);
1559 
1560  /* no need to do this for h-like */
1561  if( iso_ctrl.lgCS_Seaton[ipISO] && ipISO > ipH_LIKE )
1562  {
1563  ASSERT( l != 0 );
1564  ASSERT( lp != 0 );
1565  }
1566 
1567  if( !iso_ctrl.lgCS_therm_ave[ipISO] )
1568  {
1569  /* Do NOT average over Maxwellian */
1570  /* Compensate for the exponential weighting factor applied in
1571  * the integrand */
1572  coll_str = func(1.0)*exp(1.0);
1573  }
1574  else
1575  {
1576  /* DO average over Maxwellian */
1577  if (1)
1578  {
1579  Integrator<integType, Gaussian32> VF01_E(func);
1580 
1581  coll_str = VF01_E.sum( 0.0, 1.0 );
1582  coll_str += VF01_E.sum( 1.0, 10.0 );
1583  }
1584  else
1585  {
1586  double xn=0., xl = 10.;
1587  if (func(xl) != 0.)
1588  {
1589  while (xl < 100. && func(xl/2.) == 0.)
1590  {
1591  xn = xl/2.;
1592  xl *= 2.;
1593  }
1594  }
1595 
1596  class integrate::Romberg<integrate::Midpoint<integType> >
1597  VF01_ER(integrate::Midpoint<integType>(func,xn,xl));
1598  VF01_ER.update(1e-3);
1599  coll_str = VF01_ER.sum();
1600 
1601  if (0)
1602  {
1603  if (0)
1604  {
1605  typedef my_Integrand_VF01_E_log<P> integLogType;
1606  integLogType funcl(ipISO,nelem,n,l,lp,s,gLo,tauLo,IP_Ryd_Hi,IP_Ryd_Lo,Collider,temp);
1607  double x1n = 0.0, x1x = 1.0;
1608  do {
1609  double x1m = 0.5*(x1n+x1x);
1610  if (funcl(x1m) > 0.)
1611  x1n = x1m;
1612  else
1613  x1x = x1m;
1614  } while ((x1x-x1n) > 0.001);
1615  class integrate::Simple<integrate::Midpoint<integLogType> >
1616  VF01_ESL(integrate::Midpoint<integLogType>(funcl,0.0,x1x));
1617  VF01_ESL.update(1e-3);
1618 
1619  Integrator<integType, Gaussian32> VF01_E(func);
1620  double coll_strg = VF01_E.sum( 0.0, 1.0 );
1621  coll_strg += VF01_E.sum( 1.0, 10.0 );
1622 
1623  class integrate::Romberg<integrate::Midpoint<integLogType> >
1624  VF01_ERL(integrate::Midpoint<integLogType>(funcl,0.0,x1x));
1625  VF01_ERL.update(1e-3);
1626  fprintf(ioQQQ,"Int %ld %ld->%ld ER %g(%g) %ld ESL %g(%g) %ld ERL %g(%g) %ld G %g S %g\n",
1627  n,l,lp,
1628  VF01_ER.sum(),VF01_ER.error(),VF01_ER.evals(),
1629  VF01_ESL.sum(),VF01_ESL.error(),VF01_ESL.evals(),
1630  VF01_ERL.sum(),VF01_ERL.error(),VF01_ERL.evals(),
1631  coll_strg,func(1.0)
1632  );
1633  fflush(ioQQQ);
1634  }
1635  else
1636  {
1637  fprintf(ioQQQ,"Int %ld %ld->%ld ER %g(%g) %ld\n",
1638  n,l,lp,
1639  coll_str,VF01_ER.error(),VF01_ER.evals()
1640  );
1641  fflush(ioQQQ);
1642  }
1643  }
1644  }
1645  }
1646 
1647  return coll_str;
1648 }
1649 
1650 double CS_l_mixing_VF01(long ipISO,
1651  long nelem,
1652  long n,
1653  long l,
1654  long lp,
1655  long s,
1656  long gLo,
1657  double tauLo,
1658  double IP_Ryd_Hi,
1659  double IP_Ryd_Lo,
1660  double temp,
1661  long Collider )
1662 {
1663 
1664  return CS_l_mixing<StarkCollTransProb_VF01>(ipISO,nelem,n,l,lp,s,gLo,tauLo,
1665  IP_Ryd_Hi,IP_Ryd_Lo,temp,Collider);
1666 
1667 }
1668 
1669 double CS_l_mixing_VOS12QM(long ipISO,
1670  long nelem,
1671  long n,
1672  long l,
1673  long lp,
1674  long s,
1675  long gLo,
1676  double tauLo,
1677  double IP_Ryd_Hi,
1678  double IP_Ryd_Lo,
1679  double temp,
1680  long Collider )
1681 {
1682  return CS_l_mixing<StarkCollTransProb_VOS12QM>(ipISO,nelem,n,l,lp,s,gLo,
1683  tauLo,IP_Ryd_Hi,IP_Ryd_Lo,temp,Collider);
1684 }
1685 
1686 namespace
1687 {
1688  class Chi_VF01
1689  {
1690  double m_alpha, m_alphamin, m_sinfac;
1691  public:
1692  Chi_VF01(double alpha, double alphamin) : m_alpha(alpha), m_alphamin(alphamin)
1693  {
1694  ASSERT( m_alpha >= 0. );
1695  /* deltaPhi is the effective angle swept out by the projector as viewed by the target. */
1696  /* deltaPhi is PI for a complete collision in a straight trajectory
1697  * Kazansky et al. JPhysB: At. Mol. Opt. Phys. 29, 3651 proposed a core effect correction
1698  * This has an unwanted effect in the probability as it decreases at high impact parameters
1699  * That shouldn't happen as deflection should be stronger at low impact parameters*/
1700 
1701  //double b_over_bmax = m_alphamin / m_alpha;
1702  double deltaPhi = -PI;//2.*asin(b_over_bmax)-PI;
1703  m_sinfac = pow2(sin(0.5*deltaPhi*sqrt(1+m_alpha*m_alpha)));
1704  }
1705  double sinChiOver2sq() const
1706  {
1707  /* >>refer He l-mixing Kazansky, A. K. & Ostrovsky, V. N. 1996, JPhysB: At. Mol. Opt. Phys. 29, 3651 */
1708  if( m_alpha <= m_alphamin)
1709  {
1710  return 0.;
1711  }
1712  double denom = (1.+m_alpha*m_alpha), ratio = m_alpha/denom;
1713  return 4.*ratio*ratio*m_sinfac*(denom-m_alpha*m_alpha*m_sinfac);
1714  }
1715  double cosChi() const
1716  {
1717  double alphasq = pow2(m_alpha);
1718  return 1.-2.*alphasq*m_sinfac/(1.+alphasq);
1719  }
1720  };
1721 
1722  class StarkCollTransProb_VF01
1723  {
1724  long int n, l, lp;
1725  double cosU1, cosU2, sinU1sinU2;
1726  public:
1727  StarkCollTransProb_VF01( long int n1, long int l1, long int lp1) : n(n1), l(l1), lp(lp1)
1728  {
1729  /* These are defined on page 11 of VF01 */
1730  cosU1 = 2.*POW2((double)l/(double)n) - 1.;
1731  cosU2 = 2.*POW2((double)lp/(double)n) - 1.;
1732  sinU1sinU2 = sqrt( ( 1. - cosU1*cosU1 )*( 1. - cosU2*cosU2 ) );
1733  }
1734  double operator()(double alpha, double alphamin) const;
1735  double bfun(double alpha, double alphamin) const
1736  {
1737  double sinChiOver2sq = Chi_VF01(alpha, alphamin).sinChiOver2sq();
1738  double cosChi = 1. - 2*sinChiOver2sq;
1739  return (cosU1*cosU2 + sinU1sinU2 - cosChi);
1740  }
1741  };
1742 
1743  double StarkCollTransProb_VF01::operator() (double alpha, double alphamin) const
1744  {
1745  DEBUG_ENTRY( "StarkCollTransProb_VF01::operator()()" );
1746  double probability;
1747  ASSERT( alpha >= 1e-30 );
1748 
1749  double sinChiOver2sq = Chi_VF01(alpha, alphamin).sinChiOver2sq();
1750 
1751  /*if ( lp == 0 )
1752  // Refer to Vrinceanu and Flanery (2001)PRA 63, 032701
1753  // before eq. (35) "for l->lp=0 transitions, the probability is zero."
1754  probability = 0.;
1755  else */
1756  if( l == 0 || lp == 0 )
1757  {
1758  long lf = max(l,lp);
1759  if( n*n*sinChiOver2sq <= lf*lf ) //lf*lf ) //lp*lp )
1760  {
1761  probability = 0.;
1762  }
1763  else
1764  {
1765  /* Here is the initial state S case */
1766  /*ASSERT( sinChiOver2sq > 0. );
1767  ASSERT( n*n*sinChiOver2sq > lp*lp ); //lf*lf ); //lp*lp );*/
1768  /* This is equation 35 of VF01. There is a factor of hbar missing in the denominator of the
1769  * paper, but it's okay if you use atomic units (hbar = 1). */
1770  //probability = lp+0.5 / sqrt( n*n*sinChiOver2sq * (n*n*sinChiOver2sq - lp*lp) );
1771  probability = (lf+0.5) / sqrt( n*n*sinChiOver2sq * (n*n*sinChiOver2sq - lf*lf) );
1772  }
1773 
1774  if (lp == 0)
1775  {
1776  probability /= (2.*l+1);
1777  }
1778  }
1779  else
1780  {
1781 #if 0
1782  // Rearrangement of cosChi above means that this is no longer necessary.
1783  if( OneMinusCosChi == 0. )
1784  {
1785  double hold = sin( deltaPhi / 2. );
1786  /* From approximation at bottom of page 10 of VF01. */
1787  OneMinusCosChi = 8. * alpha * alpha * POW2( hold );
1788  }
1789 #endif
1790 
1791  if( sinChiOver2sq == 0. )
1792  {
1793  probability = 0.;
1794  }
1795  else
1796  {
1797  double cosChi = 1. - 2*sinChiOver2sq;
1798  /* Here is the general case, with factor of OneMinusCosChi cancelled for efficiency */
1799  double A = (cosU1*cosU2 - sinU1sinU2 - cosChi);
1800  double B = (cosU1*cosU2 + sinU1sinU2 - cosChi);
1801 
1802  /* These are the three cases of equation 34. */
1803  if( B <= 0. )
1804  {
1805  probability = 0.;
1806  }
1807  else
1808  {
1809  ASSERT( sinChiOver2sq > 0. );
1810 
1811  probability = 2.*lp/(PI* /*H_BAR* */ n*n* sinChiOver2sq );
1812 
1813  /*lp factor in eq (34) of VF 2001 has been changed to lp+1/2
1814  * in order to cope with statistical factors
1815  */
1816  //probability = (2.*lp+1)/(PI* /*H_BAR* */ n*n* sinChiOver2sq );
1817 
1818  if( A < 0. )
1819  {
1820  ASSERT( B > A );
1821  probability *= ellpk( -A/(B-A) ) * sqrt( 2*sinChiOver2sq / (B - A) );
1822  }
1823  else
1824  {
1825  probability *= ellpk( A/B ) * sqrt( 2*sinChiOver2sq / B );
1826  }
1827  }
1828  }
1829  }
1830 
1831  return probability;
1832  }
1833 
1834  template<class P>
1835  class my_Integrand_VF01_alpha : public D_fp
1836  {
1837  P s;
1838  double alphamin;
1839  public:
1840  my_Integrand_VF01_alpha(long int n, long int l, long int lp, double alphamin) :
1841  s( n, l, lp), alphamin(alphamin) {}
1842 
1843  double operator() (double alpha) const
1844  {
1845  ASSERT( alpha >= 1.e-30 );
1846 
1847  return s( alpha, alphamin )/(alpha*alpha*alpha);
1848  }
1849  double bfun (double alpha) const
1850  {
1851  return s.bfun(alpha,alphamin);
1852  }
1853  };
1854 
1855  template<class P>
1856  class my_Integrand_VF01_beta : public D_fp
1857  {
1858  P s;
1859  double alphamin;
1860  public:
1861  // Change variables to beta = log(alpha) to improve accuracy & convergence
1862  // rate of Romberg integration
1863  my_Integrand_VF01_beta(long int n, long int l, long int lp, double alphamin) :
1864  s( n, l, lp), alphamin(alphamin) {}
1865 
1866  double operator() (double beta) const
1867  {
1868  double alpha = exp(beta);
1869  ASSERT( alpha >= 1.e-30 );
1870 
1871  return s( alpha, alphamin )/(alpha*alpha);
1872  }
1873  double bfun (double alpha) const
1874  {
1875  return s.bfun(alpha,alphamin);
1876  }
1877  };
1878 
1879  class StarkCollTransProb_VOS12QM
1880  {
1881  long int n, l, lp;
1882  double cosU1, cosU2, sinU1sinU2;
1883  vector<double> common;
1884  public:
1885  StarkCollTransProb_VOS12QM( long int n1, long int l1, long int lp1) : n(n1), l(l1), lp(lp1), common(n)
1886  {
1887  /* These are defined on page 11 of VF01 */
1888  cosU1 = 2.*POW2((double)l/(double)n) - 1.;
1889  cosU2 = 2.*POW2((double)lp/(double)n) - 1.;
1890 
1891  sinU1sinU2 = sqrt( ( 1. - cosU1*cosU1 )*( 1. - cosU2*cosU2 ) );
1892  long Lmax = MIN2(n-1,lp+l), Lmin = abs(lp-l);
1893  double ffac = 1./sqrt(double(n));
1894  long nsq = n*n;
1895  for (long L=1; L < Lmin; ++L)
1896  {
1897  ffac *= L/sqrt(double(nsq-L*L));
1898  }
1899  if (ffac <= 0.)
1900  {
1901  fprintf(ioQQQ," PROBLEM: Catastrophic underflow in VOS12 QM transition probability\n"
1902  " Try reducing resolved levels\n");
1903  cdEXIT(EXIT_FAILURE);
1904  }
1905 
1906  const int sjtyp = 1;
1907  if (sjtyp == 3)
1908  {
1909  double l1min,l1max,lmatch;
1910  long ier;
1911  rec6j(get_ptr(common)+Lmin,lp,l,0.5*(n-1),0.5*(n-1),0.5*(n-1),
1912  &l1min,&l1max,&lmatch,n-Lmin,&ier);
1913  ASSERT(Lmin == int(l1min) && ier >= 0);
1914  }
1915  for (long L=Lmin; L <= Lmax; ++L)
1916  {
1917  if (L > 0)
1918  ffac *= L/sqrt(double(nsq-L*L));
1919  if (0)
1920  {
1921  double factL = factorial(L);
1922  fprintf(ioQQQ,"TEST %ld %g %g\n",
1923  L,ffac,factL*factL*factorial(n-L-1)/factorial(n+L));
1924  }
1925  double sixj1;
1926  // pass 2*j arguments to allow for half-integer values
1927  if (sjtyp == 1)
1928  sixj1 = sjs(2*lp,2*l,2*L,n-1,n-1,n-1);
1929  else if (sjtyp == 2)
1930  sixj1 = SixJFull(2*lp,2*l,2*L,n-1,n-1,n-1);
1931  else
1932  sixj1 = common[L];
1933  common[L] = sixj1*ffac*sqrt(2.*double(L)+1.);
1934  }
1935  }
1936  double operator()(double alpha, double alphamin) const;
1937  double bfun (double, double ) const
1938  {
1939  return 0.;
1940  }
1941  };
1942 
1943  double StarkCollTransProb_VOS12QM::operator() (double alpha, double alphamin) const
1944  {
1945  DEBUG_ENTRY( "StarkCollTransProb_VOS12QM::operator()()" );
1946  ASSERT( alpha >= 1e-30 );
1947 
1948  Chi_VF01 chi(alpha, alphamin);
1949  double sinChiOver2sq = chi.sinChiOver2sq();
1950 
1951  // Implement equation (2) of VOS12 -- QM probability
1952  // Definition of chi appears inconsistent between VF01 &
1953  // VOS12. Compare the last display expression in the LH
1954  // column of VF01 032701-10 with VOS12 eq (3)
1955  //ASSERT (sinChiOver2sq <= 1.);
1956 
1957  double cosChiVOS12QM = chi.cosChi();
1958  ASSERT( cosChiVOS12QM <= 1);
1959  //double sinChiOver2sq = 1 - pow2(cosChiVOS12QM);
1960  long Lmax = MIN2(n-1,lp+l);
1961  UltraGen u(n,cosChiVOS12QM);
1962  for (long L=n-1; L > Lmax; --L)
1963  {
1964  u.step();
1965  }
1966  double frac = 4*sinChiOver2sq;
1967  double pqm = 0.;
1968  for (long L=Lmax; L >= abs(lp-l); --L)
1969  {
1970  double gegen = u.val();
1971  u.step();
1972  if (0)
1973  {
1974  // Test relative accuracy of different
1975  // Gegenbauer/ultraspherical implementations
1976  double gegen1 = gegenbauer(n-L-1,L+1,cosChiVOS12QM);
1977  if (! fp_equal_tol(gegen1,gegen,
1978  1000*DBL_EPSILON+1e-10*fabs(gegen1)))
1979  fprintf(ioQQQ,"DEBUG %ld %ld %g: %g %g %g\n",n,L,cosChiVOS12QM,
1980  gegen1,gegen,gegen/gegen1-1.);
1981  }
1982  pqm *= frac;
1983  double fac = common[L]*gegen;
1984  pqm += fac*fac;
1985  }
1986  pqm *= powi(frac,abs(lp-l))*(2*lp+1);
1987  ASSERT (pqm >= 0);
1988 
1989  return pqm;
1990  }
1991 
1992  void compare_VOS12()
1993  {
1994  //This function allows to compare probabilities of VOS12 approaches
1995  DEBUG_ENTRY("compare_VOS12()");
1996  //VOS12, Fig 1.
1997  {
1998  long n=30, l = 4, lp = 3, npt = 1000;
1999  double alphamin=0.;
2000  StarkCollTransProb_VOS12QM qm( n, l, lp);
2001  StarkCollTransProb_VF01 cl( n, l, lp);
2002  double v = 0.25;
2003  for (long i=0; i<npt; ++i)
2004  {
2005  double ban = 900*(i+0.5)/double(npt);
2006  double alpha = 1.5/(ban*v);
2007  fprintf(ioQQQ,"%g %g %g\n",ban,
2008  cl(alpha,alphamin),qm(alpha,alphamin));
2009  }
2010  }
2011  fprintf(ioQQQ,"\n\n");
2012  //VF01, Fig 16. -- one panel
2013  {
2014  long n=28, l = 5, lp = 13, npt = 1000;
2015  double alphamin=0.;
2016  StarkCollTransProb_VOS12QM qm( n, l, lp);
2017  StarkCollTransProb_VF01 cl( n, l, lp);
2018  for (long i=0; i<npt; ++i)
2019  {
2020  double alpha = (i+0.5)/double(npt);
2021  fprintf(ioQQQ,"%g %g %g\n",alpha,
2022  cl(alpha,alphamin),qm(alpha,alphamin));
2023  }
2024  }
2025  fprintf(ioQQQ,"\n\n");
2026 
2027  {
2028  long n=8, l = 5, lp = 6, npt = 1000;
2029  double alphamin=0.;
2030  vector<StarkCollTransProb_VOS12QM> qm;
2031  for (long i=0; i<4; ++i)
2032  {
2033  long scale = 1<<i;
2034  qm.push_back(StarkCollTransProb_VOS12QM(
2035  scale*n, scale*l, scale*lp));
2036  }
2037  StarkCollTransProb_VF01 cl( n, l, lp);
2038  for (long i=0; i<npt; ++i)
2039  {
2040  double alpha = (i+0.5)/double(npt);
2041  fprintf(ioQQQ,"%g %g",alpha,cl(alpha,alphamin));
2042  long scale = 1;
2043  for (vector<StarkCollTransProb_VOS12QM>::iterator it=qm.begin();
2044  it != qm.end(); ++it)
2045  {
2046  fprintf(ioQQQ," %g",scale*(*it)(alpha,alphamin));
2047  scale <<= 1;
2048  }
2049  fprintf(ioQQQ,"\n");
2050  }
2051  }
2052  fprintf(ioQQQ,"\n\n");
2053  //Compare all rates out of level
2054  {
2055  long n=30, l = 29, npt = 10000;
2056  double alphamin=0.;
2057  vector<StarkCollTransProb_VOS12QM> qm;
2058  for (long lp=0; lp<n; ++lp)
2059  qm.push_back(StarkCollTransProb_VOS12QM( n, l, lp));
2060  double v = 0.25;
2061  for (long i=0; i<npt; ++i)
2062  {
2063  double ban = 400*(i+0.5)/double(npt);
2064  double alpha = 1.5/(ban*v);
2065  fprintf(ioQQQ,"PC %g",ban);
2066  for (long lp=0; lp<n; ++lp)
2067  fprintf(ioQQQ," %g",
2068  qm[lp](alpha,alphamin));
2069  fprintf(ioQQQ,"\n");
2070  }
2071  }
2072  }
2073 }
2074 
2078 double CS_l_mixing_VOS12(long n, long l, long lp,
2079  long nelem, double gLo, long Ztarget, long Collider, double sqrte)
2080 {
2081  // Equation (9) of Vrinceau etal ApJ 747, 56 (2012)
2082  long lmin = (lp < l) ? lp : l, dl = abs(l-lp);
2083  double nlfac = double(n*n*(n*n*(l+lp)-lmin*lmin*(l+lp+2*dl)))/
2084  ((l+0.5)*dl*dl*dl);
2085  long Z = ColliderCharge[Collider];
2086  double massratio = reduced_amu(nelem,Collider)/ELECTRON_MASS;
2087  double rate = 1.294e-5*sqrt(massratio)*Z*Z/(Ztarget*Ztarget*sqrte) * nlfac;
2088  double cs = rate / ( COLL_CONST * powpq(massratio, -3, 2) ) * sqrte * gLo;
2089  return cs;
2090 }
2091 
2092 template<class P>
2093 STATIC double CSIntegral_QG32(const my_Integrand_VF01_alpha<P>& func,
2094  double alphamin_int, double alphamax_int)
2095 {
2096  if (alphamin_int >= alphamax_int)
2097  return 0.;
2098  double CSIntegral = 0.;
2099  Integrator<my_Integrand_VF01_alpha<P>, Gaussian32> VF01_alpha(func);
2100  double step = (alphamax_int - alphamin_int)/5.;
2101  double alpha1 = alphamin_int;
2102  CSIntegral += VF01_alpha.sum( alpha1, alpha1+step );
2103  CSIntegral += VF01_alpha.sum( alpha1+step, alpha1+4.*step );
2104  return CSIntegral;
2105 }
2106 
2107 template<class P>
2108 STATIC double CSIntegral_Romberg(long ipISO, const my_Integrand_VF01_beta<P>& func,
2109  double alphamin_int, double alphamax_int, double eps)
2110 {
2111  if (alphamin_int >= alphamax_int)
2112  return 0.;
2113 
2114  // Bisection search for lower limit to integral, should be
2115  // improved to e.g. secant method
2116  double alphalo = alphamin_int;
2117  if (! iso_ctrl.lgCS_VOS12QM[ipISO] && func.bfun(alphalo) <= 0.)
2118  {
2119  double alphahi = alphamax_int;
2120  while (alphahi-alphalo > 0.5*eps*(alphahi+alphalo))
2121  {
2122  double alphamid = 0.5*(alphahi+alphalo);
2123  if (func.bfun(alphamid) <= 0.)
2124  alphalo = alphamid;
2125  else
2126  alphahi = alphamid;
2127  }
2128  // Move edge *just* inside integration range
2129  alphalo = alphahi;
2130  }
2131  class integrate::Romberg<integrate::Midpoint<my_Integrand_VF01_beta<P> > >
2132  VF01_alphaR(integrate::Midpoint<my_Integrand_VF01_beta<P> >(func,log(alphalo),log(alphamax_int)));
2133  VF01_alphaR.update(eps);
2134  return VF01_alphaR.sum();
2135 }
2136 
2137 template<class P>
2138 STATIC double CSIntegral_Romberg_alpha(long ipISO, const my_Integrand_VF01_alpha<P>& func,
2139  double alphamin_int, double alphamax_int, double eps)
2140 {
2141  if (alphamin_int >= alphamax_int)
2142  return 0.;
2143  // Bisection search for lower limit to integral, should be
2144  // improved to e.g. secant method
2145  double alphalo = alphamin_int;
2146  if (! iso_ctrl.lgCS_VOS12QM[ipISO] && func.bfun(alphalo) <= 0.)
2147  {
2148  double alphahi = alphamax_int;
2149  while (alphahi-alphalo > 0.5*eps*(alphahi+alphalo))
2150  {
2151  double alphamid = 0.5*(alphahi+alphalo);
2152  if (func.bfun(alphamid) <= 0.)
2153  alphalo = alphamid;
2154  else
2155  alphahi = alphamid;
2156  }
2157  // Move edge *just* inside integration range
2158  alphalo = alphahi;
2159  }
2160  class integrate::Romberg<integrate::Midpoint<my_Integrand_VF01_alpha<P> > >
2161  VF01_alphaR(integrate::Midpoint<my_Integrand_VF01_alpha<P> >(func,alphalo,alphamax_int));
2162  VF01_alphaR.update(eps);
2163  return VF01_alphaR.sum();
2164 }
2165 
2166 template<class P>
2167 STATIC double collision_strength_VF01( long ipISO, double E_Proj_Ryd,
2168  const my_Integrand_VF01_E<P>& vf)
2169 {
2170 
2171  DEBUG_ENTRY( "collision_strength_VF01()" );
2172 
2173  double reduced_b_max1;
2174  /* The proton mass is needed here because the ColliderMass array is
2175  * expressed in units of the proton mass, but here we need absolute mass. */
2176  double reduced_vel = sqrt( 2.*E_Proj_Ryd*EN1RYD/vf.reduced_mass )/vf.RMSv;
2177 
2178  /* put limits on the reduced velocity. These limits should be more than fair. */
2179  ASSERT( reduced_vel > 1.e-10 );
2180  ASSERT( reduced_vel < 1.e10 );
2181 
2182  /* Reduced here means in units of vf.aveRadius: */
2183  //double reduced_b_min = 1.5 * vf.ColliderCharge / reduced_vel;
2184  /* alphamax has been increased in order to account for low cut-offs
2185  * needed at low temperatures T<100 and high densities
2186  */
2187  //if (iso_ctrl.lgCS_VOS12QM[ipISO])
2188  double alphamax = 15000./pow2(vf.n);
2189  /* alphamax can be fixed in different ways. Next is done assuming that cross section is zero at scaled
2190  * temperatures T=1K are zero so the lower limit= higher limit*/
2191  /*double alphamax = 1.5*vf.ColliderCharge/(reduced_vel * vf.reduced_b_max)*
2192  sqrt(E_Proj_Ryd*phycon.te/BOLTZMANN);*/
2193  // before was 1.5*ColliderCharge/(reduced_vel * reduced_b_min); but above
2194  // definition of reduced_b_min meant this is always 1. to machine
2195  // accuracy
2196  // this alphamax was not enough to cover low temperature cross sections
2197  if (iso_ctrl.lgCS_Vrinceanu[ipISO])
2198  alphamax = 1.;
2199  //
2200 
2201 
2202  /*reduced_b_min is used only to be compared to reduced_b_max*/
2203  double reduced_b_min = 1.5 * vf.ColliderCharge /reduced_vel/alphamax;
2204  ASSERT( reduced_b_min > 1.e-10 );
2205  ASSERT( reduced_b_min < 1.e10 );
2206 
2207  /* Here cut-off, set for temperature, is scaled to real energy for Maxwell calculation */
2208  // Debye Radius, reduced in units of aveRadius
2209  double reduced_debye = sqrt( BOLTZMANN*vf.temp/vf.ColliderCharge/dense.eden/PI )
2210  / (2.*ELEM_CHARGE_ESU)/vf.aveRadius;
2211  if (vf.reduced_b_max < reduced_debye && iso_ctrl.lgCS_VOS12QM[ipISO])
2212  reduced_b_max1 = vf.reduced_b_max*sqrt(E_Proj_Ryd*EN1RYD/BOLTZMANN/vf.temp);
2213  else
2214  reduced_b_max1 = vf.reduced_b_max;
2215 
2216  // Return strictly zero when there is no allowed range: falling
2217  // through with MAX2 generates junk values due to rounding error
2218  if (reduced_b_max1 <= reduced_b_min)
2219  return 0.;
2220 
2221  reduced_b_max1 = MAX2( reduced_b_max1, reduced_b_min );
2222  ASSERT( reduced_b_max1 > 0. );
2223 
2224  double alphamin = 1.5*vf.ColliderCharge/(reduced_vel * reduced_b_max1);
2225  ASSERT( alphamin > 0. );
2226  ASSERT( alphamax > 0. );
2227 
2228  // set up the integrator.
2229 
2230  double alphamin_int = MAX2( alphamin, 1.e-30 );
2231  double alphamax_int = MAX2( alphamax, 1.e-20 );
2232 
2233  double CSIntegral = 0.;
2234 
2235  if (0)
2236  {
2237  my_Integrand_VF01_alpha<P> func(vf.n, vf.l, vf.lp, alphamin);
2238  CSIntegral = CSIntegral_QG32(func,alphamin_int,alphamax_int);
2239  }
2240  else
2241  {
2242  my_Integrand_VF01_beta<P> funcb(vf.n, vf.l, vf.lp, alphamin);
2243  if (1)
2244  {
2245  CSIntegral = CSIntegral_Romberg(ipISO,funcb,alphamin_int,alphamax_int,1e-4);
2246  }
2247  else
2248  {
2249  double epsabs=0., epsrel=1e-5, abserr, qqq;
2250  const long limit=25, lenw=4*limit;
2251  long neval, ier,last, iwork[limit];
2252  double work[lenw];
2253  double lalphamin_int = log(alphamin_int),
2254  lalphamax_int = log(alphamax_int);
2255  dqags_(funcb,&lalphamin_int,&lalphamax_int,&epsabs,&epsrel,&qqq,
2256  &abserr,&neval,&ier,&limit,&lenw,&last,iwork,work);
2257  CSIntegral = qqq;
2258  }
2259 
2260  if (0)
2261  {
2262  double qqq;
2263  my_Integrand_VF01_alpha<P> func(vf.n, vf.l, vf.lp, alphamin);
2264  const char *ccc;
2265  if (0)
2266  {
2267  ccc = "G32";
2268  qqq = CSIntegral_QG32(func,alphamin_int,alphamax_int);
2269  }
2270  else if (0)
2271  {
2272  ccc = "Unirom";
2273  qqq = CSIntegral_Romberg_alpha(ipISO,func,alphamin_int,alphamax_int,1e-5);
2274  }
2275  else
2276  {
2277  ccc = "QAGS";
2278  double epsabs=0., epsrel=1e-5, abserr;
2279  const long limit=25, lenw=4*limit;
2280  long neval, ier,last, iwork[limit];
2281  double work[lenw];
2282  double lalphamin_int = log(alphamin_int),
2283  lalphamax_int = log(alphamax_int);
2284  dqags_(funcb,&lalphamin_int,&lalphamax_int,&epsabs,&epsrel,&qqq,
2285  &abserr,&neval,&ier,&limit,&lenw,&last,iwork,work);
2286  fprintf(ioQQQ,"Check VF QAGS1 err=%g neval=%ld ier=%ld\n",abserr,neval,ier);
2287  }
2288  if (CSIntegral > 0. || qqq > 0.)
2289  {
2290  fprintf(ioQQQ,"Check VF[%ld %ld->%ld %g %g]: %s %g var %g err=%g\n",
2291  vf.n,vf.l,vf.lp,E_Proj_Ryd,alphamin,
2292  ccc,qqq,CSIntegral,(qqq-CSIntegral)/SDIV(CSIntegral));
2293  }
2294  }
2295  }
2296 
2297 
2298  /* Factors outside integral */
2299  double ConstantFactors = 4.5*PI*POW2(vf.ColliderCharge*vf.aveRadius/reduced_vel);
2300 
2301  /* Calculate cross section */
2302  double cross_section = ConstantFactors * CSIntegral;
2303 
2304  /* convert to collision strength, cross section is already in cm^2 */
2305  double coll_str = ConvCrossSect2CollStr( cross_section, vf.gLo, E_Proj_Ryd, vf.reduced_mass );
2306 
2307  coll_str = MAX2( SMALLFLOAT, coll_str);
2308 
2309  return coll_str;
2310 }
2311 
my_Integrand_VF01_E::operator()
double operator()(double EOverKT) const
Definition: helike_cs.cpp:208
my_Integrand_VF01_E::lp
long lp
Definition: helike_cs.cpp:64
my_Integrand_S62
Definition: helike_cs.cpp:998
D_fp::operator()
virtual double operator()(double) const =0
my_Integrand_VF01_E::gLo
long gLo
Definition: helike_cs.cpp:64
integrate::Romberg::evals
long evals() const
Definition: integrate.h:196
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:765
my_Integrand_S62::deltaE
double deltaE
Definition: helike_cs.cpp:1000
avx_ptr
Definition: vectorize.h:170
TransitionProxy::EnergyErg
realnum EnergyErg() const
Definition: transition.h:90
CSTemp
vector< double > CSTemp
Definition: helike_cs.cpp:25
t_iso_sp::st
qList st
Definition: iso.h:482
t_isoCTRL::lgCS_PSClassic
bool lgCS_PSClassic[NISO]
Definition: iso.h:408
avx_ptr::ptr0
T * ptr0()
Definition: vectorize.h:240
exp10
double exp10(double x)
Definition: cddefines.h:1368
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:65
my_Integrand_VF01_E_log
Definition: helike_cs.cpp:221
my_Integrand_VF01_E::RMSv
double RMSv
Definition: helike_cs.cpp:68
get_ptr
T * get_ptr(T *v)
Definition: cddefines.h:1109
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
multi_arr::alloc
void alloc()
Definition: container_classes.h:1113
BIGDOUBLE
const double BIGDOUBLE
Definition: cpu.h:249
my_Integrand_S62::my_Integrand_S62
my_Integrand_S62(double deltaE, double osc_strength, double temp, double I_energy_eV)
Definition: helike_cs.cpp:1002
my_Integrand_VF01_E::l
long l
Definition: helike_cs.cpp:64
CS_l_mixing_S62
STATIC double CS_l_mixing_S62(double deltaE_eV, double IP_Ryd_ground, long gLo, double Aul, long nelem, long Collider, double temp)
Definition: helike_cs.cpp:1010
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:246
iso_ctrl
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
e1
double e1(double x)
Definition: thirdparty.cpp:3258
my_Integrand_VF01_E
Definition: helike_cs.cpp:38
ipHe2p3P0
const int ipHe2p3P0
Definition: iso.h:48
t_isoCTRL::lgCS_therm_ave
bool lgCS_therm_ave[NISO]
Definition: iso.h:408
CS_l_mixing_VF01
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
SixJFull
double SixJFull(long j1, long j2, long j3, long j4, long j5, long j6)
Definition: thirdparty.cpp:3558
integrate::Romberg::update
void update(double eps)
Definition: integrate.h:203
t_iso_sp::nCollapsed_max
long int nCollapsed_max
Definition: iso.h:518
my_Integrand_VF01_E::Collider
long Collider
Definition: helike_cs.cpp:66
CSIntegral_Romberg
STATIC double CSIntegral_Romberg(long ipISO, const my_Integrand_VF01_beta< P > &func, double alphamin_int, double alphamax_int, double eps)
Definition: helike_cs.cpp:2108
GetHelikeCollisionStrength
realnum GetHelikeCollisionStrength(long nelem, long Collider, long nHi, long lHi, long sHi, long jHi, long gHi, double IP_Ryd_Hi, long nLo, long lLo, long sLo, long jLo, long gLo, double IP_Ryd_Lo, double Aul, double tauLo, double EnerWN, double EnerErg, const char **where)
Definition: helike_cs.cpp:434
D_fp
Definition: thirdparty_quadpack.h:26
phycon
t_phycon phycon
Definition: phycon.cpp:6
fp_equal_tol
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
Definition: cddefines.h:904
HeCollidSetup
void HeCollidSetup(void)
Definition: helike_cs.cpp:241
my_Integrand_VF01_E::reduced_mass
double reduced_mass
Definition: helike_cs.cpp:68
t_isoCTRL::lgCS_None
bool lgCS_None[NISO]
Definition: iso.h:408
CS_l_mixing_PS64
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_isoCTRL::lgColl_l_mixing
bool lgColl_l_mixing[NISO]
Definition: iso.h:359
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
my_Integrand_S62::osc_strength
double osc_strength
Definition: helike_cs.cpp:1000
helike_cs.h
sjs
double sjs(long j1, long j2, long j3, long l1, long l2, long l3)
Definition: thirdparty.cpp:3494
hydro_vs_rates.h
integrate::Romberg::sum
double sum() const
Definition: integrate.h:188
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:481
MIN2
#define MIN2(a, b)
Definition: cddefines.h:803
HeCS
multi_arr< realnum, 3 > HeCS
Definition: helike_cs.cpp:27
t_isoCTRL::lgCS_Vrinceanu
bool lgCS_Vrinceanu[NISO]
Definition: iso.h:408
HeCSInterp
realnum HeCSInterp(long nelem, long ipHi, long ipLo, long Collider)
Definition: helike_cs.cpp:390
my_Integrand_VF01_E::ipISO
long ipISO
Definition: helike_cs.cpp:64
dense
t_dense dense
Definition: global.cpp:15
bessel_k0_k1
void bessel_k0_k1(double x, double *k0val, double *k1val)
Definition: thirdparty.cpp:2792
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
vexp
void vexp(const double x[], double y[], long nlo, long nhi)
Definition: vectorize_exp.cpp:283
t_isoCTRL::SmallA
realnum SmallA
Definition: iso.h:391
reduced_amu
double reduced_amu(long nelem, long Collider)
Definition: helike_cs.cpp:47
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
my_Integrand_VF01_E_log::operator()
double operator()(double y) const
Definition: helike_cs.cpp:231
my_Integrand_VF01_E::IP_Ryd_Hi
double IP_Ryd_Hi
Definition: helike_cs.cpp:65
trace
t_trace trace
Definition: trace.cpp:5
lines_service.h
my_Integrand_VF01_E::n
long n
Definition: helike_cs.cpp:64
t_isoCTRL::lgColl_excite
bool lgColl_excite[NISO]
Definition: iso.h:362
my_Integrand_VF01_E::aveRadius
double aveRadius
Definition: helike_cs.cpp:68
my_Integrand_S62::temp
double temp
Definition: helike_cs.cpp:1000
colliders
ColliderList colliders(dense)
CS_l_mixing_VOS12
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
t_isoCTRL::nCS_new
int nCS_new[NISO]
Definition: iso.h:419
CS_VS80
double CS_VS80(long nHi, long gHi, double IP_Ryd_Hi, long nLo, long gLo, double IP_Ryd_Lo, double Aul, long nelem, long Collider, double temp)
Definition: hydro_vs_rates.cpp:49
my_Integrand_VF01_E_log::my_Integrand_VF01_E_log
my_Integrand_VF01_E_log(long ipISO_1, long nelem_1, long n_1, long l_1, long lp_1, long s, long gLo_1, double tauLo_1, double IP_Ryd_Hi_1, double IP_Ryd_Lo_1, long Collider_1, double temp_1)
Definition: helike_cs.cpp:226
TransitionProxy::EnergyWN
realnum & EnergyWN() const
Definition: transition.h:477
VecIntegrator::sum
double sum(double min, double max) const
Definition: integrate.h:68
POW2
#define POW2
Definition: cddefines.h:979
collision_strength_VF01
STATIC double collision_strength_VF01(long ipISO, double velOrEner, const my_Integrand_VF01_E< P > &vf)
Definition: helike_cs.cpp:2167
multi_arr< realnum, 3 >
STATIC
#define STATIC
Definition: cddefines.h:118
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:447
t_isoCTRL::lgCS_PS64
bool lgCS_PS64[NISO]
Definition: iso.h:408
realnum
float realnum
Definition: cddefines.h:124
dqags_
void dqags_(const D_fp &f, const double *a, const double *b, const double *epsabs, const double *epsrel, double *result, double *abserr, long *neval, long *ier, const long *limit, const long *lenw, long *last, long *iwork, double *work)
Definition: thirdparty_quadpack.cpp:3124
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:168
my_Integrand_VF01_E_log::data
my_Integrand_VF01_E< P > data
Definition: helike_cs.cpp:224
max
long max(int a, long b)
Definition: cddefines.h:817
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:482
powi
double powi(double, long int)
Definition: service.cpp:594
my_Integrand_VF01_E::IP_Ryd_Lo
double IP_Ryd_Lo
Definition: helike_cs.cpp:65
min
long min(int a, long b)
Definition: cddefines.h:762
strchr_s
const char * strchr_s(const char *s, int c)
Definition: cddefines.h:1310
S62BesselInvert
STATIC double S62BesselInvert(double zOverB2)
Definition: helike_cs.cpp:1054
my_Integrand_VF01_E::ColliderMass
double ColliderMass
Definition: helike_cs.cpp:68
Integrator::sum
double sum(double min, double max) const
Definition: integrate.h:44
my_Integrand_S62::I_energy_eV
double I_energy_eV
Definition: helike_cs.cpp:1000
my_Integrand_VF01_E::temp
double temp
Definition: helike_cs.cpp:67
t_iso_sp::QuantumNumbers2Index
multi_arr< long, 3 > QuantumNumbers2Index
Definition: iso.h:490
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:473
t_dense::AtomicWeight
realnum AtomicWeight[LIMELM]
Definition: dense.h:80
factorial
double factorial(long n)
Definition: thirdparty.cpp:356
CS_l_mixing_VOS12QM
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
t_isoCTRL::lgCS_Seaton
bool lgCS_Seaton[NISO]
Definition: iso.h:408
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:613
HydroEinstA
double HydroEinstA(long int n1, long int n2)
Definition: hydroeinsta.cpp:13
my_Integrand_VF01_E::ColliderCharge
double ColliderCharge
Definition: helike_cs.cpp:67
t_isoCTRL::lgCS_VOS12
bool lgCS_VOS12[NISO]
Definition: iso.h:408
multi_arr::reserve
void reserve(size_type i1)
Definition: container_classes.h:1077
ColliderList::list
vector< t_collider > list
Definition: collision.h:41
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:64
LIMELM
const int LIMELM
Definition: cddefines.h:308
pow2
T pow2(T a)
Definition: cddefines.h:981
my_Integrand_S62::operator()
void operator()(const double proj_energy_overKT[], double res[], long n) const
Definition: helike_cs.cpp:1123
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
powpq
double powpq(double x, int p, int q)
Definition: service.cpp:630
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:350
my_Integrand_VF01_E::nelem
long nelem
Definition: helike_cs.cpp:64
ColliderCharge
static const double ColliderCharge[4]
Definition: helike_cs.cpp:45
t_isoCTRL::lgCS_VOS12QM
bool lgCS_VOS12QM[NISO]
Definition: iso.h:408
t_dense::eden
double eden
Definition: dense.h:201
cross_section
STATIC double cross_section(double EgammaRyd, double EthRyd, long nelem, long n, long l, long s)
Definition: helike_recom.cpp:66
ConvCrossSect2CollStr
double ConvCrossSect2CollStr(double CrsSectCM2, double gLo, double E_ProjectileRyd, double reduced_mass_grams)
Definition: lines_service.cpp:686
rec6j
void rec6j(double *sixcof, double l2, double l3, double l4, double l5, double l6, double *l1min, double *l1max, double *lmatch, long ndim, long *ier)
Definition: thirdparty.cpp:3608
MAX2
#define MAX2(a, b)
Definition: cddefines.h:824
my_Integrand_VF01_E::reduced_b_max
double reduced_b_max
Definition: helike_cs.cpp:68
ellpk
double ellpk(double x)
Definition: thirdparty.cpp:1676
fprintf
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
chLine_LENGTH
#define chLine_LENGTH
my_Integrand_VF01_E::my_Integrand_VF01_E
my_Integrand_VF01_E(long ipISO_1, long nelem_1, long n_1, long l_1, long lp_1, long s, long gLo_1, double tauLo_1, double IP_Ryd_Hi_1, double IP_Ryd_Lo_1, long Collider_1, double temp_1)
Definition: helike_cs.cpp:70
SDIV
sys_float SDIV(sys_float x)
Definition: cddefines.h:1002
t_phycon::alogte
double alogte
Definition: phycon.h:92
S
#define S(I_, J_)
Definition: optimize_subplx.cpp:1826
t_isoCTRL::lgCS_B72
bool lgCS_B72[NISO]
Definition: iso.h:408
t_iso_sp::numLevels_max
long int numLevels_max
Definition: iso.h:524
read_whole_line
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
CSIntegral_Romberg_alpha
STATIC double CSIntegral_Romberg_alpha(long ipISO, const my_Integrand_VF01_alpha< P > &func, double alphamin_int, double alphamax_int, double eps)
Definition: helike_cs.cpp:2138
HeCSTableInterp
STATIC realnum HeCSTableInterp(long nelem, long Collider, long nHi, long lHi, long sHi, long jHi, long nLo, long lLo, long sLo, long jLo)
Definition: helike_cs.cpp:923
t_phycon::sqrte
double sqrte
Definition: phycon.h:58
fixit
#define fixit(a)
Definition: cddefines.h:417
gegenbauer
double gegenbauer(long n, double al, double x)
Definition: thirdparty.cpp:3464
integrate::Romberg::error
double error() const
Definition: integrate.h:192
POW3
#define POW3
Definition: cddefines.h:986
t_phycon::te
double te
Definition: phycon.h:21
frac
double frac(double d)
Definition: thirdparty.cpp:3597
CS_l_mixing_PS64_expI
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
CSIntegral_QG32
STATIC double CSIntegral_QG32(const my_Integrand_VF01_alpha< P > &func, double alphamin_int, double alphamax_int)
Definition: helike_cs.cpp:2093
my_Integrand_VF01_E::tauLo
double tauLo
Definition: helike_cs.cpp:65
CS_l_mixing
double CS_l_mixing(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:1540
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:690
thirdparty_quadpack.h
COLLISMAGIC
#define COLLISMAGIC
Definition: helike.h:24
t_isoCTRL::lgCS_PSdeg
bool lgCS_PSdeg[NISO]
Definition: iso.h:408
t_isoCTRL::lgCS_Vriens
bool lgCS_Vriens[NISO]
Definition: iso.h:408
FFmtRead
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:381
Generated by   doxygen 1.8.5