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 /*HydroT2Low called to do hydrogenic level populations when temp too low for matrix, 00004 * forces total of level populations to add up to iso.xIonSimple[ipISO][nelem], 00005 * so results of this routine will always agree with that value */ 00006 #include "cddefines.h" 00007 #include "taulines.h" 00008 #include "iso.h" 00009 #include "secondaries.h" 00010 #include "ionbal.h" 00011 #include "dense.h" 00012 #include "trace.h" 00013 #include "phycon.h" 00014 #include "hydrogenic.h" 00015 00016 void HydroT2Low( long int ipISO, long int nelem ) 00017 { 00018 long int iup, 00019 level, 00020 low; 00021 double entries, 00022 exits, 00023 sum, 00024 renorm, 00025 collider; 00026 00027 DEBUG_ENTRY( "HydroT2Low()" ); 00028 00029 /* check that we were called with valid charge */ 00030 ASSERT( nelem >= 0); 00031 ASSERT( nelem < LIMELM ); 00032 /* >>chng 05 aug 17, use eden as collider for H and corrected eden for heavier 00033 * nuclei - in hot PDR H is collisionally ionized, but not by other H0 since 00034 * collision is homonuclear */ 00035 if( nelem==ipHYDROGEN ) 00036 { 00037 /* special version for H0 onto H0 */ 00038 collider = dense.EdenHontoHCorr; 00039 } 00040 else 00041 { 00042 collider = dense.EdenHCorr; 00043 } 00044 00045 /* do radiative and downward collisional cascades */ 00046 /* 06 aug 28, from numLevels_max to _local. */ 00047 for( level=iso.numLevels_local[ipISO][nelem]-1; level >= 0; level-- ) 00048 { 00049 /* captures into this level from the continuum, 00050 * radiative rec plus three body rec, units s-1 */ 00051 entries = iso.RateCont2Level[ipISO][nelem][level]; 00052 00053 /* sum over all levels between current level and top of atom */ 00054 /* 06 aug 28, from numLevels_max to _local. */ 00055 for( iup=level + 1; iup < iso.numLevels_local[ipISO][nelem]; iup++ ) 00056 { 00057 /* radiative decays from upper level to here, s-1 */ 00058 double RadDecay; 00059 00060 RadDecay = MAX2( iso.SmallA, 00061 Transitions[ipISO][nelem][iup][level].Emis->Aul* 00062 (Transitions[ipISO][nelem][iup][level].Emis->Pesc + 00063 Transitions[ipISO][nelem][iup][level].Emis->Pelec_esc + 00064 Transitions[ipISO][nelem][iup][level].Emis->Pdest) ); 00065 00066 /* rate upper level (rel to continuum) decays to here, units s-1 */ 00067 entries += StatesElem[ipISO][nelem][iup].Pop*(RadDecay + 00068 Transitions[ipISO][nelem][iup][level].Coll.ColUL*collider); 00069 } 00070 00071 /* now find rate this levels decays to all lower levels and continuum 00072 * continuum ionization rate, s-1, first */ 00073 exits = iso.gamnc[ipISO][nelem][level]; 00074 for( low=ipH1s; low <= (level - 1); low++ ) 00075 { 00076 /* radiative decays to all lower levels - s-1 */ 00077 double RadDecay; 00078 00079 RadDecay = MAX2( iso.SmallA, 00080 Transitions[ipISO][nelem][level][low].Emis->Aul* 00081 (Transitions[ipISO][nelem][level][low].Emis->Pesc + 00082 Transitions[ipISO][nelem][level][low].Emis->Pelec_esc + 00083 Transitions[ipISO][nelem][level][low].Emis->Pdest) ); 00084 00085 /* add rad and collisional decays */ 00086 exits += RadDecay + 00087 Transitions[ipISO][nelem][level][low].Coll.ColUL*collider; 00088 } 00089 if( exits > 1e-25 ) 00090 { 00091 StatesElem[ipISO][nelem][level].Pop = entries/exits; 00092 00093 } 00094 else 00095 { 00096 StatesElem[ipISO][nelem][level].Pop = 0.; 00097 } 00098 } 00099 00100 /* =================================================================== 00101 * 00102 * at this point all destruction processes have been established 00103 * 00104 * ==================================================================== */ 00105 00106 /* do simple sums for atom to ion, then fill in for rest */ 00107 StatesElem[ipISO][nelem][ipH1s].Pop = 00108 ionbal.RateRecomTot[nelem][nelem-ipISO]/ 00109 iso.RateLevel2Cont[ipISO][nelem][ipH1s]; 00110 00111 if( ipISO==ipH_LIKE ) 00112 { 00113 00114 StatesElem[ipISO][nelem][ipH2s].Pop = 00115 (iso.RadRec_caseB[ipISO][nelem]*0.33*dense.eden + 00116 StatesElem[ipISO][nelem][ipH1s].Pop*secondaries.Hx12[ipISO][nelem][ipH2s] + 00117 StatesElem[ipISO][nelem][ipH2p].Pop* 00118 Transitions[ipISO][nelem][ipH2p][ipH2s].Coll.ColUL*collider)/ 00119 (Transitions[ipISO][nelem][ipH2s][ipH1s].Emis->Aul + 00120 Transitions[ipISO][nelem][ipH2s][ipH1s].Coll.ColUL*collider + 00121 Transitions[ipISO][nelem][ipH2p][ipH2s].Coll.ColUL*collider*6. ); 00122 00123 StatesElem[ipISO][nelem][ipH2p].Pop = 00124 (iso.RadRec_caseB[ipISO][nelem]*0.67*dense.eden + 00125 StatesElem[ipISO][nelem][ipH1s].Pop*secondaries.Hx12[ipISO][nelem][ipH2p] + 00126 StatesElem[ipISO][nelem][ipH2s].Pop* 00127 Transitions[ipISO][nelem][ipH2p][ipH2s].Coll.ColUL*collider*6.)/ 00128 (Transitions[ipISO][nelem][ipH2p][ipH1s].Emis->Aul* 00129 (Transitions[ipISO][nelem][ipH2p][ipH1s].Emis->Pesc + 00130 Transitions[ipISO][nelem][ipH2p][ipH1s].Emis->Pelec_esc + 00131 Transitions[ipISO][nelem][ipH2p][ipH1s].Emis->Pdest) + 00132 Transitions[ipISO][nelem][ipH2p][ipH1s].Coll.ColUL * collider + 00133 Transitions[ipISO][nelem][ipH2p][ipH2s].Coll.ColUL * collider); 00134 } 00135 00136 00137 /* now renormalize to simple ionization ratio calculated in hydrolevel */ 00138 sum = 0.; 00139 /* 06 aug 28, from numLevels_max to _local. */ 00140 for( level=0; level < iso.numLevels_local[ipISO][nelem]; level++ ) 00141 { 00142 sum += StatesElem[ipISO][nelem][level].Pop; 00143 } 00144 00145 /* rescale so that inverse of sum is iso.xIonSimple[ipISO][nelem] */ 00146 renorm = 1. / SDIV( iso.xIonSimple[ipISO][nelem]) / SDIV( sum); 00147 00148 for( level=0; level < iso.numLevels_local[ipISO][nelem]; level++ ) 00149 { 00150 StatesElem[ipISO][nelem][level].Pop *= renorm; 00151 ASSERT( StatesElem[ipISO][nelem][level].Pop<BIGFLOAT ); 00152 if( iso.PopLTE[ipISO][nelem][level] > 1e-25 ) 00153 { 00154 iso.DepartCoef[ipISO][nelem][level] = 00155 StatesElem[ipISO][nelem][level].Pop/ 00156 (iso.PopLTE[ipISO][nelem][level]*dense.eden); 00157 } 00158 else 00159 { 00160 iso.DepartCoef[ipISO][nelem][level] = 0.; 00161 } 00162 } 00163 00164 if( trace.lgHBug && trace.lgTrace ) 00165 { 00166 fprintf( ioQQQ, " LOW TE,=%10.3e HN(1)=%10.3e rec=%10.3e Hgamnc(1s)=%10.3e\n", 00167 phycon.te, StatesElem[ipISO][nelem][ipH1s].Pop, ionbal.RateRecomTot[nelem][nelem-ipISO], 00168 iso.gamnc[ipISO][nelem][ipH1s] ); 00169 } 00170 00171 if( trace.lgTrace ) 00172 { 00173 fprintf( ioQQQ, 00174 " HydroT2Low return, LOW TE used, HII/HI=%.3e simple=%.3e IonRate=%.3e RecCo=%.3e\n", 00175 iso.pop_ion_ov_neut[ipISO][nelem], 00176 iso.xIonSimple[ipISO][nelem], 00177 iso.RateLevel2Cont[ipISO][nelem][ipH1s], 00178 ionbal.RateRecomTot[nelem][nelem-ipISO] ); 00179 } 00180 return; 00181 }