cloudy trunk
Loading...
Searching...
No Matches
ion_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/*ion_photo fill array PhotoRate with photoionization rates for heavy elements */
4#include "cddefines.h"
5#include "yield.h"
6#include "heavy.h"
7#include "opacity.h"
8#include "dense.h"
9#include "thermal.h"
10#include "conv.h"
11#include "grainvar.h"
12#include "elementnames.h"
13#include "gammas.h"
14#include "ionbal.h"
15#include "ca.h"
16#include "mole.h"
17#include "phycon.h"
18#include "hmi.h"
19#include "rfield.h"
20#include "atoms.h"
21#include "iso.h"
22#include "oxy.h"
23#include "atmdat.h"
24#include "fe.h"
25
27 /* nlem is atomic number on C scale, 0 for H */
28 long int nelem ,
29 /* debugging flag to turn on print */
30 bool lgPrintIt )
31{
32 long int ion,
33 iphi,
34 iplow,
35 ipop,
36 limit_hi,
37 limit_lo,
38 ns;
39
40 DEBUG_ENTRY( "ion_photo()" );
41
42 /* IonLow(nelem) and IonHigh(nelem) are bounds for evaluation*/
43
44 /*begin sanity checks */
45 ASSERT( nelem < LIMELM );
46 ASSERT( dense.IonLow[nelem] >= 0 );
47 ASSERT( dense.IonLow[nelem] <= dense.IonHigh[nelem] );
48 ASSERT( dense.IonHigh[nelem] <= nelem + 1);
49 /*end sanity checks */
50
51 /* NB - in following, nelem is on c scale, so is 29 for Zn */
52
53 /* min since iso-like atom rates produced in iso_photo.
54 * IonHigh and IonLow range from 0 to nelem+1, but for photo rates
55 * we want 0 to nelem since cannot destroy fully stripped ion.
56 * since iso-seq done elsewhere, want to actually do IonHigh-xxx.*/
57 /* >>chng 00 dec 07, logic on limit_hi now precisely identical to ion_solver */
58 /* >>chng 02 mar 31, change limit_hi to < in loop (had been <=) and
59 * also coded for arbitrary number of iso sequences */
60 /*limit_hi = MIN2(nelem,dense.IonHigh[nelem]);
61 limit_hi = MIN2(nelem-2,dense.IonHigh[nelem]-1);*/
62 limit_hi = MIN2( dense.IonHigh[nelem] , nelem+1-NISO );
63
64 /* >>chng 03 sep 26, always do atom itself since may be needed for molecules */
65 limit_hi = MAX2( 1 , limit_hi );
66
67 /* when grains are present want to use atoms as lower bound to number of stages of ionization,
68 * since atomic rates needed for species within grains */
69 if( !conv.nPres2Ioniz && gv.lgDustOn() )
70 {
71 limit_lo = 0;
72 }
73 else
74 {
75 limit_lo = dense.IonLow[nelem];
76 }
77
78 /* >>chng 01 dec 11, lower bound now limit_lo */
79 /* loop over all ions for this element */
80 /* >>chng 02 mar 31, now ion < limit_hi not <= */
81 for( ion=limit_lo; ion < limit_hi; ion++ )
82 {
83 /* loop over all shells for this ion */
84 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
85 {
86 /* always reevaluate the outer shell, and all shells if lgRedoStatic is set */
87 if( (ns==(Heavy.nsShells[nelem][ion]-1) || opac.lgRedoStatic) )
88 {
89 /* option to redo the rates only on occasion */
90 iplow = opac.ipElement[nelem][ion][ns][0];
91 iphi = opac.ipElement[nelem][ion][ns][1];
92 ipop = opac.ipElement[nelem][ion][ns][2];
93
94 t_phoHeat photoHeat;
95
96 /* compute the photoionization rate, ionbal.lgPhotoIoniz_On is 1, set 0
97 * with "no photoionization" command */
98 ionbal.PhotoRate_Shell[nelem][ion][ns][0] =
99 GammaK(iplow,iphi,
100 ipop,t_yield::Inst().elec_eject_frac(nelem,ion,ns,0),
101 &photoHeat )*ionbal.lgPhotoIoniz_On;
102
103 /* these three lines must be kept parallel with the lines
104 * in GammaK ion*/
105
106 /* the heating rate */
107 ionbal.PhotoRate_Shell[nelem][ion][ns][1] = photoHeat.HeatLowEnr*ionbal.lgPhotoIoniz_On;
108 ionbal.PhotoRate_Shell[nelem][ion][ns][2] = photoHeat.HeatHiEnr*ionbal.lgPhotoIoniz_On;
109 }
110 }
111
112 /* add on compton recoil ionization for atoms to outer shell */
113 /* >>chng 02 mar 24, moved here from ion_solver */
114 /* this is the outer shell */
115 ns = (Heavy.nsShells[nelem][ion]-1);
116 /* this must be moved to photoionize and have code parallel to iso_photo code */
117 ionbal.PhotoRate_Shell[nelem][ion][ns][0] += ionbal.CompRecoilIonRate[nelem][ion];
118 /* add the heat as secondary-ionization capable heating */
119 ionbal.PhotoRate_Shell[nelem][ion][ns][2] += ionbal.CompRecoilHeatRate[nelem][ion];
120 }
121
122 /* option to print information about these rates for this element */
123 if( lgPrintIt )
124 {
125 /* option to print rates for particular shell */
126 ns = 5;
127 ion = 1;
128 GammaPrt(
129 opac.ipElement[nelem][ion][ns][0],
130 opac.ipElement[nelem][ion][ns][1],
131 opac.ipElement[nelem][ion][ns][2],
132 ioQQQ, /* io unit we will write to */
133 ionbal.PhotoRate_Shell[nelem][ion][ns][0],
134 0.05);
135
136 /* outer loop is from K to most number of shells present in atom */
137 for( ns=0; ns < Heavy.nsShells[nelem][0]; ns++ )
138 {
139 fprintf( ioQQQ, "\n %s", elementnames.chElementNameShort[nelem] );
140 fprintf( ioQQQ, " %s" , Heavy.chShell[ns]);
141 /* MB hydrogenic photo rate may not be included in beow */
142 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
143 {
144 if( Heavy.nsShells[nelem][ion] > ns )
145 {
146 fprintf( ioQQQ, " %8.1e", ionbal.PhotoRate_Shell[nelem][ion][ns][0] );
147 }
148 else
149 {
150 break;
151 }
152 }
153 }
154 fprintf(ioQQQ,"\n");
155 }
156 /* >>chng 11 may 06, moved from ion_calci.cpp */
157 /* Ly-alpha photoionization of Ca+
158 * valence shell is reevaluated by ion_photo on every call, so this does not double count */
159 if( nelem == ipCALCIUM )
160 {
161 long ns = 6, ion = 1;
162 ionbal.PhotoRate_Shell[nelem][ion][ns][0] += ca.dstCala;
163 }
164 if( nelem == ipCARBON )
165 {
166 /* >>chng 05 aug 10, add Leiden hack here to get same C0 photo rate
167 * as UMIST - negates difference in grain opacities */
168 if(mole_global.lgLeidenHack)
169 {
170 int nelem=ipCARBON , ion=0 , ns=2;
171 ionbal.PhotoRate_Shell[nelem][ion][ns][0] =
172 (HMRATE((1e-10)*3.0,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*
173 exp(-(3.0*rfield.extin_mag_V_point))/1.66));
174 /* heating rates */
175 ionbal.PhotoRate_Shell[nelem][ion][ns][1] = 0.;
176 ionbal.PhotoRate_Shell[nelem][ion][ns][2] = 0.;
177 }
178 }
179 if( nelem == ipNITROGEN )
180 {
181 /* photoionization from 2D of NI, is atomic nitrogen present? */
182 if( dense.xIonDense[ipNITROGEN][0] > 0. )
183 {
184 t_phoHeat photoHeat;
185 // photo rate, population atoms.p2nit evaluated in cooling
186 atoms.d5200r = (realnum)GammaK(opac.in1[0],opac.in1[1],opac.in1[2],1.,&photoHeat);
187 /* valence shell photoionization, followed by heating; [0][6] => atomic nitrogen */
188 ionbal.PhotoRate_Shell[ipNITROGEN][0][2][0] = ionbal.PhotoRate_Shell[ipNITROGEN][0][2][0]*
189 (1. - atoms.p2nit) + atoms.p2nit*atoms.d5200r;
190 ionbal.PhotoRate_Shell[ipNITROGEN][0][2][1] = ionbal.PhotoRate_Shell[ipNITROGEN][0][2][1]*
191 (1. - atoms.p2nit) + photoHeat.HeatNet*atoms.p2nit;
192 }
193 else
194 {
195 atoms.p2nit = 0.;
196 atoms.d5200r = 0.;
197 }
198 }
199 if ( nelem == ipMAGNESIUM )
200 {
201 if( dense.IonLow[ipMAGNESIUM] <= 1 )
202 {
203 t_phoHeat dummy;
204 /* photoionization from excited upper state of 2798 */
205 realnum rmg2l = (realnum)GammaK(opac.ipmgex,
206 iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon,opac.ipOpMgEx,1., &dummy );
207 ionbal.PhotoRate_Shell[ipMAGNESIUM][1][3][0] += rmg2l*atoms.popmg2;
208
209 if( nzone <= 1 )
210 {
211 atoms.xMg2Max = 0.;
212 }
213 else if( ionbal.PhotoRate_Shell[ipMAGNESIUM][1][3][0] > 1e-30 )
214 {
215 /* remember max relative photoionization rate for possible comment */
216 atoms.xMg2Max = (realnum)(MAX2(atoms.xMg2Max,rmg2l*atoms.popmg2/
217 ionbal.PhotoRate_Shell[ipMAGNESIUM][1][3][0]));
218 }
219 }
220 else
221 {
222 atoms.xMg2Max = 0.;
223 }
224 }
225
226 if ( nelem == ipOXYGEN )
227 {
228
229 t_phoHeat dummy;
230 /* oxygen, atomic number 8 */
231 if( !dense.lgElmtOn[ipOXYGEN] )
232 {
233 oxy.poiii2Max = 0.;
234 oxy.poiii3Max = 0.;
235 oxy.r4363Max = 0.;
236 oxy.r5007Max = 0.;
237 oxy.poiii2 = 0.;
238 oxy.AugerO3 = 0.;
239 oxy.s3727 = 0.;
240 oxy.s7325 = 0.;
241 thermal.heating[7][9] = 0.;
242 oxy.poimax = 0.;
243 return;
244 }
245 else
246 {
247 double aeff;
248
249 /* photoionization from O++ 1D
250 *
251 * estimate gamma function by assuming no frequency dependence
252 * betwen 1D and O++3P edge */
253 /* destroy upper level of OIII 5007*/
254 oxy.d5007r = (realnum)(GammaK(opac.ipo3exc[0],opac.ipo3exc[1],
255 opac.ipo3exc[2] , 1., &dummy ));
256
257 /* destroy upper level of OIII 4363*/
258 oxy.d4363 = (realnum)(GammaK(opac.ipo3exc3[0],opac.ipo3exc3[1],
259 opac.ipo3exc3[2] , 1., &dummy ));
260
261 /* destroy upper level of OI 6300*/
262 oxy.d6300 = (realnum)(GammaK(opac.ipo1exc[0],opac.ipo1exc[1],
263 opac.ipo1exc[2] , 1., &dummy ));
264
265 /* A21 = 0.0263 */
266 aeff = 0.0263 + oxy.d5007r;
267
268 /* 1. as last arg makes this the relative population */
269 oxy.poiii2 = (realnum)(atom_pop2(2.5,9.,5.,aeff,2.88e4,1.)/aeff);
270 {
271 enum {DEBUG_LOC=false};
272 if( DEBUG_LOC )
273 {
274 fprintf(ioQQQ,"pop rel %.1e rate %.1e grnd rate %.1e\n",
275 oxy.poiii2 , oxy.d5007r ,ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] );
276 }
277 }
278
279 /* photoionization from excited states */
280 if( nzone > 0 )
281 {
282 /* neutral oxygen destruction */
283 ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0] = ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0]*
284 (1. - oxy.poiexc) + oxy.d6300*oxy.poiexc;
285
286 /* doubly ionized oxygen destruction */
287 ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] = ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0]*
288 (1. - oxy.poiii2 - oxy.poiii3) + oxy.d5007r*oxy.poiii2 +
289 oxy.d4363*oxy.poiii3;
290
291 if( ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] > 1e-30 && dense.IonLow[ipOXYGEN] <= 2 )
292 {
293 if( (oxy.d5007r*oxy.poiii2 + oxy.d4363*oxy.poiii3)/
294 ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] > (oxy.r4363Max +
295 oxy.r5007Max) )
296 {
297 oxy.poiii2Max = (realnum)(oxy.d5007r*oxy.poiii2/ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0]);
298 oxy.poiii3Max = (realnum)(oxy.d4363*oxy.poiii3/ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0]);
299 }
300 oxy.r4363Max = (realnum)(MAX2(oxy.r4363Max,oxy.d4363));
301 oxy.r5007Max = (realnum)(MAX2(oxy.r5007Max,oxy.d5007r));
302 }
303
304 /* ct into excited states */
305 if( dense.IonLow[ipOXYGEN] <= 0 && (ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0] +
306 atmdat.CharExcIonOf[ipHYDROGEN][ipOXYGEN][0]*dense.xIonDense[ipHYDROGEN][1]) > 1e-30 )
307 {
308 oxy.poimax = (realnum)(MAX2(oxy.poimax,oxy.d6300*oxy.poiexc/
309 (ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0]+
310 atmdat.CharExcIonOf[ipHYDROGEN][ipOXYGEN][0]* dense.xIonDense[ipHYDROGEN][1])));
311 }
312 }
313 else
314 {
315 oxy.poiii2Max = 0.;
316 oxy.poiii3Max = 0.;
317 oxy.r4363Max = 0.;
318 oxy.r5007Max = 0.;
319 oxy.poimax = 0.;
320 }
321 }
322 long int iup;
323 /* save atomic oxygen photodistruction rate for 3727 creation */
324 if( dense.IonLow[ipOXYGEN] == 0 && oxy.i2d < rfield.nflux )
325 {
326 oxy.s3727 = (realnum)(GammaK(oxy.i2d,oxy.i2p,opac.iopo2d , 1., &dummy ));
327
328 iup = MIN2(iso_sp[ipH_LIKE][1].fb[0].ipIsoLevNIonCon,rfield.nflux);
329 oxy.s7325 = (realnum)(GammaK(oxy.i2d,iup,opac.iopo2d , 1., &dummy ));
330
331 oxy.s7325 -= oxy.s3727;
332 oxy.s3727 = oxy.s3727 + oxy.s7325;
333
334 /* ratio of cross sections */
335 oxy.s7325 *= 0.66f;
336 }
337 else
338 {
339 oxy.s3727 = 0.;
340 oxy.s7325 = 0.;
341 }
342
343 oxy.AugerO3 = (realnum)ionbal.PhotoRate_Shell[ipOXYGEN][0][0][0];
344
345 oxy.s3727 *= dense.xIonDense[ipOXYGEN][0];
346 oxy.s7325 *= dense.xIonDense[ipOXYGEN][0];
347 oxy.AugerO3 *= dense.xIonDense[ipOXYGEN][0];
348 }
349 if( nelem == ipIRON )
350 {
351 if( !dense.lgElmtOn[ipIRON] )
352 {
353 fe.fekcld = 0.;
354 fe.fekhot = 0.;
355 fe.fegrain = 0.;
356 }
357 else
358 {
359 const int NDIM = ipIRON+1;
360
361 static const double fyield[NDIM+1] = {.34,.34,.35,.35,.36,.37,.37,.38,.39,.40,
362 .41,.42,.43,.44,.45,.46,.47,.47,.48,.48,.49,.49,.11,.75,0.,0.,0.};
363
364 long int i, limit, limit2;
365 /* now find total Auger yield of K-alphas
366 * "cold" iron has M-shell electrons, up to Fe 18 */
367 fe.fekcld = 0.;
368 limit = MIN2(18,dense.IonHigh[ipIRON]);
369
370 for( i=dense.IonLow[ipIRON]; i < limit; i++ )
371 {
372 ASSERT( i < NDIM + 1 );
373 fe.fekcld +=
374 (realnum)(ionbal.PhotoRate_Shell[ipIRON][i][0][0]*dense.xIonDense[ipIRON][i]*
375 fyield[i]);
376 }
377
378 /* same sum for hot iron */
379 fe.fekhot = 0.;
380 limit = MAX2(18,dense.IonLow[ipIRON]);
381
382 limit2 = MIN2(ipIRON+1,dense.IonHigh[ipIRON]);
383 ASSERT( limit2 <= LIMELM + 1 );
384
385 for( i=limit; i < limit2; i++ )
386 {
387 ASSERT( i < NDIM + 1 );
388 fe.fekhot +=
389 (realnum)(ionbal.PhotoRate_Shell[ipIRON][i][0][0]*dense.xIonDense[ipIRON][i]*
390 fyield[i]);
391 }
392
393 /* Fe Ka from grains - Fe in grains assumed to be atomic
394 * gv.elmSumAbund[ipIRON] is number density of iron added over all grain species */
395 i = 0;
396 /* fyield is 0.34 for atomic fe */
397 fe.fegrain = ( gv.lgWD01 ) ? 0.f : (realnum)(ionbal.PhotoRate_Shell[ipIRON][i][0][0]*fyield[i]*
398 gv.elmSumAbund[ipIRON]);
399 }
400 }
401
402 return;
403}
t_atmdat atmdat
Definition atmdat.cpp:6
double atom_pop2(double omega, double g1, double g2, double a21, double bltz, double abund)
Definition atom_pop2.cpp:9
t_atoms atoms
Definition atoms.cpp:5
t_ca ca
Definition ca.cpp:5
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
const int ipNITROGEN
Definition cddefines.h:311
const int ipIRON
Definition cddefines.h:330
const int ipCARBON
Definition cddefines.h:310
#define MIN2
Definition cddefines.h:761
const int ipOXYGEN
Definition cddefines.h:312
const int LIMELM
Definition cddefines.h:258
const int ipMAGNESIUM
Definition cddefines.h:316
const int NISO
Definition cddefines.h:261
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
const int ipCALCIUM
Definition cddefines.h:324
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
#define HMRATE(a, b, c)
Definition cddefines.h:1046
static t_yield & Inst()
Definition cddefines.h:175
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)
t_conv conv
Definition conv.cpp:5
t_dense dense
Definition dense.cpp:24
t_elementnames elementnames
t_fe fe
Definition fe.cpp:5
GrainVar gv
Definition grainvar.cpp:5
t_Heavy Heavy
Definition heavy.cpp:5
t_hmi hmi
Definition hmi.cpp:5
void ion_photo(long int nelem, bool lgPrintIt)
Definition ion_photo.cpp:26
t_ionbal ionbal
Definition ionbal.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH_LIKE
Definition iso.h:62
t_mole_global mole_global
Definition mole.cpp:6
t_opac opac
Definition opacity.cpp:5
t_oxy oxy
Definition oxy.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_thermal thermal
Definition thermal.cpp:5