00001 /* This file is part of Cloudy and is copyright (C)1978-2007 by Gary J. Ferland 00002 * 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 "grainvar.h" 00010 #include "hmi.h" 00011 #include "mole.h" 00012 #include "secondaries.h" 00013 #include "iso.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]; 00039 atmdat.HCharExcIonTotal += one; 00040 } 00041 } 00042 00043 /* >>chng 01 may 07, add this in */ 00044 /* charge transfer recombination of hydrogen, 00045 * which is ionization of the heavy element */ 00046 atmdat.HCharExcRecTotal = 0.; 00047 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem) 00048 { 00049 /* this is ion on the abundances scale, 1 is atom, so goes up to nelem+1, 00050 * for helium nelem=1, ion must go up to 3 */ 00051 for( ion=0; ion<=nelem; ++ion ) 00052 { 00053 /* option to skip Oo => O+ */ 00054 /* >>chng 05 aug 17, add logic on !ionbal.lgHO_ct_chem 00055 * this is flag that is true when H O charge transfer is done in 00056 * chemistry, not here. set false with set HO charge transfer ionization 00057 * command, so we do it here - introduced by NA for the Leiden hacks */ 00058 if( (nelem == ipOXYGEN) && (ion == 0) && ionbal.lgHO_ct_chem ) 00059 continue; 00060 /* charge transfer ionization of H, recombination for other species */ 00061 atmdat.HCharExcRecTotal += 00062 atmdat.HCharExcIonOf[nelem][ion]*dense.xIonDense[nelem][ion]; 00063 } 00064 } 00065 00066 /* >>logic checked 02 may 02, think it's right */ 00067 /* add on charge transfer ionization of helium, 00068 * recombination for other element is ionization of helium */ 00069 atmdat.HeCharExcIonTotal = 0.; 00070 /* loop up from Li, the next heavier element */ 00071 for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem) 00072 { 00073 double hold_one = 0.; 00074 /* check that array bounds not exceeded */ 00075 /* ion on the C scale, 0 is atom, so goes up to nelem+1, 00076 * for helium nelem=1, ions go from 1 to 2 */ 00077 for( ion=1; ion<=nelem+1; ++ion ) 00078 { 00079 /* charge transfer ionization of He, recombination for other species */ 00080 hold_one = atmdat.HeCharExcRecTo[nelem][ion-1]*dense.xIonDense[nelem][ion]; 00081 atmdat.HeCharExcIonTotal += hold_one; 00082 } 00083 } 00084 00085 /* >>chng 04 jul 02, 00086 * add on charge transfer reactions of He-H */ 00087 atmdat.HeCharExcIonTotal += atmdat.HCharExcIonOf[ipHELIUM][0]*dense.xIonDense[ipHYDROGEN][1]; 00088 00089 /* charge transfer recombination of helium, 00090 * which is ionization of the heavy element */ 00091 atmdat.HeCharExcRecTotal = 0.; 00092 for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem) 00093 { 00094 /* this is ion on the physics scale, 1 is atom, so goes up to nelem+1, 00095 * for helium nelem=1, ion must go up to 3 */ 00096 for( ion=0; ion<=nelem; ++ion ) 00097 { 00098 /* charge transfer recombination of He, ionization for other species */ 00099 atmdat.HeCharExcRecTotal += 00100 atmdat.HeCharExcIonOf[nelem][ion]*dense.xIonDense[nelem][ion]; 00101 } 00102 } 00103 /* >>chng 04 jul 02 00104 * Add on charge transfer reactions of He+ +H0 -> He0 + H+ */ 00105 atmdat.HeCharExcRecTotal += atmdat.HCharExcRecTo[ipHELIUM][0]*dense.xIonDense[ipHYDROGEN][0]; 00106 00107 return; 00108 } 00109 00110 /*iso_ionize_recombine find state specific creation and destruction rates for iso sequences */ 00111 void iso_ionize_recombine( 00112 /* iso sequence, hydrogen or helium for now */ 00113 long ipISO , 00114 /* the chemical element, 0 for hydrogen */ 00115 long int nelem ) 00116 { 00117 long int level, 00118 ion; 00119 double Recom3Body, 00120 sum; 00121 00122 DEBUG_ENTRY( "iso_ionize_recombine()" ); 00123 00124 /* get total charge transfer ionization rate if this is hydrogen itself, 00125 * routine also does CT for helium but all models must have hydrogen, so 00126 * he is done at this stage */ 00127 if( ipISO==ipH_LIKE && nelem == ipHYDROGEN ) 00128 ChargeTransferUpdate(); 00129 00130 /* find recombination and ionization elements, will use to get simple estimate 00131 * of ionization ratio below */ 00132 /* >>chng 06 jul 20, level should go to numLevels_local instead of numLevels_max */ 00133 for( level=ipH1s; level< iso.numLevels_local[ipISO][nelem]; ++level) 00134 { 00135 /* all process moving level to continuum, units s-1 */ 00136 iso.RateLevel2Cont[ipISO][nelem][level] = iso.gamnc[ipISO][nelem][level] + 00137 iso.ColIoniz[ipISO][nelem][level]* dense.EdenHCorr + 00138 /* >>chng 04 nov 30, also turn off cr ionization with atom h-like collisions off */ 00139 /* >>chng 05 aug 17, do not turn off cosmic rays with no collisional ionization 00140 * simply leave out crs if this is desired - this allows leiden hacks to work */ 00141 secondaries.csupra[nelem][nelem-ipISO]/* *iso.lgColl_ionize[ipISO]*/; 00142 00143 /* all processes from continuum to level n, units s-1 */ 00144 iso.RateCont2Level[ipISO][nelem][level] = ( 00145 /* radiative recombination */ 00146 iso.RadRecomb[ipISO][nelem][level][ipRecRad]* 00147 iso.RadRecomb[ipISO][nelem][level][ipRecNetEsc] + 00148 00149 /* induced recombination */ 00150 iso.RecomInducRate[ipISO][nelem][level]*iso.PopLTE[ipISO][nelem][level] + 00151 00152 /* collisional or three body recombination */ 00153 /* PopLTE(level,nelem) is only LTE pop when mult by n_e n_H */ 00154 iso.ColIoniz[ipISO][nelem][level]*dense.EdenHCorr*iso.PopLTE[ipISO][nelem][level] 00155 ) * dense.eden; 00156 } 00157 00158 /* all following go into or out of ground state */ 00159 level = 0; 00160 00161 /* grain charge transfer recombination and ionization, 00162 * assume goes into and comes from ground state */ 00163 /* first grain surface losses - both to lower and higher stages of ionization */ 00164 sum = 0.; 00165 for(ion=0; ion<nelem+1; ++ion ) 00166 if( ion!=nelem-ipISO ) 00167 sum += gv.GrainChTrRate[nelem][nelem-ipISO][ion]; 00168 /*>>chng 04 sep 08, replaced grain ct rate from old sum with current sum */ 00169 iso.RateLevel2Cont[ipISO][nelem][level] += sum; 00170 00171 /* >>chng 04 apr 10, add formally correct grain surface ionization and recombination - 00172 * this HAD NOT been included for H-like species */ 00173 /* >>chng 04 sep 08, use explicit term from higher stage of ionization - note that added 00174 * GrainCreat[nelem][nelem-ipISO] for case of atomic He incorrectly added two-stage recombination 00175 * to existing matrix logic */ 00176 /*iso.RateCont2Level[ipISO][nelem][level] += ionbal.GrainCreat[nelem][nelem-ipISO];*/ 00177 iso.RateCont2Level[ipISO][nelem][level] += gv.GrainChTrRate[nelem][nelem+1-ipISO][nelem-ipISO]; 00178 00179 /* now charge transfer - all into/from ground, two cases, H and not H */ 00180 if( nelem==ipHYDROGEN ) 00181 { 00182 /* charge transfer, hydrogen onto everything else */ 00183 /* charge exchange ionization */ 00184 iso.RateLevel2Cont[ipISO][nelem][ipH1s] += atmdat.HCharExcIonTotal; 00185 /* charge exchange recombination */ 00186 iso.RateCont2Level[ipISO][nelem][ipH1s] += atmdat.HCharExcRecTotal; 00187 } 00188 else if( nelem==ipHELIUM && ipISO==ipHE_LIKE ) 00189 { 00190 /* add CO here, only for he itself, 00191 * destruction of CO but recombination for he */ 00192 /* >>chng 04 feb 14, update date to 1.4ee-9 from 1.1e-9, 00193 * make var so auto keep parallel with co.c */ 00194 00195 /* The following terms in co.c destroy He+ and thereby need to be included 00196 * here. Also included in the sum of co.hep_destroy are the contributions 00197 * to the recombination due to hmole_step.*/ 00198 00199 co.hep_destroy = (float) CO_sink_rate("He+") + 00200 hmi.H2_total*hmi.rheph2hpheh + hmi.heph2heh2p*hmi.H2_total; 00201 00202 iso.RateCont2Level[ipISO][nelem][ipHe1s1S] += co.hep_destroy; 00203 /*fprintf(ioQQQ,"DEBUG co hep destroy %e\n",co.hep_destroy );*/ 00204 00205 /* this is ioniz of He due to ct with all other gas constituents */ 00206 iso.RateLevel2Cont[ipISO][nelem][ipHe1s1S] += atmdat.HeCharExcIonTotal; 00207 /* this is recom of He due to ct with all other gas constituents */ 00208 iso.RateCont2Level[ipISO][nelem][ipHe1s1S] += atmdat.HeCharExcRecTotal; 00209 } 00210 else 00211 { 00212 iso.RateCont2Level[ipISO][nelem][level] += 00213 atmdat.HeCharExcRecTo[nelem][nelem-ipISO]*iso.Pop2Ion[ipHE_LIKE][ipHELIUM][ipHe1s1S]*dense.xIonDense[ipHELIUM][1] + 00214 atmdat.HCharExcRecTo[nelem][nelem-ipISO]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]*dense.xIonDense[ipHYDROGEN][1]; 00215 00216 iso.RateLevel2Cont[ipISO][nelem][level] += 00217 atmdat.HeCharExcIonOf[nelem][nelem-ipISO]*dense.xIonDense[ipHELIUM][1] + 00218 atmdat.HCharExcIonOf[nelem][nelem-ipISO]*dense.xIonDense[ipHYDROGEN][1]; 00219 } 00220 00221 00222 /* now create sums of recombination and ionization rates for ISO species */ 00223 ionbal.RateRecomTot[nelem][nelem-ipISO] = 0.; 00224 ionbal.RR_rate_coef_used[nelem][nelem-ipISO] = 0.; 00225 Recom3Body = 0.; 00226 /* >>chng 06 jul 20, level should go to numLevels_local instead of numLevels_max */ 00227 for( level=ipH1s; level< iso.numLevels_local[ipISO][nelem]; ++level) 00228 { 00229 00230 /* units of ionbal.RateRecomTot are s-1, 00231 * equivalent ionization term is done after level populations are known */ 00232 ionbal.RateRecomTot[nelem][nelem-ipISO] += iso.RateCont2Level[ipISO][nelem][level]; 00233 00234 /* just the radiative recombination rate coef, cm3 s-1 */ 00235 ionbal.RR_rate_coef_used[nelem][nelem-ipISO] += iso.RadRecomb[ipISO][nelem][level][ipRecRad]* 00236 iso.RadRecomb[ipISO][nelem][level][ipRecNetEsc]; 00237 /*if( nelem==ipHELIUM && ipISO==ipH_LIKE ) 00238 fprintf(ioQQQ,"DEBUG rad rec %li %.4e %.4e\n", 00239 level, 00240 iso.RadRecomb[ipISO][nelem][level][ipRecRad], 00241 iso.RadRecomb[ipISO][nelem][level][ipRecNetEsc]);*/ 00242 00243 /* >>chng 05 jul 11, from > to >=, some very opt thick sims did block escape to zero */ 00244 ASSERT( ionbal.RR_rate_coef_used[nelem][nelem-ipISO]>= 0. ); 00245 00246 /* this is three-body recombination rate coef by itself - 00247 * need factor of eden to become rate */ 00248 Recom3Body += iso.ColIoniz[ipISO][nelem][level]*dense.EdenHCorr*iso.PopLTE[ipISO][nelem][level]; 00249 } 00250 /*if( nelem==ipHELIUM && ipISO==ipH_LIKE ) 00251 fprintf(ioQQQ,"DEBUG rad rec %.4e\n",ionbal.RR_rate_coef_used[nelem][nelem-ipISO]);*/ 00252 00253 /* fraction of total recombs due to three body - when most are due to this 00254 * small changes in temperature can produce large changes in rec coef, 00255 * and so in ionization */ 00256 iso.RecomCollisFrac[ipISO][nelem] = (float)(Recom3Body* dense.eden / ionbal.RateRecomTot[nelem][nelem-ipISO] ); 00257 00258 /* get simple estimate of level of ionization */ 00259 if( ionbal.RateRecomTot[nelem][nelem-ipISO] > 0. ) 00260 { 00261 # if 0 00262 /* >>chng 06 jul 23, this was a temporary hack to force the lot T solver, which forces the 00263 * atom to agree with iso.xIonSimple, to synch up with full solln from matrix - this caused 00264 * many problems since using ionbal.RateIonizTot includes excited level ionization processes 00265 * so matrix soln was called rather than simple. this cause neg pops in ~half a dozen sims */ 00266 if( ionbal.RateIonizTot[nelem][nelem-ipISO] >= 0. ) 00267 { 00268 /* >>chng 06 jul 21, use effective rate, which includes ioniz from excited states, if defined */ 00269 iso.xIonSimple[ipISO][nelem] = ionbal.RateIonizTot[nelem][nelem-ipISO]/ionbal.RateRecomTot[nelem][nelem-ipISO]; 00270 } 00271 else 00272 # endif 00273 { 00274 iso.xIonSimple[ipISO][nelem] = iso.RateLevel2Cont[ipISO][nelem][ipH1s]/ionbal.RateRecomTot[nelem][nelem-ipISO]; 00275 } 00276 /*fprintf(ioQQQ,"DEBUG iso ion rec %li %li ion %.2e ioniztot %.2e \n" , 00277 ipISO, nelem , 00278 iso.RateLevel2Cont[ipISO][nelem][ipH1s] , 00279 ionbal.RateIonizTot[nelem][nelem-ipISO] );*/ 00280 } 00281 else 00282 { 00283 iso.xIonSimple[ipISO][nelem] = 0.; 00284 } 00285 00286 DEBUG_EXIT( "iso_ionize_recombine()" ); 00287 00288 return; 00289 } 00290 /*lint +e778 constant expression evaluates to zero - in HMRATE macro */