00001 /* This file is part of Cloudy and is copyright (C)1978-2010 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*iso_ionize_recombine find state specific creation and destruction rates for iso sequences */ 00004 /*ChargeTransferUpdate update rate of ct ionization and recombination for H atoms */ 00005 #include "cddefines.h" 00006 #include "ionbal.h" 00007 #include "conv.h" 00008 #include "atmdat.h" 00009 #include "dense.h" 00010 #include "hmi.h" 00011 #include "mole.h" 00012 #include "secondaries.h" 00013 #include "iso.h" 00014 #include "trace.h" 00015 /*lint -e778 constant expression evaluates to zero - in HMRATE macro */ 00016 /*iso_charge_transfer_update update rate of ct ionization and recombination for H atoms */ 00017 void iso_charge_transfer_update(void) 00018 { 00019 long int ion; 00020 long int nelem; 00021 00022 DEBUG_ENTRY( "iso_charge_transfer_update()" ); 00023 00024 /* find total rate for charge transfer ionization of hydrogen, 00025 * recombination for other element is ionization of hydrogen */ 00026 atmdat.HCharExcIonTotal = 0.; 00027 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem) 00028 { 00029 /* ion on the C scale, 0 is atom, so goes up to nelem+1, 00030 * for helium nelem=1, ions go from 1 to 2 */ 00031 for( ion=1; ion<=nelem+1; ++ion ) 00032 { 00033 double one; 00034 /* we intentionally skip CT with O+ since this is in hmole */ 00035 /* >>chng 05 aug 17, add logic on !ionbal.lgHO_ct_chem 00036 * this is flag that is true when H O charge transfer is done in 00037 * chemistry, not here. set false with set HO charge transfer ionization 00038 * command, so we do it here - introduced by NA for the Leiden hacks */ 00039 if( (nelem == ipOXYGEN) && (ion == 1) && ionbal.lgHO_ct_chem ) 00040 continue; 00041 /* charge transfer ionization of H, recombination for other species */ 00042 one = atmdat.HCharExcRecTo[nelem][ion-1]*dense.xIonDense[nelem][ion]; //units s-1 00043 /* this, when multiplied by dense.xIonDense[ipHYDROGEN][0], is the charge transfer 00044 * ionization RATE of neutral hydrogen, units cm-3 s-1 */ 00045 atmdat.HCharExcIonTotal += one; 00046 } 00047 } 00048 00049 /* >>chng 01 may 07, add this in */ 00050 /* charge transfer recombination of hydrogen, 00051 * which is ionization of the heavy element */ 00052 atmdat.HCharExcRecTotal = 0.; 00053 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem) 00054 { 00055 /* this is ion on the abundances scale, 1 is atom, so goes up to nelem+1, 00056 * for helium nelem=1, ion must go up to 3 */ 00057 for( ion=0; ion<=nelem; ++ion ) 00058 { 00059 /* option to skip Oo => O+ */ 00060 /* >>chng 05 aug 17, add logic on !ionbal.lgHO_ct_chem 00061 * this is flag that is true when H O charge transfer is done in 00062 * chemistry, not here. set false with set HO charge transfer ionization 00063 * command, so we do it here - introduced by NA for the Leiden hacks */ 00064 if( (nelem == ipOXYGEN) && (ion == 0) && ionbal.lgHO_ct_chem ) 00065 continue; 00066 /* charge transfer ionization of H, recombination for other species */ 00067 /* this, when multiplied by dense.xIonDense[ipHYDROGEN][1], is the charge transfer 00068 * recombination RATE of neutral hydrogen, units cm-3 s-1 */ 00069 atmdat.HCharExcRecTotal += atmdat.HCharExcIonOf[nelem][ion]*dense.xIonDense[nelem][ion]; 00070 } 00071 } 00072 00073 /* >>logic checked 02 may 02, think it's right */ 00074 /* add on charge transfer ionization of helium, 00075 * recombination for other element is ionization of helium */ 00076 atmdat.HeCharExcIonTotal = 0.; 00077 /* loop up from Li, the next heavier element */ 00078 for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem) 00079 { 00080 double hold_one = 0.; 00081 /* check that array bounds not exceeded */ 00082 /* ion on the C scale, 0 is atom, so goes up to nelem+1, 00083 * for helium nelem=1, ions go from 1 to 2 */ 00084 for( ion=1; ion<=nelem+1; ++ion ) 00085 { 00086 /* charge transfer ionization of He, recombination for other species */ 00087 hold_one = atmdat.HeCharExcRecTo[nelem][ion-1]*dense.xIonDense[nelem][ion]; 00088 atmdat.HeCharExcIonTotal += hold_one; 00089 } 00090 } 00091 00092 /* >>chng 04 jul 02, 00093 * add on charge transfer reactions of He-H */ 00094 atmdat.HeCharExcIonTotal += atmdat.HCharExcIonOf[ipHELIUM][0]*dense.xIonDense[ipHYDROGEN][1]; 00095 00096 /* charge transfer recombination of helium, 00097 * which is ionization of the heavy element */ 00098 atmdat.HeCharExcRecTotal = 0.; 00099 for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem) 00100 { 00101 /* this is ion on the physics scale, 1 is atom, so goes up to nelem+1, 00102 * for helium nelem=1, ion must go up to 3 */ 00103 for( ion=0; ion<=nelem; ++ion ) 00104 { 00105 /* charge transfer recombination of He, ionization for other species */ 00106 atmdat.HeCharExcRecTotal += 00107 atmdat.HeCharExcIonOf[nelem][ion]*dense.xIonDense[nelem][ion]; 00108 } 00109 } 00110 /* >>chng 04 jul 02 00111 * Add on charge transfer reactions of He+ +H0 -> He0 + H+ */ 00112 atmdat.HeCharExcRecTotal += atmdat.HCharExcRecTo[ipHELIUM][0]*dense.xIonDense[ipHYDROGEN][0]; 00113 00114 return; 00115 } 00116 00117 /*iso_ionize_recombine find state specific creation and destruction rates for iso sequences */ 00118 void iso_ionize_recombine( 00119 /* iso sequence, hydrogen or helium for now */ 00120 long ipISO , 00121 /* the chemical element, 0 for hydrogen */ 00122 long int nelem ) 00123 { 00124 long int level; 00125 double Recom3Body; 00126 00127 DEBUG_ENTRY( "iso_ionize_recombine()" ); 00128 00129 ASSERT( ipISO >= 0 && ipISO < NISO ); 00130 ASSERT( nelem >= 0 && nelem < LIMELM ); 00131 00132 /* find recombination and ionization elements, will use to get simple estimate 00133 * of ionization ratio below */ 00134 /* >>chng 06 jul 20, level should go to numLevels_local instead of numLevels_max */ 00135 00136 fixit(); /* must apply error to iso.gamnc */ 00137 00138 for( level=ipH1s; level< iso.numLevels_local[ipISO][nelem]; ++level) 00139 { 00140 /* all process moving level to continuum, units s-1 */ 00141 iso.RateLevel2Cont[ipISO][nelem][level] = iso.gamnc[ipISO][nelem][level] + 00142 iso.ColIoniz[ipISO][nelem][level]* dense.EdenHCorr + 00143 secondaries.csupra[nelem][nelem-ipISO]*iso.lgColl_ionize[ipISO]; 00144 00145 /* all processes from continuum to level n, units s-1 */ 00146 iso.RateCont2Level[ipISO][nelem][level] = ( 00147 /* radiative recombination */ 00148 iso.RadRecomb[ipISO][nelem][level][ipRecRad]* 00149 iso.RadRecomb[ipISO][nelem][level][ipRecNetEsc] + 00150 00151 /* dielectronic recombination */ 00152 iso.DielecRecomb[ipISO][nelem][level] + 00153 00154 /* induced recombination */ 00155 iso.RecomInducRate[ipISO][nelem][level]*iso.PopLTE[ipISO][nelem][level] + 00156 00157 /* collisional or three body recombination */ 00158 /* PopLTE(level,nelem) is only LTE pop when mult by n_e n_H */ 00159 iso.ColIoniz[ipISO][nelem][level]*dense.EdenHCorr*iso.PopLTE[ipISO][nelem][level] 00160 ) * dense.eden; 00161 00162 if( iso.lgRandErrGen[ipISO] ) 00163 { 00164 iso.RateCont2Level[ipISO][nelem][level] *= 00165 iso.ErrorFactor[ipISO][nelem][ iso.numLevels_max[ipISO][nelem] ][level][IPRAD]; 00166 } 00167 } 00168 00169 /* now charge transfer - all into/from ground, two cases, H and not H */ 00170 if( nelem==ipHYDROGEN ) 00171 { 00172 /* charge transfer, hydrogen onto everything else */ 00173 /* charge exchange ionization, these come out of every level */ 00174 for( level=0; level< iso.numLevels_local[ipISO][nelem]; ++level) 00175 { 00176 iso.RateLevel2Cont[ipISO][nelem][level] += atmdat.HCharExcIonTotal; 00177 } 00178 /* charge exchange recombination */ 00179 iso.RateCont2Level[ipISO][nelem][0] += atmdat.HCharExcRecTotal; 00180 } 00181 else if( nelem==ipHELIUM && ipISO==ipHE_LIKE ) 00182 { 00183 /* add CO here, only for he itself, 00184 * destruction of CO but recombination for he */ 00185 00186 /* The following terms destroy He+ and thereby need to be included 00187 * here. Also included are the contributions 00188 * to the recombination due to hmole_step.*/ 00189 00190 fixit(); // should these chemistry rates be accounted for in mole.source and mole.sink? 00191 double COsource = CO_source_rate("He+"); 00192 for( level=0; level< iso.numLevels_local[ipISO][nelem]; ++level) 00193 { 00194 iso.RateLevel2Cont[ipISO][nelem][level] += COsource; 00195 fixit(); //need reverse rates for hmi rates just below 00196 } 00197 iso.RateCont2Level[ipISO][nelem][0] += CO_sink_rate("He+") + 00198 hmi.H2_total*hmi.rheph2hpheh + hmi.heph2heh2p*hmi.H2_total; 00199 00200 /* this is ioniz of He due to ct with all other gas constituents */ 00201 for( level=0; level< iso.numLevels_local[ipISO][nelem]; ++level) 00202 { 00203 iso.RateLevel2Cont[ipISO][nelem][level] += atmdat.HeCharExcIonTotal; 00204 } 00205 /* this is recom of He due to ct with all other gas constituents */ 00206 iso.RateCont2Level[ipISO][nelem][0] += atmdat.HeCharExcRecTotal; 00207 } 00208 else 00209 { 00210 for( level=0; level< iso.numLevels_local[ipISO][nelem]; ++level) 00211 { 00212 iso.RateLevel2Cont[ipISO][nelem][level] += 00213 atmdat.HeCharExcIonOf[nelem][nelem-ipISO]*dense.xIonDense[ipHELIUM][1] + 00214 atmdat.HCharExcIonOf[nelem][nelem-ipISO]*dense.xIonDense[ipHYDROGEN][1]; 00215 } 00216 00217 iso.RateCont2Level[ipISO][nelem][0] += 00218 atmdat.HeCharExcRecTo[nelem][nelem-ipISO]*dense.xIonDense[ipHELIUM][0] + 00219 atmdat.HCharExcRecTo[nelem][nelem-ipISO]*dense.xIonDense[ipHYDROGEN][0]; 00220 } 00221 00222 00223 /* now create sums of recombination and ionization rates for ISO species */ 00224 ionbal.RateRecomTot[nelem][nelem-ipISO] = 0.; 00225 ionbal.RR_rate_coef_used[nelem][nelem-ipISO] = 0.; 00226 Recom3Body = 0.; 00227 /* >>chng 06 jul 20, level should go to numLevels_local instead of numLevels_max */ 00228 for( level=0; level< iso.numLevels_local[ipISO][nelem]; ++level) 00229 { 00230 00231 /* units of ionbal.RateRecomTot are s-1, 00232 * equivalent ionization term is done after level populations are known */ 00233 ionbal.RateRecomTot[nelem][nelem-ipISO] += iso.RateCont2Level[ipISO][nelem][level]; 00234 00235 /* just the radiative recombination rate coef, cm3 s-1 */ 00236 ionbal.RR_rate_coef_used[nelem][nelem-ipISO] += iso.RadRecomb[ipISO][nelem][level][ipRecRad]* 00237 iso.RadRecomb[ipISO][nelem][level][ipRecNetEsc]; 00238 00239 /* >>chng 05 jul 11, from > to >=, some very opt thick sims did block escape to zero */ 00240 ASSERT( ionbal.RR_rate_coef_used[nelem][nelem-ipISO]>= 0. ); 00241 00242 /* this is three-body recombination rate coef by itself - 00243 * need factor of eden to become rate */ 00244 Recom3Body += iso.ColIoniz[ipISO][nelem][level]*dense.EdenHCorr*iso.PopLTE[ipISO][nelem][level]; 00245 } 00246 00247 /* fraction of total recombinations due to three body - when most are due to this 00248 * small changes in temperature can produce large changes in recombination coefficient, 00249 * and so in ionization */ 00250 iso.RecomCollisFrac[ipISO][nelem] = Recom3Body* dense.eden / ionbal.RateRecomTot[nelem][nelem-ipISO]; 00251 00252 /* very first pass through here rate RateIoniz not yet evaluated */ 00253 if( conv.nTotalIoniz==0 ) 00254 ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] = iso.RateLevel2Cont[ipISO][nelem][0]; 00255 00256 /* get simple estimate of atom to ion ratio */ 00257 if( ionbal.RateRecomTot[nelem][nelem-ipISO] > 0. ) 00258 { 00259 iso.xIonSimple[ipISO][nelem] = ionbal.RateIonizTot(nelem,nelem-ipISO)/ionbal.RateRecomTot[nelem][nelem-ipISO]; 00260 } 00261 else 00262 { 00263 iso.xIonSimple[ipISO][nelem] = 0.; 00264 } 00265 00266 if( trace.lgTrace && (trace.lgHBug||trace.lgHeBug) ) 00267 { 00268 fprintf( ioQQQ, " iso_ionize_recombine iso=%2ld Z=%2ld Level2Cont[0] %10.2e RateRecomTot %10.2e xIonSimple %10.2e\n", 00269 ipISO, nelem, iso.RateLevel2Cont[ipISO][nelem][0], ionbal.RateRecomTot[nelem][nelem-ipISO], iso.xIonSimple[ipISO][nelem] ); 00270 } 00271 00272 return; 00273 } 00274 /*lint +e778 constant expression evaluates to zero - in HMRATE macro */