cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
iso_photo.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 /*iso_photo do photoionization rates for element nelem on the ipISO isoelectronic sequence */
4 #include "cddefines.h"
5 #include "hydrogenic.h"
6 #include "rfield.h"
7 #include "opacity.h"
8 #include "trace.h"
9 #include "ionbal.h"
10 #include "thermal.h"
11 #include "gammas.h"
12 #include "iso.h"
13 
14 void iso_photo(
15  /* iso sequence, hydrogen or helium for now */
16  long ipISO ,
17  /* the chemical element, 0 for hydrogen */
18  long int nelem)
19 {
20  long int limit ,
21  n;
22 
23  DEBUG_ENTRY( "iso_photo()" );
24 
25  /* check that we were called with valid charge */
26  ASSERT( nelem >= 0 && nelem < LIMELM );
27  ASSERT( ipISO < NISO );
28 
29  /* do photoionization rates */
30  /* induced recombination; FINDUC is integral of
31  * pho rate times EXP(-hn/kt) for induc rec
32  * CIND is this times hnu-hnu0 to get ind rec cooling
33  * ionbal.lgPhotoIoniz_On is 1, set to 0 with 'no photoionization' command
34  * ipSecIon points to 7.353 Ryd, lowest energy where secondary ioniz
35  * of hydrogen is possible */
36 
37  /* photoionization of ground, this is different from excited states because
38  * there will eventually be more than one shell, (when li-like done)
39  * and because upper limit is high-energy limit to code, so secondaries are
40  * included */
41  iso.gamnc[ipISO][nelem][0] = GammaBn(iso.ipIsoLevNIonCon[ipISO][nelem][0],
42  rfield.nflux,
43  iso.ipOpac[ipISO][nelem][0],
44  iso.xIsoLevNIonRyd[ipISO][nelem][0],
45  &iso.RecomInducRate[ipISO][nelem][0],
46  &iso.RecomInducCool_Coef[ipISO][nelem][0])*
48 
49  /* heating due to photo of ground */
51 
52  /* save these rates into ground photo rate vector */
53  ionbal.PhotoRate_Shell[nelem][nelem-ipISO][0][0] = iso.gamnc[ipISO][nelem][ipH1s];
54  ionbal.PhotoRate_Shell[nelem][nelem-ipISO][0][1] = thermal.HeatLowEnr*ionbal.lgPhotoIoniz_On;
55  ionbal.PhotoRate_Shell[nelem][nelem-ipISO][0][2] = thermal.HeatHiEnr*ionbal.lgPhotoIoniz_On;
56 
57  /* CompRecoilIonRate is direct photioniz rate due to
58  * bound compton scattering of very hard x-rays+Compton scat */
59  /* now add in compton recoil, to ground state, save heating as high energy heat */
60  ASSERT( ionbal.CompRecoilIonRate[nelem][nelem-ipISO]>=0. &&
61  ionbal.CompRecoilHeatRate[nelem][nelem-ipISO]>= 0. );
62  iso.gamnc[ipISO][nelem][0] += ionbal.CompRecoilIonRate[nelem][nelem-ipISO];
63  iso.PhotoHeat[ipISO][nelem][0] += ionbal.CompRecoilHeatRate[nelem][nelem-ipISO];
64 
65  /* now add bound compton scattering to ground term photoionization rate */
66  ionbal.PhotoRate_Shell[nelem][nelem-ipISO][0][0] += ionbal.CompRecoilIonRate[nelem][nelem-ipISO];
67  /* add heat to high energy heating term */
68  ionbal.PhotoRate_Shell[nelem][nelem-ipISO][0][2] += ionbal.CompRecoilHeatRate[nelem][nelem-ipISO];
69 
70  /* option to print ground state photoionization rates */
71  if( trace.lgTrace && trace.lgIsoTraceFull[ipISO] && (nelem == trace.ipIsoTrace[ipISO]) )
72  {
73  GammaPrt(iso.ipIsoLevNIonCon[ipISO][nelem][0],
74  rfield.nflux,
75  iso.ipOpac[ipISO][nelem][0],
76  ioQQQ,
77  iso.gamnc[ipISO][nelem][0],
78  iso.gamnc[ipISO][nelem][0]*0.05);
79  }
80 
81  /* for excited states upper limit to integration is ground threshold */
82  limit = iso.ipIsoLevNIonCon[ipISO][nelem][0]-1;
83  /* >>chng 06 aug 17, to numLevels_local instead of _max. */
84  for( n=1; n < iso.numLevels_local[ipISO][nelem]; n++ )
85  {
86  /* continuously update rates for n <=6, but only update
87  * rates for higher levels when redoing static opacities */
88  if( StatesElem[ipISO][nelem][n].n>6 && !opac.lgRedoStatic )
89  break;
94  if( hydro.lgHInducImp )
95  {
96  iso.gamnc[ipISO][nelem][n] =
97  GammaBn(
98  iso.ipIsoLevNIonCon[ipISO][nelem][n],
99  limit,
100  iso.ipOpac[ipISO][nelem][n],
101  iso.xIsoLevNIonRyd[ipISO][nelem][n],
102  &iso.RecomInducRate[ipISO][nelem][n],
103  &iso.RecomInducCool_Coef[ipISO][nelem][n])*
105  }
106  else
107  {
108  iso.gamnc[ipISO][nelem][n] =
109  GammaK(iso.ipIsoLevNIonCon[ipISO][nelem][n],
110  limit,
111  iso.ipOpac[ipISO][nelem][n],1.)*
113 
114  /* these are zero */
115  iso.RecomInducRate[ipISO][nelem][n] = 0.;
116  iso.RecomInducCool_Coef[ipISO][nelem][n] = 0.;
117  }
118  iso.PhotoHeat[ipISO][nelem][n] = thermal.HeatNet*ionbal.lgPhotoIoniz_On;
119 
120  ASSERT( iso.gamnc[ipISO][nelem][n]>= 0. );
121  ASSERT( iso.PhotoHeat[ipISO][nelem][n]>= 0. );
122  /* this loop only has excited states */
123  }
124 
125  {
126  /*@-redef@*/
127  enum {DEBUG_LOC=false};
128  /*@+redef@*/
129  if( DEBUG_LOC )
130  {
131  if( nelem==ipHYDROGEN )
132  {
133  fprintf(ioQQQ," buggbugg hphotodebugg%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
134  nzone,
135  iso.gamnc[ipISO][nelem][0],
136  iso.gamnc[ipISO][nelem][1],
137  iso.gamnc[ipISO][nelem][3],
138  iso.gamnc[ipISO][nelem][4],
139  iso.gamnc[ipISO][nelem][5],
140  iso.gamnc[ipISO][nelem][6]);
141  }
142  }
143  }
144 
145  /* >>chng 02 jan 19, kill excited state photoionization with case b no photo */
146  /* option for case b conditions, kill all excited state photoionizations */
147  if( opac.lgCaseB_no_photo )
148  {
149  for( n=1; n < iso.numLevels_max[ipISO][nelem]; n++ )
150  {
151  iso.gamnc[ipISO][nelem][n] = 0.;
152  iso.RecomInducRate[ipISO][nelem][n] = 0.;
153  iso.RecomInducCool_Coef[ipISO][nelem][n] = 0.;
154  }
155  }
156  {
157  /* this block turns off induced recom for some element */
158  /*@-redef@*/
159  enum {DEBUG_LOC=false};
160  /*@+redef@*/
161  if( DEBUG_LOC && ipISO==1 && nelem==5)
162  {
163  /* this debugging block is normally not active */
164  for( n=0; n < iso.numLevels_max[ipISO][nelem]; n++ )
165  {
166  iso.RecomInducRate[ipISO][nelem][n] = 0.;
167  }
168  }
169  }
170 
171  if( trace.lgTrace )
172  {
173  fprintf( ioQQQ, " iso_photo, ipISO%2ld nelem%2ld low, hi=",ipISO,nelem);
174  fprintf( ioQQQ,PrintEfmt("%9.2e", iso.gamnc[ipISO][nelem][ipH1s]));
175  ASSERT(nelem>=ipISO);
176  fprintf( ioQQQ,PrintEfmt("%9.2e", ionbal.CompRecoilIonRate[nelem][nelem-ipISO]));
177  fprintf( ioQQQ, " total=");
178  fprintf( ioQQQ,PrintEfmt("%9.2e",iso.gamnc[ipISO][nelem][ipH1s] ));
179  fprintf( ioQQQ, "\n");
180  }
181  return;
182 }

Generated for cloudy by doxygen 1.8.3.1