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 /*highen do high energy radiation field - gas interaction, Compton scattering, etc */ 00004 #include "cddefines.h" 00005 #include "physconst.h" 00006 #include "trace.h" 00007 #include "heavy.h" 00008 #include "radius.h" 00009 #include "magnetic.h" 00010 #include "hextra.h" 00011 #include "thermal.h" 00012 #include "dense.h" 00013 #include "doppvel.h" 00014 #include "ionbal.h" 00015 #include "rfield.h" 00016 #include "opacity.h" 00017 #include "gammas.h" 00018 #include "highen.h" 00019 00020 void highen( void ) 00021 { 00022 long int i, 00023 ion, 00024 nelem; 00025 00026 double CosRayDen, 00027 crsphi, 00028 heatin, 00029 sqthot; 00030 00031 double hin; 00032 00033 DEBUG_ENTRY( "highen()" ); 00034 00035 00036 /********************************************************************** 00037 * 00038 * COMPTON RECOIL IONIZATION 00039 * 00040 * bound electron scattering of >2.3 kev photons if neutral 00041 * comoff usually 1, 0 if "NO COMPTON EFFECT" command given 00042 * lgCompRecoil usually t, f if "NO RECOIL IONIZATION" command given 00043 * 00044 **********************************************************************/ 00045 00046 /* comoff turned off with no compton command, 00047 * lgCompRecoil - no recoil ionization */ 00048 if( rfield.comoff > .0 && ionbal.lgCompRecoil ) 00049 { 00050 for( nelem=0; nelem<LIMELM; ++nelem ) 00051 { 00052 for( ion=0; ion<nelem+1; ++ion ) 00053 { 00054 ionbal.CompRecoilIonRate[nelem][ion] = 0.; 00055 ionbal.CompRecoilHeatRate[nelem][ion] = 0.; 00056 if( dense.xIonDense[nelem][ion] > 0. ) 00057 { 00058 /* recoil ionization starts at 194 Ryd = 2.6 keV */ 00059 /* this is the ionization potential of the valence shell */ 00060 /* >>chng 02 may 27, lower limit is now 1 beyond actual threshold 00061 * since recoil energy at threshold was very small, sometimes negative */ 00062 for( i=ionbal.ipCompRecoil[nelem][ion]; i < rfield.nflux; ++i) 00063 { 00064 double recoil_energy; 00065 crsphi = opac.OpacStack[i-1+opac.iopcom] * rfield.SummedCon[i] * 00066 /* this is number of electrons in valence shell of this ion */ 00067 ionbal.nCompRecoilElec[nelem-ion]; 00068 00069 /* direct hydrogen ionization due to compton scattering, 00070 * does not yet include secondaries, 00071 * last term accounts for number of valence electrons that contribute */ 00072 ionbal.CompRecoilIonRate[nelem][ion] += crsphi; 00073 00074 /* recoil energy in Rydbergs 00075 * heating modified for suprathermal secondaries below; ANU2=ANU**2 */ 00076 /* >>chng 02 mar 27, use general formula for recoil energy */ 00077 /*energy = 2.66e-5*rfield.anu2[i] - 1.;*/ 00078 recoil_energy = rfield.anu2[i] / 00079 ( EN1RYD * ELECTRON_MASS * SPEEDLIGHT * SPEEDLIGHT) * EN1RYD* EN1RYD - 00080 Heavy.Valence_IP_Ryd[nelem][ion]; 00081 00082 /* heating is in rydbergs because SecIon2PrimaryErg, SecExcitLya2PrimaryErg, HeatEfficPrimary in ryd */ 00083 ionbal.CompRecoilHeatRate[nelem][ion] += crsphi*recoil_energy; 00084 } 00085 /* net heating rate, per atom, convert ryd/sec/cm3 to ergs/sec/atom */ 00086 ionbal.CompRecoilHeatRate[nelem][ion] *= EN1RYD; 00087 00088 ASSERT( ionbal.CompRecoilHeatRate[nelem][ion] >= 0.); 00089 ASSERT( ionbal.CompRecoilIonRate[nelem][ion] >= 0.); 00090 } 00091 } 00092 } 00093 } 00094 else 00095 { 00096 for( nelem=0; nelem<LIMELM; ++nelem ) 00097 { 00098 for( ion=0; ion<nelem+1; ++ion ) 00099 { 00100 ionbal.CompRecoilIonRate[nelem][ion] = 0.; 00101 ionbal.CompRecoilHeatRate[nelem][ion] = 0.; 00102 } 00103 } 00104 } 00105 00106 /********************************************************************** 00107 * 00108 * COSMIC RAYS 00109 * 00110 * heating and ionization by cosmic rays, other relativistic particles 00111 * CRYDEN=density (1/CM**3), neutral rate assumes 15ev total 00112 * energy loss, 13.6 into ionization, 1.4 into heating 00113 * units erg/sec/cm**3 00114 * iff not specified, CRTEMP is 2.6E9K 00115 * 00116 **********************************************************************/ 00117 00118 if( hextra.cryden > 0. ) 00119 { 00120 ASSERT( hextra.crtemp > 0. ); 00121 /* this is current cosmic ray density, as determined from original density times 00122 * possible dependence on radius */ 00123 if( hextra.lg_CR_B_equipartition ) 00124 { 00125 /* >>chng 06 jun 02, add this option 00126 * this is case where cr are in equipartition with magnetic field - 00127 * set with COSMIC RAY EQUIPARTITION command */ 00128 CosRayDen = hextra.background_density * 00129 /* ratio of energy density in current B to typical galactic 00130 * galactic background energy density of 1.8 eV cm-3 is from 00131 *>>refer cr background Webber, W.R. 1998, ApJ, 506, 329 */ 00132 magnetic.energydensity / 00133 (CR_EDEN_GAL_BACK_EV_CMM3/*1.8eV cm-3*/ * EN1EV/*erg eV-1*/ ); 00134 } 00135 else 00136 { 00137 /* this is usual case, CR density may depend on radius, usually does not */ 00138 CosRayDen = hextra.cryden*pow(radius.Radius/radius.rinner,(double)hextra.crpowr); 00139 } 00140 00141 /* cosmic ray energy density rescaled by ratio to background ion rate and B field */ 00142 hextra.cr_energydensity = CosRayDen/hextra.background_density * 00143 (CR_EDEN_GAL_BACK_EV_CMM3/*1.8eV cm-3*/ * EN1EV/*erg eV-1*/ ); 00144 00145 /* related to current temperature, when thermal */ 00146 sqthot = sqrt(hextra.crtemp); 00147 00148 /* rate hot electrons heat cold ones, Balbus and McKee 1982 00149 * units erg sec^-1 cm^-3, 00150 * in sumheat we will multipy this rate by sum of neturals, but for this 00151 * term we actually want eden, so mult by eden over sum of neut */ 00152 ionbal.CosRayHeatThermalElectrons = 5.5e-14/sqthot*CosRayDen; 00153 00154 /* ionization rate; Balbus and McKee */ 00155 ionbal.CosRayIonRate = 1.22e-4/sqthot* 00156 log(2.72*pow(hextra.crtemp/1e8,0.097))*CosRayDen; 00157 00158 /* option for thermal CRs, first is the usual (and default) relativistic case */ 00159 if( hextra.crtemp > 2e9 ) 00160 { 00161 /* usual circumstance; relativistic cosmic rays, 00162 * cosmic ray ionization rate s-1 cm-3; ext rel limit */ 00163 ionbal.CosRayIonRate *= 3.; 00164 00165 } 00166 else 00167 { 00168 /* option for thermal cosmic rays */ 00169 ionbal.CosRayIonRate *= 10.; 00170 } 00171 /* >>chng 04 jan 27, from 0.093 to 2.574 as per following */ 00172 /* cr heating from Table 1 of 00173 *>>refer cr heating Wolfire et al.1995, ApJ, 443, 152 00174 * For every ionization due to cosmic rays, ~35eV of heat is added 00175 * to the system. This manifests itself in the ionbal.CosRayHeatNeutralParticles term 00176 * by the 2.574*EN1RYD term, which is just the energy in ergs in 35 eV. 00177 * Change made by Nick Abel and Gargi Shaw, 04 Jan 27. In heatsum 00178 * we multiply by the number of secondaries that occur */ 00179 ionbal.CosRayHeatNeutralParticles = ionbal.CosRayIonRate*2.574*EN1RYD; 00180 00181 if( trace.lgTrace ) 00182 { 00183 fprintf( ioQQQ, " highen: cosmic ray density;%10.2e CRion rate;%10.2e CR heat rate=;%10.2e CRtemp;%10.2e\n", 00184 CosRayDen, ionbal.CosRayIonRate, ionbal.CosRayHeatNeutralParticles, hextra.crtemp ); 00185 } 00186 } 00187 else 00188 { 00189 ionbal.CosRayIonRate = 0.; 00190 ionbal.CosRayHeatNeutralParticles = 0.; 00191 } 00192 /* >>chng 06 may 23, Penning ionization 00193 ionbal.CosRayIonRate += 1e-9 * 00194 StatesElem[ipHE_LIKE][ipHELIUM][ipHe2s3S].Pop*dense.xIonDense[ipHELIUM][1]; */ 00195 00196 /*fprintf(ioQQQ,"DEBUG cr %.2f %.3e %.3e %.3e\n", 00197 fnzone, 00198 hextra.cryden , 00199 ionbal.CosRayIonRate , 00200 ionbal.CosRayHeatNeutralParticles );*/ 00201 00202 /********************************************************************** 00203 * 00204 * add on extra heating due to turbulence, goes into [1] of [x][0][11][0] 00205 * 00206 **********************************************************************/ 00207 00208 /* TurbHeat added with hextra command, DispScale added with turbulence dissipative */ 00209 if( (hextra.TurbHeat+DoppVel.DispScale) != 0. ) 00210 { 00211 /* turbulent heating only goes into the low-energy heat of this element */ 00212 /* >>>>chng 00 apr 28, functional form of radius dependence had bee turrad/depth 00213 * and so went to infinity at the illuminated face. Changed to current form as 00214 * per Ivan Hubeny comment */ 00215 if( hextra.lgHextraDepth ) 00216 { 00217 /* if turrad is >0 then vary heating with depth */ 00218 ionbal.ExtraHeatRate = 00219 hextra.TurbHeat*sexp(radius.depth /hextra.turrad); 00220 00221 /* >>chng 00 aug 16, added option for heating from back side */ 00222 if( hextra.turback != 0. ) 00223 { 00224 ionbal.ExtraHeatRate += 00225 hextra.TurbHeat*sexp((hextra.turback - radius.depth) /hextra.turrad); 00226 } 00227 } 00228 else if( hextra.lgHextraDensity ) 00229 { 00230 /* depends on density */ 00231 ionbal.ExtraHeatRate = 00232 hextra.TurbHeat*dense.gas_phase[ipHYDROGEN]/hextra.HextraScaleDensity; 00233 } 00234 /* this is turbulence dissipate command */ 00235 else if( DoppVel.DispScale > 0. ) 00236 { 00237 double turb = DoppVel.TurbVel * sexp( radius.depth / DoppVel.DispScale ); 00238 /* if turrad is >0 then vary heating with depth */ 00239 /* >>chng 01 may 10, add extra factor of length over 1e13 cm */ 00240 ionbal.ExtraHeatRate = 3.45e-28 / 2.82824 * turb * turb * turb 00241 / (DoppVel.DispScale/1e13); 00242 } 00243 else 00244 { 00245 /* constant extra heating */ 00246 ionbal.ExtraHeatRate = hextra.TurbHeat; 00247 } 00248 } 00249 00250 else 00251 { 00252 ionbal.ExtraHeatRate = 0.; 00253 } 00254 00255 /********************************************************************** 00256 * 00257 * option to add on fast neutron heating, goes into [0] & [2] of [x][0][11][0] 00258 * 00259 **********************************************************************/ 00260 if( hextra.lgNeutrnHeatOn ) 00261 { 00262 /* hextra.totneu is energy flux erg cm-2 s-1 00263 * CrsSecNeutron is 4E-26 cm^-2, cross sec for stopping rel neutrons 00264 * this is heating erg s-1 due to fast neutrons, assumed to secondary ionize */ 00265 /* neutrons assumed to only secondary ionize */ 00266 ionbal.xNeutronHeatRate = hextra.totneu*hextra.CrsSecNeutron; 00267 } 00268 else 00269 { 00270 ionbal.xNeutronHeatRate = 0.; 00271 } 00272 00273 00274 /********************************************************************** 00275 * 00276 * pair production in elec field of nucleus 00277 * 00278 **********************************************************************/ 00279 ionbal.PairProducPhotoRate[0] = GammaK(opac.ippr,rfield.nflux,opac.ioppr,1.); 00280 ionbal.PairProducPhotoRate[1] = thermal.HeatLowEnr; 00281 ionbal.PairProducPhotoRate[2] = thermal.HeatHiEnr; 00282 00283 /********************************************************************** 00284 * 00285 * Compton energy exchange 00286 * 00287 **********************************************************************/ 00288 rfield.cmcool = 0.; 00289 rfield.cmheat = 0.; 00290 heatin = 0.; 00291 /* comoff usually 1, turns off Compton */ 00292 if( rfield.comoff >0.01 ) 00293 { 00294 for( i=0; i < rfield.nflux; i++ ) 00295 { 00296 00297 /* Compton cooling 00298 * CSIGC is Tarter expression times ANU(I)*3.858E-25 00299 * 6.338E-6 is k in inf mass Rydbergs, still needs factor of TE */ 00300 rfield.comup[i] = (double)(rfield.flux[i]+rfield.ConInterOut[i]+ 00301 rfield.outlin[i]+ rfield.outlin_noplot[i])*rfield.csigc[i]*(dense.eden*4.e0* 00302 6.338e-6*1e-15); 00303 00304 rfield.cmcool += rfield.comup[i]; 00305 00306 /* Compton heating 00307 * CSIGH is Tarter expression times ANU(I)**2 * 3.858E-25 00308 * CMHEAT is just spontaneous, HEATIN is just induced */ 00309 rfield.comdn[i] = (double)(rfield.flux[i]+rfield.ConInterOut[i]+ 00310 rfield.outlin[i]+ rfield.outlin_noplot[i])*rfield.csigh[i]*dense.eden*1e-15; 00311 00312 /* induced Compton heating */ 00313 hin = (double)(rfield.flux[i]+rfield.ConInterOut[i]+rfield.outlin[i]+rfield.outlin_noplot[i])* 00314 rfield.csigh[i]*rfield.OccNumbIncidCont[i]*dense.eden*1e-15; 00315 rfield.comdn[i] += hin; 00316 heatin += hin; 00317 00318 /* following is total compton heating */ 00319 rfield.cmheat += rfield.comdn[i]; 00320 } 00321 00322 /* remember how important induced compton heating is */ 00323 if( rfield.cmheat > 0. ) 00324 { 00325 rfield.cinrat = MAX2(rfield.cinrat,heatin/rfield.cmheat); 00326 } 00327 00328 if( trace.lgTrace && trace.lgComBug ) 00329 { 00330 fprintf( ioQQQ, " HIGHEN: COOL num=%8.2e HEAT num=%8.2e\n", 00331 rfield.cmcool, rfield.cmheat ); 00332 } 00333 } 00334 00335 /* save compton heating rate into main heating save array, 00336 * no secondary ionizations from compton heating cooling */ 00337 thermal.heating[0][19] = rfield.cmheat; 00338 00339 if( trace.lgTrace && trace.lgComBug ) 00340 { 00341 fprintf( ioQQQ, 00342 " HIGHEN finds heating fracs= frac(compt)=%10.2e " 00343 " f(pair)%10.2e totHeat=%10.2e\n", 00344 thermal.heating[0][19]/thermal.htot, 00345 thermal.heating[0][21]/thermal.htot, 00346 thermal.htot ); 00347 } 00348 return; 00349 }