cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
highen.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*highen do high energy radiation field - gas interaction, Compton scattering, etc */
4 #include "cddefines.h"
5 #include "physconst.h"
6 #include "trace.h"
7 #include "heavy.h"
8 #include "radius.h"
9 #include "magnetic.h"
10 #include "hextra.h"
11 #include "thermal.h"
12 #include "dense.h"
13 #include "doppvel.h"
14 #include "ionbal.h"
15 #include "rfield.h"
16 #include "opacity.h"
17 #include "gammas.h"
18 #include "highen.h"
19 
20 void highen( void )
21 {
22  long int i,
23  ion,
24  nelem;
25 
26  double CosRayDen,
27  crsphi,
28  heatin,
29  sqthot;
30 
31  double hin;
32 
33  DEBUG_ENTRY( "highen()" );
34 
35 
36  /**********************************************************************
37  *
38  * COMPTON RECOIL IONIZATION
39  *
40  * bound electron scattering of >2.3 kev photons if neutral
41  * comoff usually 1, 0 if "NO COMPTON EFFECT" command given
42  * lgCompRecoil usually t, f if "NO RECOIL IONIZATION" command given
43  *
44  **********************************************************************/
45 
46  /* comoff turned off with no compton command,
47  * lgCompRecoil - no recoil ionization */
48  if( rfield.comoff > .0 && ionbal.lgCompRecoil )
49  {
50  for( nelem=0; nelem<LIMELM; ++nelem )
51  {
52  for( ion=0; ion<nelem+1; ++ion )
53  {
54  ionbal.CompRecoilIonRate[nelem][ion] = 0.;
55  ionbal.CompRecoilHeatRate[nelem][ion] = 0.;
56  if( dense.xIonDense[nelem][ion] > 0. )
57  {
58  /* recoil ionization starts at 194 Ryd = 2.6 keV */
59  /* this is the ionization potential of the valence shell */
60  /* >>chng 02 may 27, lower limit is now 1 beyond actual threshold
61  * since recoil energy at threshold was very small, sometimes negative */
62  for( i=ionbal.ipCompRecoil[nelem][ion]; i < rfield.nflux; ++i)
63  {
64  double recoil_energy;
65  crsphi = opac.OpacStack[i-1+opac.iopcom] * rfield.SummedCon[i] *
66  /* this is number of electrons in valence shell of this ion */
67  ionbal.nCompRecoilElec[nelem-ion];
68 
69  /* direct hydrogen ionization due to compton scattering,
70  * does not yet include secondaries,
71  * last term accounts for number of valence electrons that contribute */
72  ionbal.CompRecoilIonRate[nelem][ion] += crsphi;
73 
74  /* recoil energy in Rydbergs
75  * heating modified for suprathermal secondaries below; ANU2=ANU**2 */
76  /* >>chng 02 mar 27, use general formula for recoil energy */
77  /*energy = 2.66e-5*rfield.anu2[i] - 1.;*/
78  recoil_energy = rfield.anu2[i] /
80  Heavy.Valence_IP_Ryd[nelem][ion];
81 
82  /* heating is in rydbergs because SecIon2PrimaryErg, SecExcitLya2PrimaryErg, HeatEfficPrimary in ryd */
83  ionbal.CompRecoilHeatRate[nelem][ion] += crsphi*recoil_energy;
84  }
85  /* net heating rate, per atom, convert ryd/sec/cm3 to ergs/sec/atom */
86  ionbal.CompRecoilHeatRate[nelem][ion] *= EN1RYD;
87 
88  ASSERT( ionbal.CompRecoilHeatRate[nelem][ion] >= 0.);
89  ASSERT( ionbal.CompRecoilIonRate[nelem][ion] >= 0.);
90  }
91  }
92  }
93  }
94  else
95  {
96  for( nelem=0; nelem<LIMELM; ++nelem )
97  {
98  for( ion=0; ion<nelem+1; ++ion )
99  {
100  ionbal.CompRecoilIonRate[nelem][ion] = 0.;
101  ionbal.CompRecoilHeatRate[nelem][ion] = 0.;
102  }
103  }
104  }
105 
106  /**********************************************************************
107  *
108  * COSMIC RAYS
109  *
110  * heating and ionization by cosmic rays, other relativistic particles
111  * CRYDEN=density (1/CM**3), neutral rate assumes 15ev total
112  * energy loss, 13.6 into ionization, 1.4 into heating
113  * units erg/sec/cm**3
114  * iff not specified, CRTEMP is 2.6E9K
115  *
116  **********************************************************************/
117 
118  if( hextra.cryden > 0. )
119  {
120  ASSERT( hextra.crtemp > 0. );
121  /* this is current cosmic ray density, as determined from original density times
122  * possible dependence on radius */
124  {
125  /* >>chng 06 jun 02, add this option
126  * this is case where cr are in equipartition with magnetic field -
127  * set with COSMIC RAY EQUIPARTITION command */
128  CosRayDen = hextra.background_density *
129  /* ratio of energy density in current B to typical galactic
130  * galactic background energy density of 1.8 eV cm-3 is from
131  *>>refer cr background Webber, W.R. 1998, ApJ, 506, 329 */
133  (CR_EDEN_GAL_BACK_EV_CMM3/*1.8eV cm-3*/ * EN1EV/*erg eV-1*/ );
134  }
135  else
136  {
137  /* this is usual case, CR density may depend on radius, usually does not */
138  CosRayDen = hextra.cryden*pow(radius.Radius/radius.rinner,(double)hextra.crpowr);
139  }
140 
141  /* cosmic ray energy density rescaled by ratio to background ion rate and B field */
143  (CR_EDEN_GAL_BACK_EV_CMM3/*1.8eV cm-3*/ * EN1EV/*erg eV-1*/ );
144 
145  /* related to current temperature, when thermal */
146  sqthot = sqrt(hextra.crtemp);
147 
148  /* rate hot electrons heat cold ones, Balbus and McKee 1982
149  * units erg sec^-1 cm^-3,
150  * in sumheat we will multipy this rate by sum of neturals, but for this
151  * term we actually want eden, so mult by eden over sum of neut */
152  ionbal.CosRayHeatThermalElectrons = 5.5e-14/sqthot*CosRayDen;
153 
154  /* ionization rate; Balbus and McKee */
155  ionbal.CosRayIonRate = 1.22e-4/sqthot*
156  log(2.72*pow(hextra.crtemp/1e8,0.097))*CosRayDen;
157 
158  /* option for thermal CRs, first is the usual (and default) relativistic case */
159  if( hextra.crtemp > 2e9 )
160  {
161  /* usual circumstance; relativistic cosmic rays,
162  * cosmic ray ionization rate s-1 cm-3; ext rel limit */
163  ionbal.CosRayIonRate *= 3.;
164 
165  }
166  else
167  {
168  /* option for thermal cosmic rays */
169  ionbal.CosRayIonRate *= 10.;
170  }
171  /* >>chng 04 jan 27, from 0.093 to 2.574 as per following */
172  /* cr heating from Table 1 of
173  *>>refer cr heating Wolfire et al.1995, ApJ, 443, 152
174  * For every ionization due to cosmic rays, ~35eV of heat is added
175  * to the system. This manifests itself in the ionbal.CosRayHeatNeutralParticles term
176  * by the 2.574*EN1RYD term, which is just the energy in ergs in 35 eV.
177  * Change made by Nick Abel and Gargi Shaw, 04 Jan 27. In heatsum
178  * we multiply by the number of secondaries that occur */
180 
181  if( trace.lgTrace )
182  {
183  fprintf( ioQQQ, " highen: cosmic ray density;%10.2e CRion rate;%10.2e CR heat rate=;%10.2e CRtemp;%10.2e\n",
185  }
186  }
187  else
188  {
189  ionbal.CosRayIonRate = 0.;
191  }
192  /* >>chng 06 may 23, Penning ionization
193  ionbal.CosRayIonRate += 1e-9 *
194  StatesElem[ipHE_LIKE][ipHELIUM][ipHe2s3S].Pop*dense.xIonDense[ipHELIUM][1]; */
195 
196  /*fprintf(ioQQQ,"DEBUG cr %.2f %.3e %.3e %.3e\n",
197  fnzone,
198  hextra.cryden ,
199  ionbal.CosRayIonRate ,
200  ionbal.CosRayHeatNeutralParticles );*/
201 
202  /**********************************************************************
203  *
204  * add on extra heating due to turbulence, goes into [1] of [x][0][11][0]
205  *
206  **********************************************************************/
207 
208  /* TurbHeat added with hextra command, DispScale added with turbulence dissipative */
209  if( (hextra.TurbHeat+DoppVel.DispScale) != 0. )
210  {
211  /* turbulent heating only goes into the low-energy heat of this element */
212  /* >>>>chng 00 apr 28, functional form of radius dependence had bee turrad/depth
213  * and so went to infinity at the illuminated face. Changed to current form as
214  * per Ivan Hubeny comment */
215  if( hextra.lgHextraDepth )
216  {
217  /* if turrad is >0 then vary heating with depth */
220 
221  /* >>chng 00 aug 16, added option for heating from back side */
222  if( hextra.turback != 0. )
223  {
226  }
227  }
228  else if( hextra.lgHextraDensity )
229  {
230  /* depends on density */
233  }
234  /* this is turbulence dissipate command */
235  else if( DoppVel.DispScale > 0. )
236  {
237  double turb = DoppVel.TurbVel * sexp( radius.depth / DoppVel.DispScale );
238  /* if turrad is >0 then vary heating with depth */
239  /* >>chng 01 may 10, add extra factor of length over 1e13 cm */
240  ionbal.ExtraHeatRate = 3.45e-28 / 2.82824 * turb * turb * turb
241  / (DoppVel.DispScale/1e13);
242  }
243  else
244  {
245  /* constant extra heating */
247  }
248  }
249 
250  else
251  {
252  ionbal.ExtraHeatRate = 0.;
253  }
254 
255  /**********************************************************************
256  *
257  * option to add on fast neutron heating, goes into [0] & [2] of [x][0][11][0]
258  *
259  **********************************************************************/
260  if( hextra.lgNeutrnHeatOn )
261  {
262  /* hextra.totneu is energy flux erg cm-2 s-1
263  * CrsSecNeutron is 4E-26 cm^-2, cross sec for stopping rel neutrons
264  * this is heating erg s-1 due to fast neutrons, assumed to secondary ionize */
265  /* neutrons assumed to only secondary ionize */
267  }
268  else
269  {
271  }
272 
273 
274  /**********************************************************************
275  *
276  * pair production in elec field of nucleus
277  *
278  **********************************************************************/
282 
283  /**********************************************************************
284  *
285  * Compton energy exchange
286  *
287  **********************************************************************/
288  rfield.cmcool = 0.;
289  rfield.cmheat = 0.;
290  heatin = 0.;
291  /* comoff usually 1, turns off Compton */
292  if( rfield.comoff >0.01 )
293  {
294  for( i=0; i < rfield.nflux; i++ )
295  {
296 
297  /* Compton cooling
298  * CSIGC is Tarter expression times ANU(I)*3.858E-25
299  * 6.338E-6 is k in inf mass Rydbergs, still needs factor of TE */
300  rfield.comup[i] = (double)(rfield.flux[i]+rfield.ConInterOut[i]+
302  6.338e-6*1e-15);
303 
304  rfield.cmcool += rfield.comup[i];
305 
306  /* Compton heating
307  * CSIGH is Tarter expression times ANU(I)**2 * 3.858E-25
308  * CMHEAT is just spontaneous, HEATIN is just induced */
309  rfield.comdn[i] = (double)(rfield.flux[i]+rfield.ConInterOut[i]+
311 
312  /* induced Compton heating */
313  hin = (double)(rfield.flux[i]+rfield.ConInterOut[i]+rfield.outlin[i]+rfield.outlin_noplot[i])*
315  rfield.comdn[i] += hin;
316  heatin += hin;
317 
318  /* following is total compton heating */
319  rfield.cmheat += rfield.comdn[i];
320  }
321 
322  /* remember how important induced compton heating is */
323  if( rfield.cmheat > 0. )
324  {
326  }
327 
328  if( trace.lgTrace && trace.lgComBug )
329  {
330  fprintf( ioQQQ, " HIGHEN: COOL num=%8.2e HEAT num=%8.2e\n",
332  }
333  }
334 
335  /* save compton heating rate into main heating save array,
336  * no secondary ionizations from compton heating cooling */
337  thermal.heating[0][19] = rfield.cmheat;
338 
339  if( trace.lgTrace && trace.lgComBug )
340  {
341  fprintf( ioQQQ,
342  " HIGHEN finds heating fracs= frac(compt)=%10.2e "
343  " f(pair)%10.2e totHeat=%10.2e\n",
344  thermal.heating[0][19]/thermal.htot,
345  thermal.heating[0][21]/thermal.htot,
346  thermal.htot );
347  }
348  return;
349 }

Generated for cloudy by doxygen 1.8.4