cloudy  trunk
 All Data Structures Namespaces 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-2022 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 "iso.h"
12 #include "gammas.h"
13 #include "freebound.h"
14 #include "cosmology.h"
15 
16 void iso_photo(
17  /* iso sequence, hydrogen or helium for now */
18  long ipISO ,
19  /* the chemical element, 0 for hydrogen */
20  long int nelem)
21 {
22  long int limit ,
23  n;
24 
25  t_phoHeat photoHeat;
26 
27  DEBUG_ENTRY( "iso_photo()" );
28 
29  /* check that we were called with valid charge */
30  ASSERT( nelem >= 0 && nelem < LIMELM );
31  ASSERT( ipISO < NISO );
32 
33  t_iso_sp* sp = &iso_sp[ipISO][nelem];
34 
35  /* do photoionization rates */
36  /* induced recombination; FINDUC is integral of
37  * pho rate times EXP(-hn/kt) for induc rec
38  * CIND is this times hnu-hnu0 to get ind rec cooling
39  * ionbal.lgPhotoIoniz_On is 1, set to 0 with 'no photoionization' command
40  * ipSecIon points to 7.353 Ryd, lowest energy where secondary ioniz
41  * of hydrogen is possible */
42 
43  // photoionization of ground, this is different from excited states because
44  // there will eventually be more than one shell, (when li-like done)
45  sp->fb[0].gamnc = GammaBn(sp->fb[0].ipIsoLevNIonCon,
46  rfield.nflux,
47  sp->fb[0].ipOpac,
48  sp->fb[0].xIsoLevNIonRyd,
49  &sp->fb[0].RecomInducRate,
50  &sp->fb[0].RecomInducCool_Coef,
51  &photoHeat)*
53 
54  if( cosmology.lgDo )
55  {
56  photoHeat.HeatNet = 0.;
57  photoHeat.HeatLowEnr = 0.;
58  photoHeat.HeatHiEnr = 0.;
59  sp->fb[0].gamnc = 0.;
60  sp->fb[0].RecomInducRate = 0.;
61  sp->fb[0].RecomInducCool_Coef = 0.;
62  }
63 
64  /* heating due to photo of ground */
65  sp->fb[0].PhotoHeat = photoHeat.HeatNet*ionbal.lgPhotoIoniz_On;
66 
67  /* save these rates into ground photo rate vector */
68  ionbal.PhotoRate_Shell[nelem][nelem-ipISO][0][0] = sp->fb[ipH1s].gamnc;
69  ionbal.PhotoRate_Shell[nelem][nelem-ipISO][0][1] = photoHeat.HeatLowEnr*ionbal.lgPhotoIoniz_On;
70  ionbal.PhotoRate_Shell[nelem][nelem-ipISO][0][2] = photoHeat.HeatHiEnr*ionbal.lgPhotoIoniz_On;
71 
72  /* CompRecoilIonRate is direct photioniz rate due to
73  * bound compton scattering of very hard x-rays+Compton scat */
74  /* now add in compton recoil, to ground state, save heating as high energy heat */
75  ASSERT( ionbal.CompRecoilIonRate[nelem][nelem-ipISO]>=0. &&
76  ionbal.CompRecoilHeatRate[nelem][nelem-ipISO]>= 0. );
77  sp->fb[0].gamnc += ionbal.CompRecoilIonRate[nelem][nelem-ipISO];
78  sp->fb[0].PhotoHeat += ionbal.CompRecoilHeatRate[nelem][nelem-ipISO];
79 
80  /* now add bound compton scattering to ground term photoionization rate */
81  ionbal.PhotoRate_Shell[nelem][nelem-ipISO][0][0] += ionbal.CompRecoilIonRate[nelem][nelem-ipISO];
82  /* add heat to high energy heating term */
83  ionbal.PhotoRate_Shell[nelem][nelem-ipISO][0][2] += ionbal.CompRecoilHeatRate[nelem][nelem-ipISO];
84 
85  /* option to print ground state photoionization rates */
86  if( trace.lgTrace && trace.lgIsoTraceFull[ipISO] && (nelem == trace.ipIsoTrace[ipISO]) )
87  {
88  GammaPrt(sp->fb[0].ipIsoLevNIonCon,
89  rfield.nflux,
90  sp->fb[0].ipOpac,
91  ioQQQ,
92  sp->fb[0].gamnc,
93  sp->fb[0].gamnc*0.05);
94  }
95 
96  limit = rfield.nflux;
97  /* >>chng 06 aug 17, to numLevels_local instead of _max. */
98  for( n=1; n < sp->numLevels_local; n++ )
99  {
100  /* continuously update rates for n <=3, but only update
101  * rates for higher levels when redoing static opacities */
102  if( 0 && sp->st[n].n()>4 && !opac.lgRedoStatic && !sp->lgMustReeval )
103  break;
104 
109  if( hydro.lgHInducImp )
110  {
111  sp->fb[n].gamnc =
112  GammaBn(
113  sp->fb[n].ipIsoLevNIonCon,
114  limit,
115  sp->fb[n].ipOpac,
116  sp->fb[n].xIsoLevNIonRyd,
117  &sp->fb[n].RecomInducRate,
118  &sp->fb[n].RecomInducCool_Coef,
119  &photoHeat)*
121  }
122  else
123  {
124  sp->fb[n].gamnc =
125  GammaK(sp->fb[n].ipIsoLevNIonCon,
126  limit,
127  sp->fb[n].ipOpac,1.,
128  &photoHeat)*
130 
131  /* these are zero */
132  sp->fb[n].RecomInducRate = 0.;
133  sp->fb[n].RecomInducCool_Coef = 0.;
134  }
135  sp->fb[n].PhotoHeat = photoHeat.HeatNet*ionbal.lgPhotoIoniz_On;
136 
137  ASSERT( sp->fb[n].gamnc>= 0. );
138  ASSERT( sp->fb[n].PhotoHeat>= 0. );
139  /* this loop only has excited states */
140  }
141 
142  {
143  /*@-redef@*/
144  enum {DEBUG_LOC=false};
145  /*@+redef@*/
146  if( DEBUG_LOC )
147  {
148  if( nelem==ipHYDROGEN )
149  {
150  fprintf(ioQQQ," buggbugg hphotodebugg%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
151  nzone,
152  sp->fb[0].gamnc,
153  sp->fb[1].gamnc,
154  sp->fb[3].gamnc,
155  sp->fb[4].gamnc,
156  sp->fb[5].gamnc,
157  sp->fb[6].gamnc);
158  }
159  }
160  }
161 
162  /* >>chng 02 jan 19, kill excited state photoionization with case b no photo */
163  /* option for case b conditions, kill all excited state photoionizations */
164  if( opac.lgCaseB_no_photo )
165  {
166  for( n=1; n < sp->numLevels_max; n++ )
167  {
168  sp->fb[n].gamnc = 0.;
169  sp->fb[n].PhotoHeat = 0.;
170  sp->fb[n].RecomInducRate = 0.;
171  sp->fb[n].RecomInducCool_Coef = 0.;
172  }
173  }
174  {
175  /* this block turns off induced recom for some element */
176  /*@-redef@*/
177  enum {DEBUG_LOC=false};
178  /*@+redef@*/
179  if( DEBUG_LOC && ipISO==1 && nelem==5)
180  {
181  /* this debugging block is normally not active */
182  for( n=0; n < sp->numLevels_max; n++ )
183  {
184  sp->fb[n].RecomInducRate = 0.;
185  }
186  }
187  }
188 
189  if( trace.lgTrace && (trace.lgHBug||trace.lgHeBug) )
190  {
191  fprintf( ioQQQ, " iso_photo, ipISO%2ld nelem%2ld low, hi=",ipISO,nelem);
192  fprintf( ioQQQ,PrintEfmt("%9.2e", sp->fb[ipH1s].gamnc));
193  ASSERT(nelem>=ipISO);
194  fprintf( ioQQQ,PrintEfmt("%9.2e", ionbal.CompRecoilIonRate[nelem][nelem-ipISO]));
195  fprintf( ioQQQ, " total=");
196  fprintf( ioQQQ,PrintEfmt("%9.2e",sp->fb[ipH1s].gamnc ));
197  fprintf( ioQQQ, "\n");
198  }
199  return;
200 }
double HeatHiEnr
Definition: thermal.h:208
bool lgCaseB_no_photo
Definition: opacity.h:182
qList st
Definition: iso.h:482
double GammaK(long int ipLoEnr, long int ipHiEnr, long int ipOpac, double yield1, t_phoHeat *photoHeat)
double ** CompRecoilIonRate
Definition: ionbal.h:162
bool lgHeBug
Definition: trace.h:79
t_opac opac
Definition: opacity.cpp:5
double ** CompRecoilHeatRate
Definition: ionbal.h:168
const int NISO
Definition: cddefines.h:311
bool lgIsoTraceFull[NISO]
Definition: trace.h:85
#define PrintEfmt(F, V)
Definition: cddefines.h:1349
bool lgPhotoIoniz_On
Definition: ionbal.h:117
long int ipIsoTrace[NISO]
Definition: trace.h:88
bool lgHInducImp
Definition: hydrogenic.h:132
FILE * ioQQQ
Definition: cddefines.cpp:7
long int nzone
Definition: cddefines.cpp:14
void iso_photo(long ipISO, long nelem)
vector< freeBound > fb
Definition: iso.h:481
bool lgDo
Definition: cosmology.h:44
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
t_trace trace
Definition: trace.cpp:5
t_ionbal ionbal
Definition: ionbal.cpp:8
const int ipH1s
Definition: iso.h:29
bool lgTrace
Definition: trace.h:12
t_rfield rfield
Definition: rfield.cpp:9
t_hydro hydro
Definition: hydrogenic.cpp:5
double HeatLowEnr
Definition: thermal.h:206
double **** PhotoRate_Shell
Definition: ionbal.h:112
Definition: iso.h:470
#define ASSERT(exp)
Definition: cddefines.h:613
bool lgHBug
Definition: trace.h:82
const int LIMELM
Definition: cddefines.h:308
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
t_cosmology cosmology
Definition: cosmology.cpp:8
double GammaBn(long int ipLoEnr, long int ipHiEnr, long int ipOpac, double thresh, double *ainduc, double *rcool, t_phoHeat *photoHeat)
Definition: cont_gammas.cpp:34
bool lgRedoStatic
Definition: opacity.h:160
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
double HeatNet
Definition: thermal.h:204
long int numLevels_max
Definition: iso.h:524
void GammaPrt(long int ipLoEnr, long int ipHiEnr, long int ipOpac, FILE *ioFILE, double total, double threshold)
const int ipHYDROGEN
Definition: cddefines.h:349
long int nflux
Definition: rfield.h:46
long int numLevels_local
Definition: iso.h:529
bool lgMustReeval
Definition: iso.h:512