cloudy  trunk
 All Data Structures 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-2008 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.h"
5 #include "conv.h"
6 #include "dense.h"
7 #include "heavy.h"
8 #include "helike_cs.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 "physconst.h"
16 #include "rfield.h"
17 #include "secondaries.h"
18 #include "trace.h"
19 #include "taulines.h"
20 
21 /* These are masses relative to the proton mass of the electron, proton, he+, and alpha particle. */
22 static double ColliderMass[4] = {ELECTRON_MASS/PROTON_MASS, 1.0, 4.0, 4.0};
23 
24 void iso_collisional_ionization( long ipISO, long nelem )
25 {
26  ASSERT( ipISO <= NISO );
27 
28  DEBUG_ENTRY( "iso_collisional_ionization()" );
29 
30  /* collisional ionization from ground */
31  iso.ColIoniz[ipISO][nelem][0] = iso.lgColl_ionize[ipISO] *
32  t_ADfA::Inst().coll_ion( nelem+1, 1+ipISO, phycon.te );
33 
34  for( long ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
35  {
36  if( nelem == ipISO )
37  {
38  /* use routine from Vriens and Smeets (1981). */
39  /* >>refer iso neutral col.ion. Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940 */
40  iso.ColIoniz[ipISO][nelem][ipHi] = hydro_vs_ioniz( ipISO, nelem, ipHi );
41  }
42  else
43  {
44  /* ions */
45  /* use hydrogenic ionization rates for ions
46  * >>refer iso ions col.ion. Allen 1973, Astro. Quan. for low Te.
47  * >>refer iso ions col.ion. Sampson and Zhang 1988, ApJ, 335, 516 for High Te.
48  * */
49  iso.ColIoniz[ipISO][nelem][ipHi] = Hion_coll_ioniz_ratecoef( ipISO, nelem, ipHi );
50  }
51 
52  /* iso.lgColl_ionize is option to turn off collisions, "atom XX-like collis off" comnd */
53  /* always leave highest level coupled to continuum */
54  if( ipHi < iso.numLevels_max[ipISO][nelem] - 1 )
55  iso.ColIoniz[ipISO][nelem][ipHi] *= iso.lgColl_ionize[ipISO];
56  }
57 
58  /* Here we arbitrarily scale the highest level ionization to account for the fact
59  * that, if the atom is not full size, this level should be interacting with higher
60  * levels and not just the continuum. We did add on collisional excitation terms instead
61  * but this caused a problem at low temperatures because the collisional ionization was
62  * a sum of terms with different Boltzmann factors, while PopLTE had just one Boltzmann
63  * factor. The result was a collisional recombination that had residual exponentials of
64  * the form exp(x/kT), which blew up at small T. */
65  if( !iso.lgLevelsLowered[ipISO][nelem] )
66  {
67  iso.ColIoniz[ipISO][nelem][iso.numLevels_max[ipISO][nelem]-1] *= 10.;/*1.E5;*/
68  iso.ColIoniz[ipISO][nelem][iso.numLevels_max[ipISO][nelem]-1] *= 10.;/*1.E5;*/
69  }
70 
71  return;
72 }
73 
74 void iso_suprathermal( long ipISO, long nelem )
75 {
76  DEBUG_ENTRY( "iso_suprathermal()" );
77 
78  /* check that we were called with valid parameters */
79  ASSERT( ipISO < NISO );
80  ASSERT( nelem >= ipISO );
81  ASSERT( nelem < LIMELM );
82 
83  /***********************************************************************
84  * *
85  * get secondary excitation by suprathermal electrons *
86  * *
87  ***********************************************************************/
88 
89  for( long i=1; i < iso.numLevels_max[ipISO][nelem]; i++ )
90  {
91  if( Transitions[ipISO][nelem][i][0].ipCont > 0 )
92  {
93  /* get secondaries for all iso lines by scaling LyA
94  * excitation by ratio of cross section (oscillator strength/energy)
95  * Born approximation or plane-wave approximation */
96  secondaries.Hx12[ipISO][nelem][i] = secondaries.x12tot *
97  (Transitions[ipISO][nelem][i][0].Emis->gf/
98  Transitions[ipISO][nelem][i][0].EnergyWN) /
99  (Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][0].Emis->gf/
100  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][0].EnergyWN);
101  }
102  else
103  secondaries.Hx12[ipISO][nelem][i] = 0.;
104  }
105 
106  return;
107 }
108 
109 /*============================*/
110 /* evaluate collisional rates */
111 void iso_collide( long ipISO, long nelem )
112 {
113  long ipHi, ipLo;
114  double factor, ConvLTEPOP, edenCorr;
115 
116  /* this array stores the last temperature at which collision data were evaluated for
117  * each species of the isoelectronic sequence. */
118  static double TeUsed[NISO][LIMELM]={
119  {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},
120  {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} };
121 
122  /* collision strengths are assumed to be roughly constant for small changes in temperature
123  * and are not recalculated as often as other data. This array stores the last temperature
124  * at which collision strengths were evaluated for each species of the isoelectronic sequence. */
125  static double TeUsedForCS[NISO][LIMELM]={
126  {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},
127  {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} };
128 
129  DEBUG_ENTRY( "iso_collide()" );
130 
131  /* check that we were called with valid parameters */
132  ASSERT( ipISO < NISO );
133  ASSERT( nelem >= ipISO );
134  ASSERT( nelem < LIMELM );
135 
136  /* skip most of this routine if temperature has not changed,
137  * the check on conv.nTotalIoniz is to make sure that we redo this
138  * on the very first call in a grid calc - it is 0 on the first call */
139  if( fp_equal( TeUsed[ipISO][nelem], phycon.te ) && conv.nTotalIoniz )
140  {
141  ASSERT( Transitions[ipISO][nelem][ iso.nLyaLevel[ipISO] ][0].Coll.ColUL >= 0. );
142 
143  if( trace.lgTrace )
144  {
145  fprintf( ioQQQ,
146  " iso_collide called %s nelem %li - no reeval Boltz fac, LTE dens\n",
147  iso.chISO[ipISO], nelem );
148  }
149  }
150  else
151  {
152  TeUsed[ipISO][nelem] = phycon.te;
153 
154  if( trace.lgTrace )
155  {
156  fprintf( ioQQQ,
157  " iso_collide called %s nelem %li - will reeval Boltz fac, LTE dens\n",
158  iso.chISO[ipISO], nelem );
159  }
160 
161  /**********************************************************
162  * *
163  * Boltzmann factors for all levels, and *
164  * collisional ionization and excitation *
165  * *
166  **********************************************************/
167 
168  for( ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
169  {
170  for( ipLo=0; ipLo<ipHi; ipLo++ )
171  {
172  iso.Boltzmann[ipISO][nelem][ipHi][ipLo] =
173  sexp( Transitions[ipISO][nelem][ipHi][ipLo].EnergyK / phycon.te );
174  }
175  }
176 
177  /* HION_LTE_POP is planck^2 / (2 pi m_e k ), raised to 3/2 next */
178  factor = HION_LTE_POP*dense.AtomicWeight[nelem]/
180 
181  /* term in () is stat weight of electron * ion */
182  ConvLTEPOP = pow(factor,1.5)/(2.*iso.stat_ion[ipISO])/phycon.te32;
183 
184  iso.lgPopLTE_OK[ipISO][nelem] = true;
185 
186  /* fully define Boltzmann factors to continuum for model levels */
187  for( ipLo=0; ipLo<iso.numLevels_max[ipISO][nelem]; ipLo++ )
188  {
189  /* this Boltzmann factor is exp( +ioniz energy / Te ) */
190  StatesElem[ipISO][nelem][ipLo].ConBoltz =
191  dsexp(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo]/phycon.te_ryd);
192 
193  /***************************************
194  * *
195  * LTE abundances for all levels *
196  * *
197  ***************************************/
198 
199  if( StatesElem[ipISO][nelem][ipLo].ConBoltz > SMALLDOUBLE )
200  {
201  /* LTE population of given level. */
202  iso.PopLTE[ipISO][nelem][ipLo] =
203  StatesElem[ipISO][nelem][ipLo].g / StatesElem[ipISO][nelem][ipLo].ConBoltz * ConvLTEPOP;
204  ASSERT( iso.PopLTE[ipISO][nelem][ipLo] < BIGDOUBLE );
205  }
206  else
207  {
208  iso.PopLTE[ipISO][nelem][ipLo] = 0.;
209  }
210 
211  /* now check for any zeros - if present then matrix cannot be used */
212  if( iso.PopLTE[ipISO][nelem][ipLo] <= 0. )
213  {
214  iso.lgPopLTE_OK[ipISO][nelem] = false;
215  }
216  }
217 
218  iso_collisional_ionization( ipISO, nelem );
219 
220  /***********************************************************
221  * *
222  * collision strengths for all lines in iso sequence *
223  * *
224  ***********************************************************/
225 
226  /* Calculate initial value of collision strengths, OR
227  * Reevaluate collision strengths if the temperature has changed by more than 15%. */
228  if( (TeUsedForCS[ipISO][nelem] == 0.) ||
229  ( TeUsedForCS[ipISO][nelem]/phycon.te > 1.15 ) ||
230  ( TeUsedForCS[ipISO][nelem]/phycon.te < 0.85 ) )
231  {
232  /* Update temperature at which collision strengths are evaluated. */
233  TeUsedForCS[ipISO][nelem] = phycon.te;
234 
235  for( ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
236  {
237  for( ipLo=0; ipLo < ipHi; ipLo++ )
238  {
239  if( ipISO == ipH_LIKE )
240  {
241  /* collision strengths due to electron impact */
242  Transitions[ipISO][nelem][ipHi][ipLo].Coll.cs =
243  HydroCSInterp( nelem , ipHi , ipLo, ipELECTRON );
244  /* proton impact */
245  Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipPROTON] =
246  HydroCSInterp( nelem , ipHi , ipLo, ipPROTON );
247  /* He+ impact */
248  Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipHE_PLUS] =
249  HydroCSInterp( nelem , ipHi , ipLo, ipHE_PLUS );
250  /* and He++ impact */
251  Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipALPHA] =
252  HydroCSInterp( nelem , ipHi , ipLo, ipALPHA );
253  }
254  else if( ipISO == ipHE_LIKE )
255  {
256  /* collision strengths due to electron impact */
257  Transitions[ipISO][nelem][ipHi][ipLo].Coll.cs =
258  HeCSInterp( nelem , ipHi , ipLo, ipELECTRON );
259  /* proton impact */
260  Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipPROTON] =
261  HeCSInterp( nelem , ipHi , ipLo, ipPROTON );
262  /* and He+ impact */
263  Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipHE_PLUS] =
264  HeCSInterp( nelem , ipHi , ipLo, ipHE_PLUS );
265 
266  Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipALPHA] = 0.;
267  }
268  else
269  TotalInsanity();
270 
271  if( N_(ipHi) != N_(ipLo) )
272  {
273  /* kill all collisional excitation if this set with
274  * atom XX-like collisional excitation */
275  Transitions[ipISO][nelem][ipHi][ipLo].Coll.cs *= iso.lgColl_excite[ipISO];
276  Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipPROTON] *= iso.lgColl_excite[ipISO];
277  Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipHE_PLUS] *= iso.lgColl_excite[ipISO];
278  Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipALPHA] *= iso.lgColl_excite[ipISO];
279 
280  }
281  else
282  {
283  Transitions[ipISO][nelem][ipHi][ipLo].Coll.cs *= iso.lgColl_l_mixing[ipISO];
284  Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipPROTON] *= iso.lgColl_l_mixing[ipISO];
285  Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipHE_PLUS] *= iso.lgColl_l_mixing[ipISO];
286  Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipALPHA] *= iso.lgColl_excite[ipISO];
287  }
288 
289  /* check for sanity */
290  ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].Coll.cs >= 0. );
291  ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipPROTON] >= 0. );
292  ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipHE_PLUS] >= 0. );
293  ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipALPHA] >= 0. );
294  }
295  }
296  }
297 
298  /* the case b hummer and storey option,
299  * this kills collisional excitation and ionization from n=1 and n=2 */
300  if( opac.lgCaseB_HummerStorey && ipISO==ipH_LIKE )
301  {
302  for( ipLo=0; ipLo<iso.numLevels_max[ipISO][nelem]-1; ipLo++ )
303  {
304  if( N_(ipLo)>=3 )
305  break;
306 
307  iso.ColIoniz[ipISO][nelem][ipLo] = 0.;
308 
309  for( ipHi=ipLo+1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
310  {
311  Transitions[ipISO][nelem][ipHi][ipLo].Coll.cs = 0.;
312  Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipPROTON] = 0.;
313  Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipHE_PLUS] = 0.;
314  Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipALPHA] = 0.;
315  }
316  }
317  }
318 
319  /***********************************************************************
320  * *
321  * collisional deexcitation for all lines *
322  * *
323  ***********************************************************************/
324  for( ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
325  {
326  for( ipLo=0; ipLo<ipHi; ipLo++ )
327  {
328  double reduced_mass_proton = dense.AtomicWeight[nelem]*ColliderMass[ipPROTON]/
330 
331  double reduced_mass_heplus = dense.AtomicWeight[nelem]*ColliderMass[ipHE_PLUS]/
333 
334  /********************************************************
335  ********************************************************
336  * NB - the collision strengths for proton and helium *
337  * ion impact are multiplied by the ratio of collider *
338  * to electron densities to precorrect for getting *
339  * multiplied by the electron density when being put *
340  * into rate matrix later. They are also multiplied *
341  * by the pow(m_e/reduced mass)^1.5 to correct for the *
342  * fact that COLL_CONST contains the electron mass. *
343  * Here, reduced mass means the reduced mass of the *
344  * collider-target system. *
345  ********************************************************
346  ********************************************************/
347  Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL = (realnum)(
348  (
349  /* due to electron impact */
350  Transitions[ipISO][nelem][ipHi][ipLo].Coll.cs +
351  /* due to proton impact */
352  Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipPROTON]*
354  pow(ELECTRON_MASS/reduced_mass_proton, 1.5) +
355  /* due to he+ impact */
356  Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipHE_PLUS]*
358  pow(ELECTRON_MASS/reduced_mass_heplus, 1.5) +
359  /* due to he++ impact */
360  Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipALPHA]*
362  pow(ELECTRON_MASS/reduced_mass_heplus, 1.5)
363  ) / phycon.sqrte*COLL_CONST/(double)StatesElem[ipISO][nelem][ipHi].g );
364 
365  if( ipISO == ipH_LIKE )
366  {
367  if( N_(ipHi) > iso.n_HighestResolved_max[ipISO][nelem] &&
368  N_(ipLo) <= iso.n_HighestResolved_max[ipISO][nelem] )
369  {
370  Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL *=
371  (8.f/3.f)*(log((realnum)N_(ipHi))+2.f);
372  }
373  }
374  else if( ipISO == ipHE_LIKE )
375  {
376  fixit();
377  /* This is intended to be a trick to get the correct collisional excitation from
378  * collapsed levels to resolved levels. Is it needed or does the stat weight used
379  * above handle it automatically? If it is needed, is this correct? */
380  if( N_(ipHi) > iso.n_HighestResolved_max[ipISO][nelem] &&
381  N_(ipLo) <= iso.n_HighestResolved_max[ipISO][nelem] )
382  {
383  Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL *=
384  (2.f/3.f)*(log((realnum)N_(ipHi))+2.f);
385  }
386  }
387  }
388  }
389 
390  if( (trace.lgTrace && trace.lgIsoTraceFull[ipISO]) && (nelem == trace.ipIsoTrace[ipISO]) )
391  {
392  fprintf( ioQQQ, " iso_collide: %s Z=%li de-excitation rates coefficients\n", iso.chISO[ipISO], nelem + 1 );
393  for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
394  {
395  fprintf( ioQQQ, " %li\t", ipHi );
396  for( ipLo=0; ipLo < ipHi; ipLo++ )
397  {
398  fprintf( ioQQQ,PrintEfmt("%10.3e", Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL ));
399  }
400  fprintf( ioQQQ, "\n" );
401  }
402 
403  fprintf( ioQQQ, " iso_collide: %s Z=%li collisional ionization coefficients\n", iso.chISO[ipISO], nelem + 1 );
404  for( ipHi=0; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
405  {
406  fprintf( ioQQQ,PrintEfmt("%10.3e", iso.ColIoniz[ipISO][nelem][ipHi] ));
407  }
408  fprintf( ioQQQ, "\n" );
409 
410  fprintf( ioQQQ, " iso_collide: %s Z=%li continuum boltzmann factor\n", iso.chISO[ipISO], nelem + 1 );
411  for( ipHi=0; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
412  {
413  fprintf( ioQQQ,PrintEfmt("%10.3e", StatesElem[ipISO][nelem][ipHi].ConBoltz ));
414  }
415  fprintf( ioQQQ, "\n" );
416 
417  fprintf( ioQQQ, " iso_collide: %s Z=%li continuum boltzmann factor\n", iso.chISO[ipISO], nelem + 1 );
418  for( ipHi=0; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
419  {
420  fprintf( ioQQQ,PrintEfmt("%10.3e", iso.PopLTE[ipISO][nelem][ipHi] ));
421  }
422  fprintf( ioQQQ, "\n" );
423  }
424  }
425 
426  iso_suprathermal( ipISO, nelem );
427 
428  /* >>chng 05 aug 17, must use real electron density for collisional ionization of H
429  * since in Leiden v4 pdr with its artificial high temperature coll ion can be important
430  * H on H is homonuclear and scaling laws for other elements does not apply
431  * next two were multiplied by dense.EdenHCorr and changed to dense.eden for H,
432  * EdenHCorr for rest */
433  if( nelem==ipHYDROGEN )
434  {
435  /* special version for H0 onto H0 */
436  edenCorr = dense.EdenHontoHCorr;
437  }
438  else
439  {
440  edenCorr = dense.EdenHCorr;
441  }
442 
443  /* this must be reevaluated every time since eden can change when Te does not */
444  /* save into main array - collisional ionization by thermal electrons */
445  ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0] =
446  iso.ColIoniz[ipISO][nelem][0]*edenCorr;
447 
448  /* cooling due to collisional ionization, which only includes thermal electrons */
449  ionbal.CollIonRate_Ground[nelem][nelem-ipISO][1] =
450  ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0]*
451  rfield.anu[Heavy.ipHeavy[nelem][nelem-ipISO]-1]*EN1RYD;
452 
453  return;
454 }

Generated for cloudy by doxygen 1.8.1.2