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