00001 /* This file is part of Cloudy and is copyright (C)1978-2013 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 "phycon.h" 00015 #include "rfield.h" 00016 #include "trace.h" 00017 #include "taulines.h" 00018 /*lint -e778 constant expression evaluates to zero - in HMRATE macro */ 00019 /*iso_charge_transfer_update update rate of ct ionization and recombination for H atoms */ 00020 void iso_charge_transfer_update(long nelem1) 00021 { 00022 DEBUG_ENTRY( "iso_charge_transfer_update()" ); 00023 00024 if( nelem1==ipHYDROGEN ) 00025 { 00026 atmdat.HCharExcIonTotal = 0.; 00027 atmdat.HCharExcRecTotal = 0.; 00028 // find total charge transfer rate coefficients of hydrogen 00029 for( long nelem=ipHELIUM; nelem<LIMELM; ++nelem) 00030 { 00031 for( long ion=0; ion<=nelem; ++ion ) 00032 { 00033 // charge transfer ionization of H, recombination for other species, cm-3 s-1 00034 atmdat.HCharExcIonTotal += atmdat.HCharExcRecTo[nelem][ion]*dense.xIonDense[nelem][ion+1]; 00035 // charge transfer recombination of H, ionization for other species, cm-3 s-1 00036 atmdat.HCharExcRecTotal += atmdat.HCharExcIonOf[nelem][ion]*dense.xIonDense[nelem][ion]; 00037 } 00038 } 00039 } 00040 else if( nelem1==ipHELIUM ) 00041 { 00042 atmdat.HeCharExcIonTotal = 0.; 00043 atmdat.HeCharExcRecTotal = 0.; 00044 // find total charge transfer rate coefficients of helium 00045 for( long nelem=ipLITHIUM; nelem<LIMELM; ++nelem) 00046 { 00047 for( long ion=0; ion<=nelem; ++ion ) 00048 { 00049 // charge transfer ionization of He, recombination for other species, cm-3 s-1 00050 atmdat.HeCharExcIonTotal += atmdat.HeCharExcRecTo[nelem][ion]*dense.xIonDense[nelem][ion+1]; 00051 // charge transfer recombination of He, ionization for other species, cm-3 s-1 00052 atmdat.HeCharExcRecTotal += atmdat.HeCharExcIonOf[nelem][ion]*dense.xIonDense[nelem][ion]; 00053 } 00054 } 00055 // add on charge transfer reactions between He and H 00056 atmdat.HeCharExcIonTotal += atmdat.HCharExcIonOf[ipHELIUM][0]*dense.xIonDense[ipHYDROGEN][1]; 00057 atmdat.HeCharExcRecTotal += atmdat.HCharExcRecTo[ipHELIUM][0]*dense.xIonDense[ipHYDROGEN][0]; 00058 } 00059 00060 return; 00061 } 00062 00063 /*iso_ionize_recombine find state specific creation and destruction rates for iso sequences */ 00064 void iso_ionize_recombine( 00065 /* iso sequence, hydrogen or helium for now */ 00066 long ipISO , 00067 /* the chemical element, 0 for hydrogen */ 00068 long int nelem ) 00069 { 00070 long int level; 00071 double Recom3Body; 00072 00073 DEBUG_ENTRY( "iso_ionize_recombine()" ); 00074 00075 ASSERT( ipISO >= 0 && ipISO < NISO ); 00076 ASSERT( nelem >= 0 && nelem < LIMELM ); 00077 00078 /* find recombination and ionization elements, will use to get simple estimate 00079 * of ionization ratio below */ 00080 /* >>chng 06 jul 20, level should go to numLevels_local instead of numLevels_max */ 00081 00082 fixit(); /* must apply error to iso.gamnc */ 00083 00084 for( level=ipH1s; level< iso_sp[ipISO][nelem].numLevels_local; ++level) 00085 { 00086 // This term is intended to capture LTE in DR-dominated plasmas. 00087 // It is scaled by Occ num / Occ num at STE to ensure correct limits w.r.t. ambient radiation. 00088 double DR_reverse = 0.; 00089 if( ipISO > ipH_LIKE && iso_sp[ipISO][nelem].fb[level].PopLTE > SMALLFLOAT ) 00090 { 00091 long indexIP = iso_sp[ipISO][nelem].fb[level].ipIsoLevNIonCon-1; 00092 double OccNum = rfield.OccNumbIncidCont[indexIP] + rfield.OccNumbContEmitOut[indexIP]; 00093 double OccNumBB = 1./(exp( iso_sp[ipISO][nelem].fb[level].xIsoLevNIonRyd / phycon.te_ryd ) - 1. ); 00094 double RelOccNum = MIN2( 1., OccNum / OccNumBB ); 00095 //fprintf( ioQQQ, "DEBUGGG DR_reverse %li\t%li\t%e\t%e\t%e\t%e\t%e\t%e\n", level, indexIP, RelOccNum, OccNum, OccNumBB, 00096 // rfield.anu[indexIP], iso_sp[ipISO][nelem].fb[level].xIsoLevNIonRyd, rfield.anu[indexIP+1] ); 00097 DR_reverse = iso_sp[ipISO][nelem].fb[level].DielecRecomb * RelOccNum / 00098 iso_sp[ipISO][nelem].fb[level].PopLTE; 00099 } 00100 00101 /* all process moving level to continuum, units s-1 */ 00102 iso_sp[ipISO][nelem].fb[level].RateLevel2Cont = iso_sp[ipISO][nelem].fb[level].gamnc + 00103 iso_sp[ipISO][nelem].fb[level].ColIoniz* dense.EdenHCorr + 00104 DR_reverse + 00105 secondaries.csupra[nelem][nelem-ipISO]*iso_ctrl.lgColl_ionize[ipISO]; 00106 00107 /* all processes from continuum to level n, units s-1 */ 00108 iso_sp[ipISO][nelem].fb[level].RateCont2Level = ( 00109 /* radiative recombination */ 00110 iso_sp[ipISO][nelem].fb[level].RadRecomb[ipRecRad]* 00111 iso_sp[ipISO][nelem].fb[level].RadRecomb[ipRecNetEsc] + 00112 00113 /* dielectronic recombination */ 00114 iso_sp[ipISO][nelem].fb[level].DielecRecomb + 00115 00116 /* induced recombination */ 00117 iso_sp[ipISO][nelem].fb[level].RecomInducRate*iso_sp[ipISO][nelem].fb[level].PopLTE + 00118 00119 /* collisional or three body recombination */ 00120 /* PopLTE(level,nelem) is only LTE pop when mult by n_e n_H */ 00121 iso_sp[ipISO][nelem].fb[level].ColIoniz*dense.EdenHCorr*iso_sp[ipISO][nelem].fb[level].PopLTE 00122 ) * dense.eden; 00123 00124 ASSERT( iso_sp[ipISO][nelem].fb[level].RadRecomb[ipRecRad] >= 0. ); 00125 ASSERT( iso_sp[ipISO][nelem].fb[level].RadRecomb[ipRecNetEsc] >= 0. ); 00126 ASSERT( iso_sp[ipISO][nelem].fb[level].DielecRecomb >= 0. ); 00127 ASSERT( iso_sp[ipISO][nelem].fb[level].RecomInducRate >= 0. ); 00128 ASSERT( iso_sp[ipISO][nelem].fb[level].PopLTE >= 0. ); 00129 ASSERT( iso_sp[ipISO][nelem].fb[level].ColIoniz >= 0. ); 00130 ASSERT( iso_sp[ipISO][nelem].fb[level].RateCont2Level >= 0. ); 00131 00132 if( iso_ctrl.lgRandErrGen[ipISO] ) 00133 { 00134 iso_sp[ipISO][nelem].fb[level].RateCont2Level *= 00135 iso_sp[ipISO][nelem].ex[ iso_sp[ipISO][nelem].numLevels_max ][level].ErrorFactor[IPRAD]; 00136 } 00137 } 00138 00139 /* now create sums of recombination and ionization rates for ISO species */ 00140 ionbal.RateRecomIso[nelem][ipISO] = 0.; 00141 ionbal.RR_rate_coef_used[nelem][nelem-ipISO] = 0.; 00142 Recom3Body = 0.; 00143 /* >>chng 06 jul 20, level should go to numLevels_local instead of numLevels_max */ 00144 for( level=0; level< iso_sp[ipISO][nelem].numLevels_local; ++level) 00145 { 00146 00147 /* units of ionbal.RateRecomTot are s-1, 00148 * equivalent ionization term is done after level populations are known */ 00149 ionbal.RateRecomIso[nelem][ipISO] += iso_sp[ipISO][nelem].fb[level].RateCont2Level; 00150 00151 /* just the radiative recombination rate coef, cm3 s-1 */ 00152 ionbal.RR_rate_coef_used[nelem][nelem-ipISO] += iso_sp[ipISO][nelem].fb[level].RadRecomb[ipRecRad]* 00153 iso_sp[ipISO][nelem].fb[level].RadRecomb[ipRecNetEsc]; 00154 00155 /* >>chng 05 jul 11, from > to >=, some very opt thick sims did block escape to zero */ 00156 ASSERT( ionbal.RR_rate_coef_used[nelem][nelem-ipISO]>= 0. ); 00157 00158 /* this is three-body recombination rate coef by itself - 00159 * need factor of eden to become rate */ 00160 Recom3Body += iso_sp[ipISO][nelem].fb[level].ColIoniz*dense.EdenHCorr*iso_sp[ipISO][nelem].fb[level].PopLTE; 00161 } 00162 00163 /* fraction of total recombinations due to three body - when most are due to this 00164 * small changes in temperature can produce large changes in recombination coefficient, 00165 * and so in ionization */ 00166 iso_sp[ipISO][nelem].RecomCollisFrac = Recom3Body* dense.eden / ionbal.RateRecomIso[nelem][ipISO]; 00167 00168 /* very first pass through here rate RateIoniz not yet evaluated */ 00169 if( conv.nTotalIoniz==0 ) 00170 ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] = iso_sp[ipISO][nelem].fb[0].RateLevel2Cont; 00171 00172 /* get simple estimate of atom to ion ratio */ 00173 if( ionbal.RateRecomIso[nelem][ipISO] > 0. ) 00174 { 00175 iso_sp[ipISO][nelem].xIonSimple = ionbal.RateIonizTot(nelem,nelem-ipISO)/ionbal.RateRecomIso[nelem][ipISO]; 00176 } 00177 else 00178 { 00179 iso_sp[ipISO][nelem].xIonSimple = 0.; 00180 } 00181 00182 if( trace.lgTrace && (trace.lgHBug||trace.lgHeBug) ) 00183 { 00184 fprintf( ioQQQ, " iso_ionize_recombine iso=%2ld Z=%2ld Level2Cont[0] %10.2e RateRecomTot %10.2e xIonSimple %10.2e\n", 00185 ipISO, nelem, iso_sp[ipISO][nelem].fb[0].RateLevel2Cont, ionbal.RateRecomIso[nelem][ipISO], iso_sp[ipISO][nelem].xIonSimple ); 00186 } 00187 00188 return; 00189 } 00190 /*lint +e778 constant expression evaluates to zero - in HMRATE macro */