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