cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
atmdat_dielsupres.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 /*atmdat_DielSupres derive scale factors for suppression of Burgess dielectronic recombination */
4 #include "cddefines.h"
5 #include "ionbal.h"
6 #include "dense.h"
7 #include "phycon.h"
8 #include "punch.h"
9 #include "atmdat.h"
10 
11 /* >>chng 07 feb 28 by Mitchell Martin, added new dr suppression routine */
12 /* This function computes the standard electron density-dependent
13  * suppression factor of the DR rate coefficient of the C3+ ion
14  * at T = 10^5 K, based on Nigel Badnell's 1993 ApJ paper.
15  * It is then scalable for other choices of ionic charge and temperature.
16  */
17 //#define USENEW
18 #ifdef USENEW
19 STATIC double dr_suppress(
20  /* This routine takes the following arguments:
21  * atomic_number = nuclear charge */
22  long int atomic_number,
23  /*ionic_charge = ionic charge*/
24  long int ionic_charge,
25  /*eden = electron density */
26  double eden,
27  /*T = temperature (K)*/
28  double T )
29 {
30 
31 
32  /* fitting constants to compute nominal suppression factor function */
33  const double A = 12.4479; /* various fitting parameters follow */
34  const double mu = 0.46665;
35  const double w = 4.96916;
36  const double y_c0 = 0.5498; /* log10 of the electron density fitting parameter for C3+ */
37 
38  /* Nigel's 1993 ApJ paper computed the DR rate as a function of log n_e
39  * at a temperature T = 100,000 K.
40  */
41  const double T_0 = 1.e5; /* the standard temperature in Nigel's original C3+ fit */
42  const double q_0 = 3.; /* the ionic charge of C3+, addressed in Nigel's paper */
43 
44  /* a fitting constant to compute the suppression factor corrected for an
45  * estimate of surviving DR based on the lowest dipole allowed core
46  * excitation energy
47  */
48  const double c = 3.0; /* smaller c means larger fraction will survive, and vice versa */
49 
50  double s, snew, y_c, E_c, E_c0, x, g_x;
51  long int iso_sequence;
52 
53  eden = log10(eden);
54  iso_sequence = atomic_number - ionic_charge; /* the isoelectronic sequence number, iso_sequence=3 for Li-like, etc */
55 
56  y_c = y_c0 + log10( pow( (ionic_charge/q_0), 7. ) * sqrt( T/T_0 ) );
57 
58  /* here we compute the standard suppression factor function, s( n_e, T, ionic_charge ) */
59  if( eden >= y_c )
60  {
61  s = (A/(PI*w)) * ( mu/( 1. + pow((eden-y_c)/w, 2.) ) +
62  (1. - mu) * sqrt(PI*LN_TWO) * exp( -LN_TWO *
63  pow((eden-y_c)/w, 2.) ) );
64  }
65  else
66  {
67  s = (A/(PI*w)) * ( mu + (1.- mu) * sqrt(PI*LN_TWO) );
68  }
69 
70  /* Now we're going to modify this standard suppression factor curve to
71  * allow for the survival of some fraction of the total DR rate at
72  * generally lower temperatures T, when appropriate.
73  */
74 
75  /* Computational estimates of lowest dipole allowed core excitation
76  * energies for various iso-electronic sequences of recombining ion;
77  * these are fits to NIST statistical weighted energies
78  */
79  if( iso_sequence == 3 ) /* Li-like ions */
80  {
81  E_c = 2.08338 + 19.1356*(ionic_charge/10.) + 0.974 *
82  pow( ionic_charge/10., 2. ) - 0.309032*pow( ionic_charge/10., 3. ) +
83  0.419951*pow( ionic_charge/10., 4. );
84  }
85  else if( iso_sequence == 4 ) /* Be-like ions */
86  {
87  E_c = 5.56828 + 34.6774*(ionic_charge/10.) + 1.005 *
88  pow( ionic_charge/10., 2. ) - 0.994177*pow( ionic_charge/10., 3. ) +
89  0.726053*pow( ionic_charge/10., 4. );
90  }
91  else if( iso_sequence == 7 ) /* N-like ions */
92  {
93  E_c = 10.88361 + 39.7851*(ionic_charge/10.) + 0.423 *
94  pow( ionic_charge/10., 2. ) - 0.310368*pow( ionic_charge/10., 3. ) +
95  0.937186*pow( ionic_charge/10., 4. );
96  }
97  else if( iso_sequence == 11 ) /* Na-like ions */
98  {
99  E_c = 2.17262 + 22.5038*(ionic_charge/10.) - 1.227*pow( ionic_charge/10., 2. ) +
100  0.801291*pow( ionic_charge/10., 3. ) +
101  0.0434168*pow( ionic_charge/10., 4. );
102  }
103  else if( iso_sequence == 1 || iso_sequence == 2 || iso_sequence == 10 )
104  {
105  /* set to a very large number to force suppression factor to 1.0
106  * for H, He, Ne-like ions */
107  E_c = 1.e10;
108  }
109  else
110  {
111  /* specifically B, C, O, or F-like ions (iso_sequence = 5, 6, 8, 9) */
112  E_c = 0.0; /* forces suppression factor to s for all T */
113  /* iso_sequence.B. ion sequences beyond Na-like, iso_sequence > 11, are currently not
114  * treated */
115  }
116 
117  /* the lowest dipole allowed energy for Li-like C3+, atomic_number = 6, iso_sequence = atomic_number-ionic_charge = 3 */
118  E_c0 = 2.08338 + 19.1356*(q_0/10.) + 0.974 *
119  pow( (q_0/10.), 2. ) - 0.309032 *
120  pow( (q_0/10.), 3. ) + 0.419951 *
121  pow( (q_0/10.), 4. );
122 
123  /* and important factor that determines what survives */
124  x = ( (E_c*EVDEGK)/(c*T) - (E_c0*EVDEGK)/(c*T_0) );
125 
126  if( x > 1 )
127  {
128  g_x = x;
129  }
130  else if( x >= 0 && x <= 1 )
131  {
132  g_x = x*x;
133  }
134  else
135  {
136  g_x = 0.0;
137  }
138 
139  /* converting the standard curve to the revised one allowing for
140  * survival at lower energies
141  */
142  snew = 1. + (s-1.)*exp(-g_x);
143 
144  return snew;
145  ASSERT( snew >=0. && snew <= 1. );
146 }
147 #endif
148 
150 {
151  long int i;
152 
153  DEBUG_ENTRY( "atmdat_DielSupres()" );
154 
155  /* dielectronic burgess recombination suppression, default is true */
156  if( ionbal.lgSupDie[0] )
157  {
158  for( i=0; i < LIMELM; i++ )
159  {
160 # ifdef USENEW
161  ionbal.DielSupprs[0][i] = (realnum)dr_suppress( i+1, 3, dense.eden, phycon.te );
162 # else
163  /* effective density for scaling from Davidson's plot
164  * first do temperature scaling, term in () is SQRT(te/15,000) */
165  double effden = dense.eden/(phycon.sqrte/122.47);
166 
167  /* this is rough charge dependence, z^7 from Davidson */
168  effden /= powi((realnum)(i+1)/3.,7);
169 
170  ionbal.DielSupprs[0][i] = (realnum)(1.-0.092*log10(effden));
171  ionbal.DielSupprs[0][i] = (realnum)MIN2(1.,ionbal.DielSupprs[0][i]);
172  ionbal.DielSupprs[0][i] = (realnum)MAX2(0.08,ionbal.DielSupprs[0][i]);
173 # endif
174  }
175  }
176 
177  else
178  {
179  for( i=0; i < LIMELM; i++ )
180  {
181  ionbal.DielSupprs[0][i] = 1.;
182  }
183  }
184 
185  /* nussbaumer and storey recombination
186  * default is this to be false */
187  if( ionbal.lgSupDie[1] )
188  {
189  for( i=0; i < LIMELM; i++ )
190  {
191  /* assume same factors as above */
192  ionbal.DielSupprs[1][i] = ionbal.DielSupprs[0][i];
193  }
194  }
195  else
196  {
197  for( i=0; i < LIMELM; i++ )
198  {
199  ionbal.DielSupprs[1][i] = 1.;
200  }
201  }
202 
203  /* option to punch recombination coefficients, set with *punch recombination
204  * coefficients* command*/
205  if( punch.lgioRecom )
206  {
207  fprintf( punch.ioRecom, " atmdat_DielSupres finds following dielectronic"
208  " recom suppression factors.\n" );
209  fprintf( punch.ioRecom, " Z fac \n" );
210  for( i=0; i < LIMELM; i++ )
211  {
212  fprintf( punch.ioRecom, "%3ld %10.3e %10.3e\n", i+1,
213  ionbal.DielSupprs[0][i], ionbal.DielSupprs[1][i] );
214  }
215  fprintf( punch.ioRecom, "\n");
216  }
217  return;
218 }

Generated for cloudy by doxygen 1.8.1.2