cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*atmdat_DielSupres derive scale factors for suppression of Burgess dielectronic recombination */ 00004 #include "cddefines.h" 00005 #include "ionbal.h" 00006 #include "dense.h" 00007 #include "phycon.h" 00008 #include "punch.h" 00009 #include "atmdat.h" 00010 00011 /* >>chng 07 feb 28 by Mitchell Martin, added new dr suppression routine */ 00012 /* This function computes the standard electron density-dependent 00013 * suppression factor of the DR rate coefficient of the C3+ ion 00014 * at T = 10^5 K, based on Nigel Badnell's 1993 ApJ paper. 00015 * It is then scalable for other choices of ionic charge and temperature. 00016 */ 00017 //#define USENEW 00018 #ifdef USENEW 00019 STATIC double dr_suppress( 00020 /* This routine takes the following arguments: 00021 * atomic_number = nuclear charge */ 00022 long int atomic_number, 00023 /*ionic_charge = ionic charge*/ 00024 long int ionic_charge, 00025 /*eden = electron density */ 00026 double eden, 00027 /*T = temperature (K)*/ 00028 double T ) 00029 { 00030 00031 00032 /* fitting constants to compute nominal suppression factor function */ 00033 const double A = 12.4479; /* various fitting parameters follow */ 00034 const double mu = 0.46665; 00035 const double w = 4.96916; 00036 const double y_c0 = 0.5498; /* log10 of the electron density fitting parameter for C3+ */ 00037 00038 /* Nigel's 1993 ApJ paper computed the DR rate as a function of log n_e 00039 * at a temperature T = 100,000 K. 00040 */ 00041 const double T_0 = 1.e5; /* the standard temperature in Nigel's original C3+ fit */ 00042 const double q_0 = 3.; /* the ionic charge of C3+, addressed in Nigel's paper */ 00043 00044 /* a fitting constant to compute the suppression factor corrected for an 00045 * estimate of surviving DR based on the lowest dipole allowed core 00046 * excitation energy 00047 */ 00048 const double c = 3.0; /* smaller c means larger fraction will survive, and vice versa */ 00049 00050 double s, snew, y_c, E_c, E_c0, x, g_x; 00051 long int iso_sequence; 00052 00053 eden = log10(eden); 00054 iso_sequence = atomic_number - ionic_charge; /* the isoelectronic sequence number, iso_sequence=3 for Li-like, etc */ 00055 00056 y_c = y_c0 + log10( pow( (ionic_charge/q_0), 7. ) * sqrt( T/T_0 ) ); 00057 00058 /* here we compute the standard suppression factor function, s( n_e, T, ionic_charge ) */ 00059 if( eden >= y_c ) 00060 { 00061 s = (A/(PI*w)) * ( mu/( 1. + pow((eden-y_c)/w, 2.) ) + 00062 (1. - mu) * sqrt(PI*LN_TWO) * exp( -LN_TWO * 00063 pow((eden-y_c)/w, 2.) ) ); 00064 } 00065 else 00066 { 00067 s = (A/(PI*w)) * ( mu + (1.- mu) * sqrt(PI*LN_TWO) ); 00068 } 00069 00070 /* Now we're going to modify this standard suppression factor curve to 00071 * allow for the survival of some fraction of the total DR rate at 00072 * generally lower temperatures T, when appropriate. 00073 */ 00074 00075 /* Computational estimates of lowest dipole allowed core excitation 00076 * energies for various iso-electronic sequences of recombining ion; 00077 * these are fits to NIST statistical weighted energies 00078 */ 00079 if( iso_sequence == 3 ) /* Li-like ions */ 00080 { 00081 E_c = 2.08338 + 19.1356*(ionic_charge/10.) + 0.974 * 00082 pow( ionic_charge/10., 2. ) - 0.309032*pow( ionic_charge/10., 3. ) + 00083 0.419951*pow( ionic_charge/10., 4. ); 00084 } 00085 else if( iso_sequence == 4 ) /* Be-like ions */ 00086 { 00087 E_c = 5.56828 + 34.6774*(ionic_charge/10.) + 1.005 * 00088 pow( ionic_charge/10., 2. ) - 0.994177*pow( ionic_charge/10., 3. ) + 00089 0.726053*pow( ionic_charge/10., 4. ); 00090 } 00091 else if( iso_sequence == 7 ) /* N-like ions */ 00092 { 00093 E_c = 10.88361 + 39.7851*(ionic_charge/10.) + 0.423 * 00094 pow( ionic_charge/10., 2. ) - 0.310368*pow( ionic_charge/10., 3. ) + 00095 0.937186*pow( ionic_charge/10., 4. ); 00096 } 00097 else if( iso_sequence == 11 ) /* Na-like ions */ 00098 { 00099 E_c = 2.17262 + 22.5038*(ionic_charge/10.) - 1.227*pow( ionic_charge/10., 2. ) + 00100 0.801291*pow( ionic_charge/10., 3. ) + 00101 0.0434168*pow( ionic_charge/10., 4. ); 00102 } 00103 else if( iso_sequence == 1 || iso_sequence == 2 || iso_sequence == 10 ) 00104 { 00105 /* set to a very large number to force suppression factor to 1.0 00106 * for H, He, Ne-like ions */ 00107 E_c = 1.e10; 00108 } 00109 else 00110 { 00111 /* specifically B, C, O, or F-like ions (iso_sequence = 5, 6, 8, 9) */ 00112 E_c = 0.0; /* forces suppression factor to s for all T */ 00113 /* iso_sequence.B. ion sequences beyond Na-like, iso_sequence > 11, are currently not 00114 * treated */ 00115 } 00116 00117 /* the lowest dipole allowed energy for Li-like C3+, atomic_number = 6, iso_sequence = atomic_number-ionic_charge = 3 */ 00118 E_c0 = 2.08338 + 19.1356*(q_0/10.) + 0.974 * 00119 pow( (q_0/10.), 2. ) - 0.309032 * 00120 pow( (q_0/10.), 3. ) + 0.419951 * 00121 pow( (q_0/10.), 4. ); 00122 00123 /* and important factor that determines what survives */ 00124 x = ( (E_c*EVDEGK)/(c*T) - (E_c0*EVDEGK)/(c*T_0) ); 00125 00126 if( x > 1 ) 00127 { 00128 g_x = x; 00129 } 00130 else if( x >= 0 && x <= 1 ) 00131 { 00132 g_x = x*x; 00133 } 00134 else 00135 { 00136 g_x = 0.0; 00137 } 00138 00139 /* converting the standard curve to the revised one allowing for 00140 * survival at lower energies 00141 */ 00142 snew = 1. + (s-1.)*exp(-g_x); 00143 00144 return snew; 00145 ASSERT( snew >=0. && snew <= 1. ); 00146 } 00147 #endif 00148 00149 void atmdat_DielSupres(void) 00150 { 00151 long int i; 00152 00153 DEBUG_ENTRY( "atmdat_DielSupres()" ); 00154 00155 /* dielectronic burgess recombination suppression, default is true */ 00156 if( ionbal.lgSupDie[0] ) 00157 { 00158 for( i=0; i < LIMELM; i++ ) 00159 { 00160 # ifdef USENEW 00161 ionbal.DielSupprs[0][i] = (realnum)dr_suppress( i+1, 3, dense.eden, phycon.te ); 00162 # else 00163 /* effective density for scaling from Davidson's plot 00164 * first do temperature scaling, term in () is SQRT(te/15,000) */ 00165 double effden = dense.eden/(phycon.sqrte/122.47); 00166 00167 /* this is rough charge dependence, z^7 from Davidson */ 00168 effden /= powi((realnum)(i+1)/3.,7); 00169 00170 ionbal.DielSupprs[0][i] = (realnum)(1.-0.092*log10(effden)); 00171 ionbal.DielSupprs[0][i] = (realnum)MIN2(1.,ionbal.DielSupprs[0][i]); 00172 ionbal.DielSupprs[0][i] = (realnum)MAX2(0.08,ionbal.DielSupprs[0][i]); 00173 # endif 00174 } 00175 } 00176 00177 else 00178 { 00179 for( i=0; i < LIMELM; i++ ) 00180 { 00181 ionbal.DielSupprs[0][i] = 1.; 00182 } 00183 } 00184 00185 /* nussbaumer and storey recombination 00186 * default is this to be false */ 00187 if( ionbal.lgSupDie[1] ) 00188 { 00189 for( i=0; i < LIMELM; i++ ) 00190 { 00191 /* assume same factors as above */ 00192 ionbal.DielSupprs[1][i] = ionbal.DielSupprs[0][i]; 00193 } 00194 } 00195 else 00196 { 00197 for( i=0; i < LIMELM; i++ ) 00198 { 00199 ionbal.DielSupprs[1][i] = 1.; 00200 } 00201 } 00202 00203 /* option to punch recombination coefficients, set with *punch recombination 00204 * coefficients* command*/ 00205 if( punch.lgioRecom ) 00206 { 00207 fprintf( punch.ioRecom, " atmdat_DielSupres finds following dielectronic" 00208 " recom suppression factors.\n" ); 00209 fprintf( punch.ioRecom, " Z fac \n" ); 00210 for( i=0; i < LIMELM; i++ ) 00211 { 00212 fprintf( punch.ioRecom, "%3ld %10.3e %10.3e\n", i+1, 00213 ionbal.DielSupprs[0][i], ionbal.DielSupprs[1][i] ); 00214 } 00215 fprintf( punch.ioRecom, "\n"); 00216 } 00217 return; 00218 }