cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
iso_collide.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 #include "cddefines.h"
4 #include "atmdat_adfa.h"
5 #include "conv.h"
6 #include "heavy.h"
7 #include "helike_cs.h"
8 #include "hydroeinsta.h"
9 #include "hydrogenic.h"
10 #include "hydro_vs_rates.h"
11 #include "ionbal.h"
12 #include "iso.h"
13 #include "opacity.h"
14 #include "phycon.h"
15 #include "rfield.h"
16 #include "secondaries.h"
17 #include "trace.h"
18 #include "freebound.h"
19 #include "dense.h"
20 #include "lines_service.h"
21 #include "vectorize.h"
22 
23 STATIC double iso_get_collision_strength( long ipISO, long nelem, long ipCollider, long ipHi, long ipLo );
24 #if 0
25 STATIC double iso_get_collision_strength_collapsed_to_collapsed( long ipISO, long nelem, long ipCollider,
26  long nHi, double IP_Ryd_Hi,
27  long nLo, double IP_Ryd_Lo,
28  double Aul, double tauLo, double EnerWN, double EnerErg );
29 #endif
30 STATIC double iso_get_collision_strength_collapsed_to_resolved( long ipISO, long nelem, long ipCollider,
31  long nHi, double IP_Ryd_Hi,
32  long nLo, long lLo, long sLo, long jLo, long gLo, double IP_Ryd_Lo,
33  double Aul, double tauLo, double EnerWN, double EnerErg );
34 STATIC double iso_get_collision_strength_resolved( long ipISO, long nelem, long ipCollider,
35  long nHi, long lHi, long sHi, long jHi, long gHi, double IP_Ryd_Hi,
36  long nLo, long lLo, long sLo, long jLo, long gLo, double IP_Ryd_Lo,
37  double Aul, double tauLo, double EnerWN, double EnerErg );
38 
39 STATIC double iso_get_collision_strength_collapsed_to_collapsed_fast( long ipISO, long nelem, long ipCollider,
40  long nHi, double IP_Ryd_Hi,
41  long nLo, double IP_Ryd_Lo,
42  double tauLo, double EnerWN, double EnerErg );
43 
44 void iso_collisional_ionization( long ipISO, long nelem )
45 {
46  ASSERT( ipISO < NISO );
47 
48  DEBUG_ENTRY( "iso_collisional_ionization()" );
49 
50  t_iso_sp* sp = &iso_sp[ipISO][nelem];
51 
52  /* collisional ionization from ground */
53  sp->fb[0].ColIoniz = iso_ctrl.lgColl_ionize[ipISO] *
54  t_ADfA::Inst().coll_ion_wrapper( nelem, nelem-ipISO, phycon.te );
55 
56  iso_put_error(ipISO,nelem,sp->numLevels_max,0,IPCOLLIS,0.20f,0.20f);
57 
58  for( long ipHi=1; ipHi<sp->numLevels_max; ipHi++ )
59  {
60  if( nelem == ipISO )
61  {
62  /* use routine from Vriens and Smeets (1981). */
63  /* >>refer iso neutral col.ion. Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940 */
64  sp->fb[ipHi].ColIoniz = hydro_vs_ioniz( sp->fb[ipHi].xIsoLevNIonRyd, phycon.te );
65  }
66  else
67  {
68  /* ions */
69  /* use hydrogenic ionization rates for ions
70  * >>refer iso ions col.ion. Allen 1973, Astro. Quan. for low Te.
71  * >>refer iso ions col.ion. Sampson and Zhang 1988, ApJ, 335, 516 for High Te.
72  * */
73  sp->fb[ipHi].ColIoniz =
74  Hion_coll_ioniz_ratecoef( ipISO, nelem, N_(ipHi), sp->fb[ipHi].xIsoLevNIonRyd, phycon.te );
75  }
76 
77  // iso_ctrl.lgColl_ionize is option to turn off collisions, "atom XX-like collis off" command
78  sp->fb[ipHi].ColIoniz *= iso_ctrl.lgColl_ionize[ipISO];
79 
80  iso_put_error(ipISO,nelem,sp->numLevels_max,ipHi,IPCOLLIS,0.20f,0.20f);
81  }
82 
83  /* Here we arbitrarily scale the highest level ionization to account for the fact
84  * that, if the atom is not full size, this level should be interacting with higher
85  * levels and not just the continuum. We did add on collisional excitation terms instead
86  * but this caused a problem at low temperatures because the collisional ionization was
87  * a sum of terms with different Boltzmann factors, while PopLTE had just one Boltzmann
88  * factor. The result was a collisional recombination that had residual exponentials of
89  * the form exp(x/kT), which blew up at small T. */
90  if( !sp->lgLevelsLowered && iso_ctrl.lgTopoff[ipISO] )
91  {
92  sp->fb[sp->numLevels_max-1].ColIoniz *= 100.;
93  /*>>chng 16 nov 13, we had tried the following which is suggested by the atomic physics.
94  * To come into LTE with the continuum the population of the highest level has to be
95  * dominated by collisional ionization / three body recombination. The power law
96  * charge dependence is that of the lifetime of the highest level. This caused
97  * an abort due to non-conservation of energy in auto / blr_n13_p22_Z20.out and
98  * slow / feii_blr_n13_p22_Z20.out. See ticket #374
99  */
100  //sp->fb[sp->numLevels_max-1].ColIoniz *= (100. * powi(nelem-ipISO+1,4) );
101  }
102 
103  return;
104 }
105 
106 void iso_suprathermal( long ipISO, long nelem )
107 {
108  DEBUG_ENTRY( "iso_suprathermal()" );
109 
110  /* check that we were called with valid parameters */
111  ASSERT( ipISO < NISO );
112  ASSERT( nelem >= ipISO );
113  ASSERT( nelem < LIMELM );
114 
115  t_iso_sp* sp = &iso_sp[ipISO][nelem];
116 
117  /***********************************************************************
118  * *
119  * get secondary excitation by suprathermal electrons *
120  * *
121  ***********************************************************************/
122 
123  for( long i=1; i < sp->numLevels_max; i++ )
124  {
125  if( sp->trans(i,0).ipCont() > 0 )
126  {
127  /* get secondaries for all iso lines by scaling LyA
128  * excitation by ratio of cross section (oscillator strength/energy)
129  * Born approximation or plane-wave approximation based on
130  *>>refer HI excitation Shemansky, D.E., et al., 1985, ApJ, 296, 774 */
132  (sp->trans(i,0).Emis().gf()/
133  sp->trans(i,0).EnergyWN()) /
134  (iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,0).Emis().gf()/
136  iso_ctrl.lgColl_excite[ipISO];;
137  }
138  else
139  sp->trans(i,0).Coll().rate_lu_nontherm_set() = 0.;
140  }
141 
142  return;
143 }
144 
145 /*============================*/
146 /* evaluate collisional rates */
147 void iso_collide( long ipISO, long nelem )
148 {
149  DEBUG_ENTRY( "iso_collide()" );
150 
151  /* this array stores the last temperature at which collision data were evaluated for
152  * each species of the isoelectronic sequence. */
153  static double TeUsed[NISO][LIMELM]={
154  {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0},
155  {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0} };
156 
157  /* check that we were called with valid parameters */
158  ASSERT( ipISO < NISO );
159  ASSERT( nelem >= ipISO );
160  ASSERT( nelem < LIMELM );
161 
162  t_iso_sp* sp = &iso_sp[ipISO][nelem];
163 
164  /* skip most of this routine if temperature has not changed,
165  * the check on conv.nTotalIoniz is to make sure that we redo this
166  * on the very first call in a grid calc - it is 0 on the first call */
167  if( fp_equal( TeUsed[ipISO][nelem], phycon.te ) && conv.nTotalIoniz )
168  {
169  ASSERT( sp->trans( iso_ctrl.nLyaLevel[ipISO], 0 ).Coll().ColUL( colliders ) >= 0. );
170 
171  if( trace.lgTrace && (trace.lgHBug||trace.lgHeBug) )
172  {
173  fprintf( ioQQQ,
174  " iso_collide called %s nelem %li - no reeval Boltz fac, LTE dens\n",
175  iso_ctrl.chISO[ipISO], nelem );
176  }
177  }
178  else
179  {
180  TeUsed[ipISO][nelem] = phycon.te;
181 
182  if( trace.lgTrace && (trace.lgHBug||trace.lgHeBug) )
183  {
184  fprintf( ioQQQ,
185  " iso_collide called %s nelem %li - will reeval Boltz fac, LTE dens\n",
186  iso_ctrl.chISO[ipISO], nelem );
187  }
188 
189  /**********************************************************
190  * *
191  * Boltzmann factors for all levels, and *
192  * collisional ionization and excitation *
193  * *
194  **********************************************************/
195 
196  /* HION_LTE_POP is planck^2 / (2 pi m_e k ), raised to 3/2 next */
197  double factor = HION_LTE_POP*dense.AtomicWeight[nelem]/
198  (dense.AtomicWeight[nelem]+ELECTRON_MASS/ATOMIC_MASS_UNIT);
199 
200  /* term in () is stat weight of electron * ion */
201  double ConvLTEPOP = powpq(factor,3,2)/(2.*iso_ctrl.stat_ion[ipISO])/phycon.te32;
202 
203  sp->lgPopLTE_OK = true;
204 
205  // this is the maximum value of iso.PopLTE (units cm^3) that corresponds
206  // to the minimum positive density values. A smaller density will be
207  // regarded as zero, and the product PopLTE*n_e*n_Z+ will also be zero.
208  const double MAX_POP_LTE = (MAX_DENSITY/dense.density_low_limit/dense.density_low_limit);
209 
211  /* this Boltzmann factor is exp( +ioniz energy / Te ) */
212  for( long ipLo=0; ipLo < sp->numLevels_max; ipLo++ )
213  arg[ipLo] = -sp->st[ipLo].energy().Ryd()/phycon.te_ryd;
214  vexp( arg.ptr0(), sp->st.Boltzmann(), 0, sp->numLevels_max );
215  for( long ipLo=0; ipLo < sp->numLevels_max; ipLo++ )
216  arg[ipLo] = -sp->fb[ipLo].xIsoLevNIonRyd/phycon.te_ryd;
217  vexp( arg.ptr0(), sp->st.ConBoltz(), 0, sp->numLevels_max );
218  /* fully define Boltzmann factors to continuum for model levels */
219  for( long ipLo=0; ipLo < sp->numLevels_max; ipLo++ )
220  {
221  /***************************************
222  * *
223  * LTE abundances for all levels *
224  * *
225  ***************************************/
226 
227  double safe_limit = SMALLDOUBLE*sp->st[ipLo].g()*max(ConvLTEPOP,1.);
228  if( sp->st[ipLo].ConBoltz() > safe_limit )
229  {
230  /* LTE population of given level. */
231  sp->fb[ipLo].PopLTE =
232  sp->st[ipLo].g() / sp->st[ipLo].ConBoltz() * ConvLTEPOP;
233  ASSERT( sp->fb[ipLo].PopLTE < BIGDOUBLE );
234  }
235  else
236  {
237  sp->fb[ipLo].PopLTE = 0.;
238  }
239 
240  sp->fb[ipLo].PopLTE = MIN2( sp->fb[ipLo].PopLTE, MAX_POP_LTE );
241 
242  /* now check for any zeros - if present then matrix cannot be used */
243  if( sp->fb[ipLo].PopLTE <= 0. )
244  {
245  sp->lgPopLTE_OK = false;
246  }
247  }
248 
249  iso_collisional_ionization( ipISO, nelem );
250 
251  /***********************************************************
252  * *
253  * collisional deexcitation for all lines in iso sequence *
254  * *
255  ***********************************************************/
256 
257  if( iso_ctrl.lgColl_excite[ipISO] )
258  {
259  vector<double> ratefac(ipALPHA+1);
260  for( long ipCollider = ipELECTRON; ipCollider <= ipALPHA; ipCollider++ )
261  {
262  double reduced_mass_collider_system = dense.AtomicWeight[nelem]*colliders.list[ipCollider].mass_amu/
263  (dense.AtomicWeight[nelem]+colliders.list[ipCollider].mass_amu)*ATOMIC_MASS_UNIT;
264  ratefac[ipCollider] = powpq(ELECTRON_MASS/reduced_mass_collider_system,3,2) * COLL_CONST/phycon.sqrte;
265  }
266 
267  for( long ipHi=1; ipHi<sp->numLevels_max; ipHi++ )
268  {
269  double overHig = 1./(double)sp->st[ipHi].g();
270  for( long ipLo=0; ipLo < ipHi; ipLo++ )
271  {
272  for( long ipCollider = ipELECTRON; ipCollider <= ipALPHA; ipCollider++ )
273  {
274  double cs_temp = iso_get_collision_strength( ipISO, nelem, ipCollider, ipHi, ipLo );
275 
276  // store electron collision strength in generic collision strength
277  if( ipCollider == ipELECTRON )
278  sp->trans(ipHi,ipLo).Coll().col_str() = (realnum) cs_temp;
279 
280  double rateCoef = cs_temp * ratefac[ipCollider] * overHig;
281 
282  sp->trans(ipHi,ipLo).Coll().rate_coef_ul_set()[ipCollider] =
283  (realnum) rateCoef;
284  }
285 
286  /* check for sanity */
287  ASSERT( sp->trans(ipHi,ipLo).Coll().rate_coef_ul()[ipELECTRON] >= 0. );
288 
289  if( N_(ipHi) <= 5 && N_(ipLo) <= 2 )
290  iso_put_error( ipISO, nelem, ipHi, ipLo, IPCOLLIS, 0.10f, 0.30f );
291  else
292  iso_put_error( ipISO, nelem, ipHi, ipLo, IPCOLLIS, 0.20f, 0.30f );
293  }
294  }
295  }
296 
297  if( (trace.lgTrace && trace.lgIsoTraceFull[ipISO]) && (nelem == trace.ipIsoTrace[ipISO]) )
298  {
299  fprintf( ioQQQ, " iso_collide: %s Z=%li de-excitation rates coefficients\n", iso_ctrl.chISO[ipISO], nelem + 1 );
300  long upper_limit = sp->numLevels_local;
301  for( long ipHi=1; ipHi < upper_limit; ipHi++ )
302  {
303  fprintf( ioQQQ, " %li\t", ipHi );
304  for( long ipLo=0; ipLo < ipHi; ipLo++ )
305  {
306  fprintf( ioQQQ,PrintEfmt("%10.3e",
307  sp->trans(ipHi,ipLo).Coll().ColUL( colliders ) / dense.EdenHCorr ));
308  }
309  fprintf( ioQQQ, "\n" );
310  }
311 
312  fprintf( ioQQQ, " iso_collide: %s Z=%li collisional ionization coefficients\n", iso_ctrl.chISO[ipISO], nelem + 1 );
313  for( long ipHi=0; ipHi < upper_limit; ipHi++ )
314  {
315  fprintf( ioQQQ,PrintEfmt("%10.3e", sp->fb[ipHi].ColIoniz ));
316  }
317  fprintf( ioQQQ, "\n" );
318 
319  fprintf( ioQQQ, " iso_collide: %s Z=%li continuum boltzmann factor\n", iso_ctrl.chISO[ipISO], nelem + 1 );
320  for( long ipHi=0; ipHi < upper_limit; ipHi++ )
321  {
322  fprintf( ioQQQ,PrintEfmt("%10.3e", sp->st[ipHi].ConBoltz() ));
323  }
324  fprintf( ioQQQ, "\n" );
325 
326  fprintf( ioQQQ, " iso_collide: %s Z=%li continuum boltzmann factor\n", iso_ctrl.chISO[ipISO], nelem + 1 );
327  for( long ipHi=0; ipHi < upper_limit; ipHi++ )
328  {
329  fprintf( ioQQQ,PrintEfmt("%10.3e", sp->fb[ipHi].PopLTE ));
330  }
331  fprintf( ioQQQ, "\n" );
332  }
333 
334  /* the case b hummer and storey option,
335  * this kills collisional excitation and ionization from n=1 and n=2 */
337  {
338  for( long ipLo=0; ipLo<sp->numLevels_max-1; ipLo++ )
339  {
340  if( N_(ipLo)>=3 )
341  break;
342 
343  sp->fb[ipLo].ColIoniz = 0.;
344 
345  for( long ipHi=ipLo+1; ipHi<sp->numLevels_max; ipHi++ )
346  {
347  /* don't disable 2-2 collisions */
348  if( N_(ipLo)==2 && N_(ipHi)==2 )
349  continue;
350 
351  sp->trans(ipHi,ipLo).Coll().col_str() = 0.;
352  for( long k=0; k<ipNCOLLIDER; ++k )
353  sp->trans(ipHi,ipLo).Coll().rate_coef_ul_set()[k] = 0.;
354  }
355  }
356  }
357  }
358 
359  iso_suprathermal( ipISO, nelem );
360 
361  /* this must be reevaluated every time since eden can change when Te does not */
362  /* save into main array - collisional ionization by thermal electrons */
363  ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0] =
364  sp->fb[0].ColIoniz*dense.EdenHCorr;
365 
366  /* cooling due to collisional ionization, which only includes thermal electrons */
367  ionbal.CollIonRate_Ground[nelem][nelem-ipISO][1] =
368  ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0]*
369  rfield.anu(Heavy.ipHeavy[nelem][nelem-ipISO]-1)*EN1RYD;
370 
371  return;
372 }
373 
374 STATIC double iso_get_collision_strength( long ipISO, long nelem, long ipCollider, long ipHi, long ipLo )
375 {
376  DEBUG_ENTRY( "iso_get_collision_strength()" );
377 
378  t_iso_sp* sp = &iso_sp[ipISO][nelem];
379  long nHi = sp->st[ipHi].n();
380  long lHi = sp->st[ipHi].l();
381  long sHi = sp->st[ipHi].S();
382  long jHi = sp->st[ipHi].j();
383  long gHi = sp->st[ipHi].g();
384  double IP_Ryd_Hi = sp->fb[ipHi].xIsoLevNIonRyd;
385  long nLo = sp->st[ipLo].n();
386  long lLo = sp->st[ipLo].l();
387  long sLo = sp->st[ipLo].S();
388  long jLo = sp->st[ipLo].j();
389  long gLo = sp->st[ipLo].g();
390  double IP_Ryd_Lo = sp->fb[ipLo].xIsoLevNIonRyd;
391  double Aul = sp->trans(ipHi,ipLo).Emis().Aul();
392  // collisions are from high to low level, then initial level lifetime is from higher level
393  double tauLo = sp->st[ipHi].lifetime();
394  double EnerWN = sp->trans(ipHi,ipLo).EnergyWN();
395  double EnerErg = sp->trans(ipHi,ipLo).EnergyErg();
396 
397  /* collision strength we derive */
398  double cs = 0.;
399 
400  if( nHi==nLo && !iso_ctrl.lgColl_l_mixing[ipISO] )
401  cs = 0.;
402  else if( nHi-nLo > 2 && ipCollider > ipELECTRON )
403  cs = 0.;
404  else if( nHi > sp->n_HighestResolved_max && nLo > sp->n_HighestResolved_max )
405  {
406  ASSERT( ipLo >= sp->numLevels_max - sp->nCollapsed_max );
407 #if 1
408  cs = iso_get_collision_strength_collapsed_to_collapsed_fast( ipISO, nelem, ipCollider,
409  nHi, IP_Ryd_Hi,
410  nLo, IP_Ryd_Lo,
411  tauLo, EnerWN, EnerErg );
412 #else
413  cs = iso_get_collision_strength_collapsed_to_collapsed( ipISO, nelem, ipCollider,
414  nHi, IP_Ryd_Hi,
415  nLo, IP_Ryd_Lo,
416  Aul, tauLo, EnerWN, EnerErg );
417 #endif
418  }
419  else if( nHi > sp->n_HighestResolved_max && nLo <= sp->n_HighestResolved_max )
420  {
421  ASSERT( ipLo < sp->numLevels_max - sp->nCollapsed_max );
422  cs = iso_get_collision_strength_collapsed_to_resolved( ipISO, nelem, ipCollider,
423  nHi, IP_Ryd_Hi,
424  nLo, lLo, sLo, jLo, gLo, IP_Ryd_Lo,
425  Aul, tauLo, EnerWN, EnerErg );
426  }
427  else
428  {
429  cs = iso_get_collision_strength_resolved( ipISO, nelem, ipCollider,
430  nHi, lHi, sHi, jHi, gHi, IP_Ryd_Hi,
431  nLo, lLo, sLo, jLo, gLo, IP_Ryd_Lo,
432  Aul, tauLo, EnerWN, EnerErg );
433  }
434 
435  cs *= iso_ctrl.lgColl_excite[ipISO];
436 
437  if( opac.lgCaseB_HummerStorey && nHi==nLo+1 && nHi<=sp->n_HighestResolved_max && abs(lHi-lLo)==1 )
438  {
439  double Aul_hydro = HydroEinstA( nLo, nHi );
440  cs *= ( sp->trans(ipHi,ipLo).Emis().gf() / gLo ) /
441  ( GetGF(Aul_hydro, EnerWN, 2.*nHi*nHi)/(2.*nLo*nLo) );
442  }
443  if (0)
444  {
445  //printing some of the rates for testing
446  double rate_test;
447  double oHi=1./(double)gHi;
448  //double ogLo=1./(double)gLo;
449  if ( nelem == 1 && nHi == nLo && nHi == 30 && (ipCollider == ipPROTON ))
450  {
451  double reduced_mass_collider_system = dense.AtomicWeight[nelem]*colliders.list[ipCollider].mass_amu/
452  (dense.AtomicWeight[nelem]+colliders.list[ipCollider].mass_amu)*ATOMIC_MASS_UNIT;
453  double ratef = powpq(ELECTRON_MASS/reduced_mass_collider_system,3,2) * COLL_CONST/phycon.sqrte;
454 
455  rate_test = cs*ratef*oHi;
456  fprintf(ioQQQ, "tauLo = %g\n",tauLo);
457  fprintf(ioQQQ,"Rates for H %ld %ld %ld %ld %ld %ld %ld %ld %g: %g %g \n",
458  nLo,lLo,lHi,sLo,sHi,jLo,jHi,nelem,phycon.te,rate_test, dense.eden);
459  /*fprintf(ioQQQ," E_lo,E_Hi,DE: %g , %g, %g \n", IP_Ryd_Hi,IP_Ryd_Lo,
460  IP_Ryd_Hi-IP_Ryd_Lo);*/
461  }
462  }
463 
464  return cs;
465 }
466 
467 STATIC double iso_get_collision_strength_collapsed_to_collapsed_fast( long ipISO, long nelem, long ipCollider,
468  long nHi, double IP_Ryd_Hi,
469  long nLo, double IP_Ryd_Lo,
470  double tauLo, double EnerWN, double EnerErg )
471 {
472  DEBUG_ENTRY( "iso_get_collision_strength_collapsed_to_collapsed_fast()" );
473 
474  t_iso_sp* sp = &iso_sp[ipISO][nelem];
475  ASSERT( nLo > sp->n_HighestResolved_max );
476  ASSERT( nLo > 2 );
477 
478  double cs = 0.;
479  double Aul_total = 0.;
480  long nAul = 0;
481  for( long lLo = 0; lLo < nLo; ++lLo )
482  {
483  long sLo=2;
484  if( ipISO==ipHE_LIKE )
485  sLo=1;
486  long sHi = sLo;
487  // NB - This is the index we would have if all n were resolved. 2nd dim of CachedAs is indexed by this.
488  long ipLoRes = sp->IndexIfAllResolved[nLo][lLo][sLo];
489  long gHi;
490 
491  if( lLo!= 0 )
492  {
493  gHi = (2*(lLo-1)+1)*sHi;
494  Aul_total += gHi * sp->CachedAs[ nHi-sp->n_HighestResolved_max-1 ][ ipLoRes ][1];
495  ++nAul;
496  }
497 
498  gHi = (2*(lLo+1)+1)*sHi;
499  Aul_total += gHi * sp->CachedAs[ nHi-sp->n_HighestResolved_max-1 ][ ipLoRes ][0];
500  ++nAul;
501 
502  // Account for triplet counterparts.
503  if( ipISO==ipHE_LIKE )
504  {
505  sLo=3;
506  sHi=sLo;
507  ipLoRes = sp->IndexIfAllResolved[nLo][lLo][sLo];
508  if( lLo!= 0 )
509  {
510  gHi = (2*(lLo-1)+1)*sHi;
511  Aul_total += gHi * sp->CachedAs[ nHi-sp->n_HighestResolved_max-1 ][ ipLoRes ][1];
512  ++nAul;
513  }
514 
515  gHi = (2*(lLo+1)+1)*sHi;
516  Aul_total += gHi * sp->CachedAs[ nHi-sp->n_HighestResolved_max-1 ][ ipLoRes ][0];
517  ++nAul;
518  }
519  }
520 
521  long g_nLo, factor;
522  if( ipISO==ipH_LIKE )
523  {
524  g_nLo = (2 * POW2(nLo));
525  factor = (nHi - 1) * g_nLo;
526  }
527  else if( ipISO==ipHE_LIKE )
528  {
529  g_nLo = (4 * POW2(nLo));
530  factor = (2 * nHi - 2) * g_nLo;
531  }
532  else
533  TotalInsanity();
534  ASSERT( factor > 0 );
535 
536  if (opac.lgCaseB_HummerStorey || (ipISO == ipH_LIKE && nelem != ipHYDROGEN ) )
537  cs = iso_get_collision_strength_resolved( ipISO, nelem, ipCollider,
538  nHi, -1, -1, -1, 1, IP_Ryd_Hi,
539  nLo, -1, -1, -1, 1, IP_Ryd_Lo,
540  Aul_total/nAul, tauLo, EnerWN, EnerErg );
541  else
542  {
543  cs = nAul * iso_get_collision_strength_resolved( ipISO, nelem, ipCollider,
544  nHi, -1, -1, -1, 1, IP_Ryd_Hi,
545  nLo, -1, -1, -1, 1, IP_Ryd_Lo,
546  Aul_total/nAul, tauLo, EnerWN, EnerErg );
547  cs += factor * iso_get_collision_strength_resolved( ipISO, nelem, ipCollider,
548  nHi, -1, -1, -1, 1, IP_Ryd_Hi,
549  nLo, -1, -1, -1, 1, IP_Ryd_Lo,
550  iso_ctrl.SmallA, tauLo, EnerWN, EnerErg );
551  }
552 
553  ASSERT( cs >= 0. );
554  return cs;
555 }
556 
557 #if 0
558 STATIC double iso_get_collision_strength_collapsed_to_collapsed( long ipISO, long nelem, long ipCollider,
559  long nHi, double IP_Ryd_Hi,
560  long nLo, double IP_Ryd_Lo,
561  double Aul, double tauLo, double EnerWN, double EnerErg )
562 {
563  DEBUG_ENTRY( "iso_get_collision_strength_collapsed_to_collapsed()" );
564 
565  t_iso_sp* sp = &iso_sp[ipISO][nelem];
566  ASSERT( nLo > sp->n_HighestResolved_max );
567 
568  vector<long> ang_mom_lo;
569  for( long i = 0; i < nLo; ++i )
570  ang_mom_lo.push_back( i );
571 
572  double cs = 0.;
573  for( vector<long>::iterator itLo = ang_mom_lo.begin(); itLo != ang_mom_lo.end(); ++itLo )
574  {
575  long lLo = *itLo;
576  long sLo = 2;
577  if( ipISO==ipHE_LIKE )
578  {
579  sLo = 1;
580  }
581  long jLo = -1;
582  long gLo = (2*lLo+1) * sLo;
583  double cs1 = iso_get_collision_strength_collapsed_to_resolved( ipISO, nelem, ipCollider,
584  nHi, IP_Ryd_Hi,
585  nLo, lLo, sLo, jLo, gLo, IP_Ryd_Lo,
586  Aul, tauLo, EnerWN, EnerErg );
587  cs += cs1;
588 
589  // Now do other spin for He-like
590  if( ipISO==ipHE_LIKE )
591  {
592  sLo = 3;
593  long gLo = (2*lLo+1) * sLo;
594  double cs1 = iso_get_collision_strength_collapsed_to_resolved( ipISO, nelem, ipCollider,
595  nHi, IP_Ryd_Hi,
596  nLo, lLo, sLo, jLo, gLo, IP_Ryd_Lo,
597  Aul, tauLo, EnerWN, EnerErg );
598  cs += cs1;
599  }
600  }
601 
602  return cs;
603 }
604 #endif
605 
606 STATIC double iso_get_collision_strength_collapsed_to_resolved( long ipISO, long nelem, long ipCollider,
607  long nHi, double IP_Ryd_Hi,
608  long nLo, long lLo, long sLo, long jLo, long gLo, double IP_Ryd_Lo,
609  double Aul, double tauLo, double EnerWN, double EnerErg )
610 {
611  DEBUG_ENTRY( "iso_get_collision_strength_collapsed_to_resolved()" );
612 
613  t_iso_sp* sp = &iso_sp[ipISO][nelem];
614  double cs;
615  long lHi = -1;
616  ASSERT( sLo > 0 );
617  long sHi = sLo;
618  long jHi = -1;
619  long gHi = 0;
620 
621  ASSERT( nHi >= 4 );
622  ASSERT( lLo >= 0 && lLo < nLo );
623 
624  // NB - This is the index we would have if all n were resolved. 2nd dim of CachedAs is indexed by this.
625  long ipLoRes = sp->IndexIfAllResolved[nLo][lLo][sLo];
626  if( (nLo==2) && (lLo==1) && (sLo==3) )
627  {
628  ASSERT( (jLo>=0) && (jLo<=2) );
629  ipLoRes -= (2 - jLo);
630  }
631 
632  double cs_less_1 = 0.;
633  if( lLo!= 0 )
634  {
635  lHi = lLo-1;
636  gHi = (2*lHi+1) * sHi;
637  Aul = sp->CachedAs[ nHi-sp->n_HighestResolved_max-1 ][ ipLoRes ][1];
638  cs_less_1 = iso_get_collision_strength_resolved( ipISO, nelem, ipCollider,
639  nHi, lHi, sHi, jHi, gHi, IP_Ryd_Hi,
640  nLo, lLo, sLo, jLo, gLo, IP_Ryd_Lo,
641  Aul, tauLo, EnerWN, EnerErg );
642  }
643 
644  lHi = lLo+1;
645  gHi = (2*lHi+1) * sHi;
646  Aul = sp->CachedAs[ nHi-sp->n_HighestResolved_max-1 ][ ipLoRes ][0];
647  double cs_plus_1 = iso_get_collision_strength_resolved( ipISO, nelem, ipCollider,
648  nHi, lHi, sHi, jHi, gHi, IP_Ryd_Hi,
649  nLo, lLo, sLo, jLo, gLo, IP_Ryd_Lo,
650  Aul, tauLo, EnerWN, EnerErg );
651 
652  if( lLo-2 >= 0 )
653  lHi = lLo-2;
654  else if( lLo+2 < nHi )
655  lHi = lLo+2;
656  else
657  // This can't be true because nHi is asserted >= 4 above.
658  TotalInsanity();
659 
660  gHi = 1;
661  Aul = iso_ctrl.SmallA;
662  double cs_other = iso_get_collision_strength_resolved( ipISO, nelem, ipCollider,
663  nHi, lHi, sHi, jHi, gHi, IP_Ryd_Hi,
664  nLo, lLo, sLo, jLo, gLo, IP_Ryd_Lo,
665  Aul, tauLo, EnerWN, EnerErg );
666 
667  long factor = nHi-2;
668  if( lLo==0 )
669  factor = nHi-1;
670  // There are twice as many transitions from He-like than from H-like.
671  if( ipISO==ipHE_LIKE )
672  factor *= 2;
673  cs = cs_less_1 + cs_plus_1 + factor * cs_other;
674 
675  return cs;
676 }
677 
678 STATIC double iso_get_collision_strength_resolved( long ipISO, long nelem, long ipCollider,
679  long nHi, long lHi, long sHi, long jHi, long gHi, double IP_Ryd_Hi,
680  long nLo, long lLo, long sLo, long jLo, long gLo, double IP_Ryd_Lo,
681  double Aul, double tauLo, double EnerWN, double EnerErg )
682 {
683  DEBUG_ENTRY( "iso_get_collision_strength_resolved()" );
684 
685  double cs;
686 
687  if( ipISO == ipH_LIKE )
688  {
689  cs = GetHlikeCollisionStrength( nelem, ipCollider,
690  nHi, lHi, sHi, gHi, IP_Ryd_Hi,
691  nLo, lLo, sLo, gLo, IP_Ryd_Lo,
692  Aul, tauLo, EnerErg );
693  }
694  else if( ipISO == ipHE_LIKE )
695  {
696  cs = GetHelikeCollisionStrength( nelem, ipCollider,
697  nHi, lHi, sHi, jHi, gHi, IP_Ryd_Hi,
698  nLo, lLo, sLo, jLo, gLo, IP_Ryd_Lo,
699  Aul, tauLo, EnerWN, EnerErg );
700  }
701  else
702  TotalInsanity();
703 
704  return cs;
705 }
706 
realnum x12tot
Definition: secondaries.h:65
STATIC double iso_get_collision_strength(long ipISO, long nelem, long ipCollider, long ipHi, long ipLo)
realnum GetHlikeCollisionStrength(long nelem, long ipCollider, long nHi, long lHi, long sHi, long gHi, double IP_Ryd_Hi, long nLo, long lLo, long sLo, long gLo, double IP_Ryd_Lo, double Aul, double tauLo, double EnerErg)
realnum EnergyErg() const
Definition: transition.h:90
qList st
Definition: iso.h:482
T * ptr0()
Definition: vectorize.h:221
const int ipHE_LIKE
Definition: iso.h:65
NORETURN void TotalInsanity(void)
Definition: service.cpp:1067
double EdenHCorr
Definition: dense.h:227
bool lgHeBug
Definition: trace.h:79
const double BIGDOUBLE
Definition: cpu.h:249
t_opac opac
Definition: opacity.cpp:5
t_Heavy Heavy
Definition: heavy.cpp:5
multi_arr< realnum, 3 > CachedAs
Definition: iso.h:596
bool lgColl_ionize[NISO]
Definition: iso.h:365
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
const int NISO
Definition: cddefines.h:310
bool lgIsoTraceFull[NISO]
Definition: trace.h:85
double hydro_vs_ioniz(double ionization_energy_Ryd, double Te)
#define PrintEfmt(F, V)
Definition: cddefines.h:1364
const double SMALLDOUBLE
Definition: cpu.h:250
long int nCollapsed_max
Definition: iso.h:518
t_conv conv
Definition: conv.cpp:5
STATIC double iso_get_collision_strength_collapsed_to_resolved(long ipISO, long nelem, long ipCollider, long nHi, 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)
t_phycon phycon
Definition: phycon.cpp:6
void iso_suprathermal(long ipISO, long nelem)
double coll_ion_wrapper(long int z, long int n, double t)
long int ipIsoTrace[NISO]
Definition: trace.h:88
bool lgColl_l_mixing[NISO]
Definition: iso.h:359
FILE * ioQQQ
Definition: cddefines.cpp:7
void iso_put_error(long ipISO, long nelem, long ipHi, long ipLo, long whichData, realnum errorOpt, realnum errorPess)
vector< freeBound > fb
Definition: iso.h:481
#define MIN2(a, b)
Definition: cddefines.h:807
double anu(size_t i) const
Definition: mesh.h:111
multi_arr< long, 3 > IndexIfAllResolved
Definition: iso.h:492
double * rate_coef_ul_set() const
Definition: collision.h:202
t_dense dense
Definition: global.cpp:15
const double MAX_DENSITY
Definition: cddefines.h:318
static t_ADfA & Inst()
Definition: cddefines.h:209
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
void vexp(const double x[], double y[], long nlo, long nhi)
bool lgPopLTE_OK
Definition: iso.h:554
realnum SmallA
Definition: iso.h:391
t_trace trace
Definition: trace.cpp:5
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:858
bool lgColl_excite[NISO]
Definition: iso.h:362
t_ionbal ionbal
Definition: ionbal.cpp:8
ColliderList colliders(dense)
realnum & gf() const
Definition: emission.h:558
long int n_HighestResolved_max
Definition: iso.h:536
realnum & EnergyWN() const
Definition: transition.h:477
double * Boltzmann()
Definition: quantumstate.h:159
#define POW2
Definition: cddefines.h:983
bool lgLevelsLowered
Definition: iso.h:505
#define STATIC
Definition: cddefines.h:118
bool lgTrace
Definition: trace.h:12
EmissionList::reference Emis() const
Definition: transition.h:447
#define N_(A_)
Definition: iso.h:22
const char * chISO[NISO]
Definition: iso.h:348
t_rfield rfield
Definition: rfield.cpp:9
#define IPCOLLIS
Definition: iso.h:89
long & ipCont() const
Definition: transition.h:489
float realnum
Definition: cddefines.h:124
long max(int a, long b)
Definition: cddefines.h:821
STATIC double iso_get_collision_strength_collapsed_to_collapsed_fast(long ipISO, long nelem, long ipCollider, long nHi, double IP_Ryd_Hi, long nLo, double IP_Ryd_Lo, double tauLo, double EnerWN, double EnerErg)
void iso_collisional_ionization(long ipISO, long nelem)
Definition: iso_collide.cpp:44
long int nTotalIoniz
Definition: conv.h:159
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:473
Definition: iso.h:470
int nLyaLevel[NISO]
Definition: iso.h:397
const int ipH2p
Definition: iso.h:31
double *** CollIonRate_Ground
Definition: ionbal.h:118
realnum AtomicWeight[LIMELM]
Definition: dense.h:80
#define ASSERT(exp)
Definition: cddefines.h:617
double HydroEinstA(long int n1, long int n2)
Definition: hydroeinsta.cpp:13
double GetGF(double trans_prob, double enercm, double gup)
vector< t_collider > list
Definition: collision.h:41
bool lgHBug
Definition: trace.h:82
const int ipH_LIKE
Definition: iso.h:64
const int LIMELM
Definition: cddefines.h:307
realnum & col_str() const
Definition: collision.h:191
CollisionProxy Coll() const
Definition: transition.h:463
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
double powpq(double x, int p, int q)
Definition: service.cpp:726
double te_ryd
Definition: phycon.h:27
realnum & rate_lu_nontherm_set() const
Definition: collision.h:211
const double * rate_coef_ul() const
Definition: collision.h:206
double eden
Definition: dense.h:201
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
realnum stat_ion[NISO]
Definition: iso.h:382
long int numLevels_max
Definition: iso.h:524
double sqrte
Definition: phycon.h:58
bool lgTopoff[NISO]
Definition: iso.h:435
t_secondaries secondaries
Definition: secondaries.cpp:5
double * ConBoltz()
Definition: quantumstate.h:151
double te
Definition: phycon.h:21
STATIC double iso_get_collision_strength_resolved(long ipISO, long nelem, long ipCollider, 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 int ipHYDROGEN
Definition: cddefines.h:348
double Hion_coll_ioniz_ratecoef(long int ipISO, long int nelem, long int n, double ionization_energy_Ryd, double Te)
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
realnum & Aul() const
Definition: emission.h:668
void iso_collide(long ipISO, long nelem)
long int numLevels_local
Definition: iso.h:529
bool lgCaseB_HummerStorey
Definition: opacity.h:177
double ColUL(const ColliderList &colls) const
Definition: collision.h:106
double te32
Definition: phycon.h:58
long int ipHeavy[LIMELM][LIMELM]
Definition: heavy.h:11
double density_low_limit
Definition: dense.h:208