cloudy trunk
Loading...
Searching...
No Matches
cont_gammas.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/* GammaBn evaluate photoionization rate for single shell with induced recomb */
4/* GammaPrt special version of gamma function to print strong contributors */
5/* GammaK evaluate photoionization rate for single shell */
6/* GammaPrtRate print photo rates for all shells of a ion and element */
7/*GammaPrtShells for the element nelem and ion, print total photo rate, subshells,
8 * and call GamaPrt for important subshells */
9#include "cddefines.h"
10#include "physconst.h"
11#include "iso.h"
12#include "thermal.h"
13#include "secondaries.h"
14#include "opacity.h"
15#include "rfield.h"
16#include "ionbal.h"
17#include "atmdat.h"
18#include "heavy.h"
19#include "gammas.h"
20
21/*
22 * these are the routines that evaluate the photoionization rates, gammas,
23 * throughout cloudy. a considerable amount of time is spent in these routines,
24 * so they should be compiled at the highest possible efficientcy.
25 */
26
27/* removed these two unused routines around r3500 (late October, 2009). */
28/* GammaBnPL evaluate photoionization rate for single shell with induced recomb */
29/* GammaPL evaluate photoionization rate for power law photo cross section */
30
31
32/*>>chng 99 apr 16, ConInterOut had not been included in the threshold point for any
33 * of these routines. added it. also moved thresholds above loop for a few */
34
35double GammaBn(
36 /* index for low energy */
37 long int ipLoEnr,
38 /* index for high energy */
39 long int ipHiEnr,
40 /* offset within the opacity stack */
41 long int ipOpac,
42 /* threshold (Ryd) for arbitrary photoionization */
43 double thresh,
44 /* induced rec rate */
45 double *ainduc,
46 /* rec cooling */
47 double *rcool,
48 t_phoHeat *photoHeat)
49{
50 long int i,
51 ilo,
52 iup,
53 limit;
54 double GamHi,
55 bnfun_v,
56 emin,
57 g,
58 phisig,
59 RateInducRec,
60 RateInducRecCool,
61 prod;
62
63 DEBUG_ENTRY( "GammaBn()" );
64
65/**********************************************************************
66 *
67 * special version of GAMFUN for arbitrary threshold, and induc rec
68 * compute ainduc, rate of inducted recombinations, rcool, induc rec cool
69 *
70 **********************************************************************/
71
72 /* special version of GAMFUN for arbitrary threshold, and induc rec
73 * compute ainduc, rate of inducted recombinations, rcool, induc rec cool */
74 if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr )
75 {
76 bnfun_v = 0.;
77 photoHeat->HeatNet = 0.;
78 photoHeat->HeatLowEnr = 0.;
79 photoHeat->HeatHiEnr = 0.;
80 *ainduc = 0.;
81 *rcool = 0.;
82 return( bnfun_v );
83 }
84
85 ASSERT( ipLoEnr >= 0 && ipHiEnr >= 0 );
86
87 /* >>chng 96 oct 25, first photo elec energy may have been negative
88 * possible for first any point to be below threshold due to
89 * finite resolution of continuum mesh */
90 emin = MAX2(thresh,rfield.anu[ipLoEnr-1]);
91 /* >>chng 99 jun 08 heating systematically too small, change to correct
92 * threshold, and protect first cell */
93 emin = thresh;
94
95 photoHeat->HeatNet = 0.;
96 g = 0.;
97 RateInducRec = 0.;
98 RateInducRecCool = 0.;
99
100 /* do first point without otscon, which may have diffuse cont,
101 * extra minus one after ipOpac is correct since not present in i = */
102 i = ipLoEnr;
103 /* >>chng 99 apr 16, add ConInterOut */
104 /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */
105 phisig = (rfield.flux[0][i-1] + rfield.otslin[i-1] +
106 rfield.ConInterOut[i-1]*rfield.lgOutOnly)*
107 opac.OpacStack[i-ipLoEnr+ipOpac-1];
108
109 g += phisig;
110 photoHeat->HeatNet += phisig*rfield.anu[i-1];
111
112 /* integral part of induced rec rate */
113 prod = phisig*rfield.ContBoltz[i-1];
114 RateInducRec += prod;
115 /* incuded rec cooling */
116 RateInducRecCool += prod*(rfield.anu[i-1] - emin);
117
118 iup = MIN2(ipHiEnr,rfield.nflux);
119 limit = MIN2(iup,secondaries.ipSecIon-1);
120
121 /* these is no -1 after the ipLoEnr in the following, due to c off by one counting */
122 for( i=ipLoEnr; i < limit; i++ )
123 {
124 /* SummedCon defined in RT_OTS_Update,
125 * includes flux, otscon, otslin, ConInterOut, outlin */
126 phisig = rfield.SummedCon[i]*opac.OpacStack[i-ipLoEnr+ipOpac];
127
128 g += phisig;
129 photoHeat->HeatNet += phisig*rfield.anu[i];
130
131 /* phi-sigma times exp(-hnu/kT) */
132 prod = phisig*rfield.ContBoltz[i];
133 /* induced recombination rate */
134 RateInducRec += prod;
135 /* incuded rec cooling */
136 RateInducRecCool += prod*(rfield.anu[i] - emin);
137 }
138
139 /* heating in Rydbergs, no secondary ionization */
140 /* >>chng 02 mar 31, added MAX2 due to roundoff error
141 * creating slightly negative number, caught downstream as insanity */
142 photoHeat->HeatNet = MAX2(0.,photoHeat->HeatNet - g*thresh);
143
144 /* we will save this heat - the part that does not secondary ionize */
145 photoHeat->HeatLowEnr = photoHeat->HeatNet;
146
147 /* now do part where secondary ionization is possible */
148 photoHeat->HeatHiEnr = 0.;
149 GamHi = 0.;
150
151 /* make sure we don't count twice when ipSecIon = ipLoEnr */
152 ilo = MAX2(ipLoEnr+1,secondaries.ipSecIon);
153 for( i=ilo-1; i < iup; i++ )
154 {
155 phisig = rfield.SummedCon[i]*opac.OpacStack[i-ipLoEnr+ipOpac];
156
157 GamHi += phisig;
158 photoHeat->HeatHiEnr += phisig*rfield.anu[i];
159
160 /* integral part of induced rec rate */
161 prod = phisig*rfield.ContBoltz[i];
162 RateInducRec += prod;
163 /* incuded rec cooling */
164 RateInducRecCool += prod*(rfield.anu[i] - emin);
165 }
166
167 /* this is total photo rate, low and high energy parts */
168 bnfun_v = g + GamHi;
169
170 /* heating in Rydbergs, uncorrected for secondary ionization */
171 photoHeat->HeatHiEnr = photoHeat->HeatHiEnr - GamHi*thresh;
172
173 /* add corrected high energy heat to total */
174 photoHeat->HeatNet += photoHeat->HeatHiEnr*secondaries.HeatEfficPrimary;
175
176 /* final product is heating in erg */
177 photoHeat->HeatNet *= EN1RYD;
178 photoHeat->HeatHiEnr *= EN1RYD;
179 photoHeat->HeatLowEnr *= EN1RYD;
180
181 /* this is an option to turn off induced processes */
182 if( rfield.lgInducProcess )
183 {
184 *rcool = RateInducRecCool*EN1RYD;
185 *ainduc = RateInducRec;
186 }
187 else
188 {
189 /* the "no induced" command was given */
190 *rcool = 0.;
191 *ainduc = 0.;
192 }
193
194 ASSERT( bnfun_v >= 0. );
195 ASSERT( photoHeat->HeatNet>= 0. );
196 return( bnfun_v );
197}
198
199/*GammaPrtShells for the element nelem and ion, print total photo rate, subshells,
200 * and call GamaPrt for important subshells */
201void GammaPrtShells( long nelem , long ion )
202{
203 double sum;
204 long int ns;
205
206 DEBUG_ENTRY( "GammaPrtShells()" );
207
208 fprintf(ioQQQ," GammaPrtShells nz\t%.2f \t%.2li %.2li ",
209 fnzone ,
210 nelem,
211 ion
212 );
213 sum = 0;
214 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ++ns )
215 {
216 sum += ionbal.PhotoRate_Shell[nelem][ion][ns][0];
217 }
218 fprintf(ioQQQ,"\ttot\t%.2e", sum);
219 /* loop over all shells for this ion */
220 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ++ns )
221 {
222 fprintf(ioQQQ,"\t%.2e", ionbal.PhotoRate_Shell[nelem][ion][ns][0] );
223# if 0
224 /* always reevaluate the outer shell, and all shells if lgRedoStatic is set */
225 if( (ns==(Heavy.nsShells[nelem][ion]-1) || opac.lgRedoStatic) )
226 {
227 /* option to redo the rates only on occasion */
228 iplow = opac.ipElement[nelem][ion][ns][0];
229 iphi = opac.ipElement[nelem][ion][ns][1];
230 ipop = opac.ipElement[nelem][ion][ns][2];
231
232 /* compute the photoionization rate, ionbal.lgPhotoIoniz_On is 1, set 0
233 * with "no photoionization" command */
234 ionbal.PhotoRate_Shell[nelem][ion][ns][0] =
235 GammaK(iplow,iphi,
236 ipop,t_yield::Inst().elec_eject_frac(nelem,ion,ns,0),
237 &phoHeat)*ionbal.lgPhotoIoniz_On;
238 }
239# endif
240 }
241 fprintf(ioQQQ,"\n");
242 return;
243}
244
245/**********************************************************************
246 *
247 * GammaPrt special version of gamma function to print strong contributors
248 * this only prints info - does not update heating rates properly since
249 * this is only a diagnostic
250 *
251 **********************************************************************/
252
254 /* first three par are identical to GammaK */
255 long int ipLoEnr,
256 long int ipHiEnr,
257 long int ipOpac,
258 /* io unit we will write to */
259 FILE * ioFILE,
260 /* total photo rate from previous call to GammaK */
261 double total,
262 /* we will print contributors that are more than this rate */
263 double threshold)
264{
265 long int i,
266 j,
267 k;
268 double flxcor,
269 phisig;
270
271 DEBUG_ENTRY( "GammaPrt()" );
272
273 /* special special version of GAMFUN to save step-by-step results */
274 /* >>chng 99 apr 02, test for lower greater than higher (when shell
275 * does not exist) added. This caused incorrect photo rates for
276 * non-existant shells */
277 if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr )
278 {
279 return;
280 }
281
282 fprintf( ioFILE, " GammaPrt %.2f from ",fnzone);
283 fprintf( ioFILE,PrintEfmt("%9.2e",rfield.anu[ipLoEnr-1]));
284 fprintf( ioFILE, " to ");
285 fprintf( ioFILE,PrintEfmt("%9.2e",rfield.anu[ipHiEnr-1]));
286 fprintf( ioFILE, "R rates >");
287 fprintf( ioFILE,PrintEfmt("%9.2e",threshold));
288 fprintf( ioFILE, " of total=");
289 fprintf( ioFILE,PrintEfmt("%9.2e",total));
290 fprintf( ioFILE, " (frac inc, otslin, otscon, ConInterOut, outlin ConOTS_local_OTS_rate ) chL, C\n");
291
292 if( threshold <= 0. || total <= 0. )
293 {
294 return;
295 }
296
297 k = ipLoEnr;
298 j = MIN2(ipHiEnr,rfield.nflux);
299
300 /* do theshold special, do not pick up otscon */
301 i = k-1;
302
303 /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */
304 phisig = (rfield.flux[0][i] + rfield.otslin[i]+ rfield.ConInterOut[i]*rfield.lgOutOnly)*
305 opac.OpacStack[i-k+ipOpac];
306 if( phisig > threshold || phisig < 0.)
307 {
308 flxcor = rfield.flux[0][i] + rfield.otslin[i] + rfield.ConInterOut[i]*rfield.lgOutOnly;
309
310 /* this really is array index on C scale */
311 fprintf( ioFILE, "[%5ld]" , i );
312 fprintf( ioFILE, PrintEfmt("%9.2e",rfield.anu[i]));
313 fprintf( ioFILE, PrintEfmt("%9.2e",phisig/total));
314 fprintf( ioFILE, "%5.2f%5.2f%5.2f%5.2f%5.2f%5.2f %4.4s %4.4s %.2e \n",
315 rfield.flux[0][i]/SDIV(flxcor),
316 rfield.otslin[i]/SDIV(flxcor),
317 /* otscon will appear below, but is not counted here, so do not print it (deceiving) */
318 0./SDIV(flxcor),
319 rfield.ConInterOut[i]*rfield.lgOutOnly/SDIV(flxcor),
320 (rfield.outlin[0][i]+rfield.outlin_noplot[i])/SDIV(flxcor),
321 rfield.ConOTS_local_OTS_rate[i]/SDIV(flxcor),
322 rfield.chLineLabel[i],
323 rfield.chContLabel[i],
324 opac.OpacStack[i-k+ipOpac]);
325 }
326
327 for( i=k; i < j; i++ )
328 {
329 phisig = rfield.SummedCon[i]*opac.OpacStack[i-k+ipOpac];
330 if( phisig > threshold || phisig < 0.)
331 {
332 /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */
333 flxcor = rfield.flux[0][i] + rfield.otslin[i] + rfield.otscon[i] +
334 rfield.outlin[0][i] + rfield.outlin_noplot[i] +rfield.ConInterOut[i]*rfield.lgOutOnly;
335
336 fprintf( ioFILE, "%5ld", i );
337 fprintf(ioFILE,PrintEfmt("%9.2e",rfield.anu[i]));
338 fprintf(ioFILE,PrintEfmt("%9.2e",phisig/total));
339 fprintf( ioFILE, "%5.2f%5.2f%5.2f%5.2f%5.2f%5.2f %4.4s %4.4s %.2e \n",
340 rfield.flux[0][i]/SDIV(flxcor),
341 rfield.otslin[i]/SDIV(flxcor),
342 rfield.otscon[i]/SDIV(flxcor),
343 rfield.ConInterOut[i]*rfield.lgOutOnly/SDIV(flxcor),
344 (rfield.outlin[0][i] + rfield.outlin_noplot[i])/SDIV(flxcor),
345 rfield.ConOTS_local_OTS_rate[i]/SDIV(flxcor),
346 rfield.chLineLabel[i],
347 rfield.chContLabel[i],
348 opac.OpacStack[i-k+ipOpac]);
349 }
350 }
351 return;
352}
353
354/*GammaK evaluate photoionization rate for single shell */
355
356/* this routine is a major pace setter for the code
357 * carefully check anything put into the following do loop */
358
359double GammaK(
360 /* ipLoEnr and ipHiEnr are pointers within frequency array to upper and
361 * lower bounds of evaluation */
362 long int ipLoEnr,
363 long int ipHiEnr,
364 /* ipOpac is offset to map onto OPSV*/
365 long int ipOpac,
366 /* yield1 is fraction of ionizations that emit 1 electron only,
367 * only used for energy balance */
368 double yield1,
369 t_phoHeat *photoHeat)
370{
371 long int i,
372 ilo,
373 ipOffset,
374 iup,
375 limit;
376 double GamHi,
377 eauger;
378
379 double gamk_v ,
380 phisig;
381
382 DEBUG_ENTRY( "GammaK()" );
383
384 /* evaluate photoioinzation rate and photoelectric heating
385 * returns photoionization rate as GAMK, heating is H */
386
387 /* >>chng 99 apr 02, test for lower greater than higher (when shell
388 * does not exist) added. This caused incorrect photo rates for
389 * non-existant shells */
390 if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr)
391 {
392 gamk_v = 0.;
393 photoHeat->HeatNet = 0.;
394 photoHeat->HeatHiEnr = 0.;
395 photoHeat->HeatLowEnr = 0.;
396 return( gamk_v );
397 }
398
399 iup = MIN2(ipHiEnr,rfield.nflux);
400
401 /* anu(i) is threshold, assume each auger electron has this energy
402 * less threshold energy, IP lost in initial photoionizaiton
403 * yield1 is the fraction of photos that emit 1 electron */
404 eauger = rfield.anu[ipLoEnr-1]*yield1;
405
406 /* low energies where secondary ionization cannot occur
407 * will do threshold point, ipLoEnr, later */
408 gamk_v = 0.;
409 photoHeat->HeatNet = 0.;
410
411 /* set up total offset for pointer addition not in loop */
412 ipOffset = ipOpac - ipLoEnr;
413
414 /* >>>chng 99 apr 16, this had followed the loop below, moved here to
415 * be like rest of gamma functions */
416 /* first do the threshold point
417 * do not include otscon, which may contain diffuse continuum */
418 /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */
419 phisig = (rfield.flux[0][ipLoEnr-1] + rfield.otslin[ipLoEnr-1]+
420 rfield.ConInterOut[ipLoEnr-1]*rfield.lgOutOnly) * opac.OpacStack[ipOpac-1];
421 gamk_v += phisig;
422 photoHeat->HeatNet += phisig*rfield.anu[ipLoEnr-1];
423
424 /* this loop only executed for energies than cannot sec ioniz */
425 limit = MIN2(iup,secondaries.ipSecIon-1);
426 for( i=ipLoEnr; i < limit; i++ )
427 {
428 phisig = rfield.SummedCon[i]*opac.OpacStack[i+ipOffset];
429 gamk_v += phisig;
430 photoHeat->HeatNet += phisig*rfield.anu[i];
431 }
432
433 /* correct heating for work function */
434 /* >>chng 02 apr 10, from first to sec, due to neg heat in blister.in */
435 /* *photoHeat->HeatNet += -gamk_v*eauger;*/
436 ASSERT( photoHeat->HeatNet >= 0. );
437 photoHeat->HeatNet = MAX2(0.,photoHeat->HeatNet - gamk_v*eauger);
438
439 /* this is the total low-energy heating, in ryd, that cannot secondary ionize */
440 photoHeat->HeatLowEnr = photoHeat->HeatNet;
441
442 /* now do part where secondary ionization possible */
443 photoHeat->HeatHiEnr = 0.;
444 GamHi = 0.;
445 /* make sure we don't count twice when ipSecIon = ipLoEnr */
446 ilo = MAX2(ipLoEnr+1,secondaries.ipSecIon);
447 for( i=ilo-1; i < iup; i++ )
448 {
449 /* produce of flux and cross section */
450 phisig = rfield.SummedCon[i]*opac.OpacStack[i+ipOffset];
451 GamHi += phisig;
452 photoHeat->HeatHiEnr += phisig*rfield.anu[i];
453 }
454
455 /* add on the high energy part */
456 gamk_v += GamHi;
457 /* correct for work function */
458 photoHeat->HeatHiEnr -= GamHi*eauger;
459
460 /* net heating include high energy heat, with correction for
461 * secondary ionization */
462 photoHeat->HeatNet += photoHeat->HeatHiEnr*secondaries.HeatEfficPrimary;
463
464 /* final product is heating in erg */
465 photoHeat->HeatNet *= EN1RYD;
466 photoHeat->HeatLowEnr *= EN1RYD;
467 photoHeat->HeatHiEnr *= EN1RYD;
468
469 ASSERT( gamk_v >= 0. );
470 ASSERT( photoHeat->HeatNet>= 0. );
471 return( gamk_v );
472}
473
474/* GammaPrtRate will print resulting rates for ion and element */
476 /* io unit we will write to */
477 FILE * ioFILE,
478 /* stage of ionization on C scale, 0 for atom */
479 long int ion ,
480 /* 0 for H, etc */
481 long int nelem,
482 /* true - then print photo sources for each shell */
483 bool lgPRT )
484{
485 long int nshell , ns;
486
487 DEBUG_ENTRY( "GammaPrtRate()" );
488
489 /* number of shells for this species */
490 nshell = Heavy.nsShells[nelem][ion];
491
492 /* now print subshell photo rate */
493 fprintf(ioFILE , "GammaPrtRate: %li %li",ion , nelem );
494 for( ns=nshell-1; ns>=0; --ns )
495 {
496 fprintf(ioFILE , " %.2e" , ionbal.PhotoRate_Shell[nelem][ion][ns][0]);
497
498 /* option to print individual contributors to each shell */
499 if( lgPRT )
500 {
501 fprintf(ioFILE , "\n");
502 GammaPrt(
503 /* first three par are identical to GammaK */
504 opac.ipElement[nelem][ion][ns][0],
505 opac.ipElement[nelem][ion][ns][1],
506 opac.ipElement[nelem][ion][ns][2],
507 /* io unit we will write to */
508 ioFILE,
509 /* total photo rate from previous call to GammaK */
510 ionbal.PhotoRate_Shell[nelem][ion][ns][0],
511 /* we will print contributors that are more than this rate */
512 ionbal.PhotoRate_Shell[nelem][ion][ns][0]*0.05);
513 }
514 }
515 fprintf(ioFILE , "\n");
516 return;
517}
FILE * ioQQQ
Definition cddefines.cpp:7
double fnzone
Definition cddefines.cpp:15
#define ASSERT(exp)
Definition cddefines.h:578
#define PrintEfmt(F, V)
Definition cddefines.h:1472
#define MIN2
Definition cddefines.h:761
#define MAX2
Definition cddefines.h:782
sys_float SDIV(sys_float x)
Definition cddefines.h:952
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
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)
void GammaPrtRate(FILE *ioFILE, long int ion, long int nelem, bool lgPRT)
double GammaK(long int ipLoEnr, long int ipHiEnr, long int ipOpac, double yield1, t_phoHeat *photoHeat)
void GammaPrtShells(long nelem, long ion)
double GammaBn(long int ipLoEnr, long int ipHiEnr, long int ipOpac, double thresh, double *ainduc, double *rcool, t_phoHeat *photoHeat)
t_Heavy Heavy
Definition heavy.cpp:5
t_ionbal ionbal
Definition ionbal.cpp:5
t_opac opac
Definition opacity.cpp:5
UNUSED const double EN1RYD
Definition physconst.h:179
t_rfield rfield
Definition rfield.cpp:8
t_secondaries secondaries
static double * g
Definition species2.cpp:28
double HeatNet
Definition thermal.h:172
double HeatHiEnr
Definition thermal.h:176
double HeatLowEnr
Definition thermal.h:174