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