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 /*IonIron ionization balance for iron */ 00004 #include "cddefines.h" 00005 #include "dense.h" 00006 #include "fe.h" 00007 #include "gammas.h" 00008 #include "opacity.h" 00009 #include "trace.h" 00010 #include "grainvar.h" 00011 #include "ionbal.h" 00012 #include "mole.h" 00013 00014 /*IonIron ionization balance for iron */ 00015 void IonIron(void) 00016 { 00017 const int NDIM = ipIRON+1; 00018 00019 static const double dicoef[2][NDIM] = { 00020 {1.58e-3,8.38e-3,1.54e-2,3.75e-2,0.117,0.254,0.291,0.150,0.140,0.100,0.200,0.240, 00021 0.260,0.190,0.120,0.350,0.066,0.10,0.13,0.23,0.14,0.11,0.041,0.747,0.519,0.}, 00022 {.456,.323,.310,.411,.359,.0975,.229,4.20,3.30,5.30,1.50,0.700,.600,.5,1.,0.,7.8, 00023 6.3,5.5,3.6,4.9,1.6,4.2,.284,.279,0.} 00024 }; 00025 static const double dite[2][NDIM] = { 00026 {6.00e4,1.94e5,3.31e5,4.32e5,6.28e5,7.50e5,7.73e5,2.62e5,2.50e5,2.57e5,2.84e5, 00027 8.69e5,4.21e5,4.57e5,2.85e5,8.18e5,1.51e6,1.30e6,1.19e6,1.09e6,9.62e5,7.23e5, 00028 4.23e5,5.84e7,6.01e7,0.}, 00029 {8.97e4,1.71e5,2.73e5,3.49e5,5.29e5,4.69e5,6.54e5,1.32e6,1.33e6,1.41e6,1.52e6, 00030 1.51e6,1.82e6,1.84e6,2.31e6,0.,9.98e6,9.98e6,1.00e7,1.10e7,8.34e6,1.01e7,1.07e7, 00031 1.17e7,9.97e6,0.} 00032 }; 00033 static const double ditcrt[NDIM] = {1e11,1e11,1e11,1e11,1e11,1e11,1e11, 00034 1e11,1e11,1e11,1e11,1e11,1e11,1e11,1e11,1e11,1e11,1e11,1e11, 00035 1e11,1e11,1e11,1e11,1e11,1e11,1e11}; 00036 static const double aa[NDIM] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 00037 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; 00038 static const double bb[NDIM] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 00039 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; 00040 static const double cc[NDIM] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 00041 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; 00042 static const double dd[NDIM] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 00043 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; 00044 /* >>chng 03 aug 15, zero out this array. having non-zero values in 00045 * this array had the effect of using 0 for the guess to the low-T dr 00046 * in makerecomb */ 00047 static const double ff[NDIM] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 00048 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; 00049 static const double fyield[NDIM+1] = {.34,.34,.35,.35,.36,.37,.37,.38,.39,.40, 00050 .41,.42,.43,.44,.45,.46,.47,.47,.48,.48,.49,.49,.11,.75,0.,0.,0.}; 00051 00052 long int i, limit, limit2; 00053 00054 DEBUG_ENTRY( "IonIron()" ); 00055 00056 /* iron nelem=26 00057 * 00058 * Fe rates from Woods et al. Ap.J. 249, 399. 00059 * rec from +23, 24 25 from Arnauld et al 85 */ 00060 00061 /* Pequignot and Aldrovandi Ast Ap 161, 169. */ 00062 00063 if( !dense.lgElmtOn[ipIRON] ) 00064 { 00065 fe.fekcld = 0.; 00066 fe.fekhot = 0.; 00067 fe.fegrain = 0.; 00068 return; 00069 } 00070 00071 ion_zero(ipIRON); 00072 00073 ion_photo(ipIRON,false); 00074 /* debugging printout for shell photo rates */ 00075 /*GammaPrtRate( ioQQQ , 0 , ipIRON , true);*/ 00076 /*GammaPrtShells( ipIRON , 13 ); 00077 GammaPrtShells( ipIRON , 12 ); 00078 GammaPrtShells( ipIRON , 0 );*/ 00079 { 00080 /*@-redef@*/ 00081 enum {DEBUG_LOC=false}; 00082 /*@+redef@*/ 00083 if( DEBUG_LOC ) 00084 { 00085 long int iplow , iphi , ipop , ns , ion; 00086 double rate; 00087 ns = 6; 00088 ion = 0; 00089 /* show what some of the ionization sources are */ 00090 iplow = opac.ipElement[ipIRON][ion][ns][0]; 00091 iphi = opac.ipElement[ipIRON][ion][ns][1]; 00092 ipop = opac.ipElement[ipIRON][ion][ns][2]; 00093 rate = ionbal.PhotoRate_Shell[ipIRON][ion][ns][0]; 00094 GammaPrt( iplow , iphi , ipop , ioQQQ, rate , rate*0.1 ); 00095 } 00096 } 00097 00098 /* find collisional ionization rates */ 00099 ion_collis(ipIRON); 00100 00101 /* get recombination coefficients */ 00102 ion_recomb(false,(const double*)dicoef,(const double*)dite,ditcrt,aa,bb,cc,dd,ff,ipIRON); 00103 00104 /* >>chng 06 Feb 28 -- NPA. Add in charge transfer ionization of Mg with S+, Si+, and C+ */ 00105 /* only include this if molecular network is enabled - otherwise no feedback onto 00106 * Si+, S+, and C+ soln */ 00107 if( !co.lgNoCOMole ) 00108 { 00109 ionbal.PhotoRate_Shell[ipIRON][0][6][0] += 00110 CO_findrk("S+,Fe=>S,Fe+")*dense.xIonDense[ipSULPHUR][1] + 00111 CO_findrk("Si+,Fe=>Si,Fe+")*dense.xIonDense[ipSILICON][1] + 00112 CO_findrk("C+,Fe=>C,Fe+")*dense.xIonDense[ipCARBON][1]; 00113 /* CO_sink_rate("Fe");*/ 00114 00115 } 00116 00117 /* solve for ionization balance */ 00118 ion_solver(ipIRON,false); 00119 00120 /* now find total Auger yield of K-alphas 00121 * "cold" iron has M-shell electrons, up to Fe 18 */ 00122 fe.fekcld = 0.; 00123 limit = MIN2(18,dense.IonHigh[ipIRON]); 00124 00125 for( i=dense.IonLow[ipIRON]; i < limit; i++ ) 00126 { 00127 ASSERT( i < NDIM + 1 ); 00128 fe.fekcld += 00129 (realnum)(ionbal.PhotoRate_Shell[ipIRON][i][0][0]*dense.xIonDense[ipIRON][i]* 00130 fyield[i]); 00131 } 00132 00133 /* same sum for hot iron */ 00134 fe.fekhot = 0.; 00135 limit = MAX2(18,dense.IonLow[ipIRON]); 00136 00137 limit2 = MIN2(ipIRON+1,dense.IonHigh[ipIRON]); 00138 ASSERT( limit2 <= LIMELM + 1 ); 00139 00140 for( i=limit; i < limit2; i++ ) 00141 { 00142 ASSERT( i < NDIM + 1 ); 00143 fe.fekhot += 00144 (realnum)(ionbal.PhotoRate_Shell[ipIRON][i][0][0]*dense.xIonDense[ipIRON][i]* 00145 fyield[i]); 00146 } 00147 00148 /* Fe Ka from grains - Fe in grains assumed to be atomic 00149 * gv.elmSumAbund[ipIRON] is number density of iron added over all grain species */ 00150 i = 0; 00151 /* fyield is 0.34 for atomic fe */ 00152 fe.fegrain = ( gv.lgWD01 ) ? 0.f : (realnum)(ionbal.PhotoRate_Shell[ipIRON][i][0][0]*fyield[i]* 00153 gv.elmSumAbund[ipIRON]); 00154 00155 if( trace.lgTrace && trace.lgHeavyBug ) 00156 { 00157 /* print densities of each stage of ionization */ 00158 fprintf( ioQQQ, " Fe" ); 00159 for( i=0; i < 15; i++ ) 00160 { 00161 fprintf( ioQQQ, "\t%.1e", dense.xIonDense[ipIRON][i]/dense.gas_phase[ipIRON] ); 00162 } 00163 fprintf( ioQQQ, "\n" ); 00164 } 00165 00166 if( trace.lgTrace && trace.lgFeBug ) 00167 { 00168 fprintf( ioQQQ, " IonIron-Abund:" ); 00169 for( i=0; i < 27; i++ ) 00170 { 00171 fprintf( ioQQQ, "%10.2e", dense.xIonDense[ipIRON][i] ); 00172 } 00173 fprintf( ioQQQ, "\n" ); 00174 00175 fprintf( ioQQQ, " IonIron - Ka(Auger)%10.2e\n", fe.fekcld + 00176 fe.fekhot ); 00177 } 00178 return; 00179 }