cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cont_gammas.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 /* GammaBn evaluate photoionization rate for single shell with induced recomb */
4 /* GammaBnPL evaluate photoionization rate for single shell with induced recomb */
5 /* GammaPrt special version of gamma function to print strong contributors */
6 /* GammaK evaluate photoionization rate for single shell */
7 /* GammaPL evaluate photoionization rate for power law photo cross section */
8 /* GammaPrtRate print photo rates for all shells of a ion and element */
9 /*GammaPrtShells for the element nelem and ion, print total photo rate, subshells,
10  * and call GamaPrt for important subshells */
11 #include "cddefines.h"
12 #include "physconst.h"
13 #include "iso.h"
14 #include "thermal.h"
15 #include "secondaries.h"
16 #include "opacity.h"
17 #include "rfield.h"
18 #include "ionbal.h"
19 #include "atmdat.h"
20 #include "heavy.h"
21 #include "gammas.h"
22 
23 /*
24  * these are the routines that evaluate the photoionization rates, gammas,
25  * throughout cloudy. a considerable amount of time is spent in these routines,
26  * so they should be compiled at the highest possible efficientcy.
27  */
28 
29 /*>>chng 99 apr 16, ConInterOut had not been included in the threshold point for any
30  * of these routines. added it. also moved thresholds above loop for a few */
31 
32 double GammaBn(
33  /* index for low energy */
34  long int ipLoEnr,
35  /* index for high energy */
36  long int ipHiEnr,
37  /* offset within the opacity stack */
38  long int ipOpac,
39  /* threshold (Ryd) for arbitrary photoionization */
40  double thresh,
41  /* induced rec rate */
42  double *ainduc,
43  /* rec cooling */
44  double *rcool)
45 {
46  long int i,
47  ilo,
48  iup,
49  limit;
50  double GamHi,
51  bnfun_v,
52  emin,
53  g,
54  phisig,
55  RateInducRec,
56  RateInducRecCool,
57  prod;
58 
59  DEBUG_ENTRY( "GammaBn()" );
60 
61 /**********************************************************************
62  *
63  * special version of GAMFUN for arbitrary threshold, and induc rec
64  * compute ainduc, rate of inducted recombinations, rcool, induc rec cool
65  *
66  **********************************************************************/
67 
68  /* special version of GAMFUN for arbitrary threshold, and induc rec
69  * compute ainduc, rate of inducted recombinations, rcool, induc rec cool */
70  if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr )
71  {
72  bnfun_v = 0.;
73  thermal.HeatNet = 0.;
74  thermal.HeatLowEnr = 0.;
75  thermal.HeatHiEnr = 0.;
76  *ainduc = 0.;
77  *rcool = 0.;
78  return( bnfun_v );
79  }
80 
81  ASSERT( ipLoEnr >= 0 && ipHiEnr >= 0 );
82 
83  /* >>chng 96 oct 25, first photo elec energy may have been negative
84  * possible for first any point to be below threshold due to
85  * finite resolution of continuum mesh */
86  emin = MAX2(thresh,rfield.anu[ipLoEnr-1]);
87  /* >>chng 99 jun 08 heating systematically too small, change to correct
88  * threshold, and protect first cell */
89  emin = thresh;
90 
91  thermal.HeatNet = 0.;
92  g = 0.;
93  RateInducRec = 0.;
94  RateInducRecCool = 0.;
95 
96  /* do first point without otscon, which may have diffuse cont,
97  * extra minus one after ipOpac is correct since not present in i = */
98  i = ipLoEnr;
99  /* >>chng 99 apr 16, add ConInterOut */
100  /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */
101  phisig = (rfield.flux[i-1] + rfield.otslin[i-1] +
103  opac.OpacStack[i-ipLoEnr+ipOpac-1];
104 
105  g += phisig;
106  thermal.HeatNet += phisig*rfield.anu[i-1];
107 
108  /* integral part of induced rec rate */
109  prod = phisig*rfield.ContBoltz[i-1];
110  RateInducRec += prod;
111  /* incuded rec cooling */
112  RateInducRecCool += prod*(rfield.anu[i-1] - emin);
113 
114  iup = MIN2(ipHiEnr,rfield.nflux);
115  limit = MIN2(iup,secondaries.ipSecIon-1);
116 
117  /* these is no -1 after the ipLoEnr in the following, due to c off by one counting */
118  for( i=ipLoEnr; i < limit; i++ )
119  {
120  /* SummedCon defined in RT_OTS_Update,
121  * includes flux, otscon, otslin, ConInterOut, outlin */
122  phisig = rfield.SummedCon[i]*opac.OpacStack[i-ipLoEnr+ipOpac];
123 
124  g += phisig;
125  thermal.HeatNet += phisig*rfield.anu[i];
126 
127  /* phi-sigma times exp(-hnu/kT) */
128  prod = phisig*rfield.ContBoltz[i];
129  /* induced recombination rate */
130  RateInducRec += prod;
131  /* incuded rec cooling */
132  RateInducRecCool += prod*(rfield.anu[i] - emin);
133  }
134 
135  /* heating in Rydbergs, no secondary ionization */
136  /* >>chng 02 mar 31, added MAX2 due to roundoff error
137  * creating slightly negative number, caught downstream as insanity */
138  thermal.HeatNet = MAX2(0.,thermal.HeatNet - g*thresh);
139 
140  /* we will save this heat - the part that does not secondary ionize */
142 
143  /* now do part where secondary ionization is possible */
144  thermal.HeatHiEnr = 0.;
145  GamHi = 0.;
146 
147  /* make sure we don't count twice when ipSecIon = ipLoEnr */
148  ilo = MAX2(ipLoEnr+1,secondaries.ipSecIon);
149  for( i=ilo-1; i < iup; i++ )
150  {
151  phisig = rfield.SummedCon[i]*opac.OpacStack[i-ipLoEnr+ipOpac];
152 
153  GamHi += phisig;
154  thermal.HeatHiEnr += phisig*rfield.anu[i];
155 
156  /* integral part of induced rec rate */
157  prod = phisig*rfield.ContBoltz[i];
158  RateInducRec += prod;
159  /* incuded rec cooling */
160  RateInducRecCool += prod*(rfield.anu[i] - emin);
161  }
162 
163  /* this is total photo rate, low and high energy parts */
164  bnfun_v = g + GamHi;
165 
166  /* heating in Rydbergs, uncorrected for secondary ionization */
167  thermal.HeatHiEnr = thermal.HeatHiEnr - GamHi*thresh;
168 
169  /* add corrected high energy heat to total */
171 
172  /* final product is heating in erg */
176 
177  /* this is an option to turn off induced processes */
178  if( rfield.lgInducProcess )
179  {
180  *rcool = RateInducRecCool*EN1RYD;
181  *ainduc = RateInducRec;
182  }
183  else
184  {
185  /* the "no induced" command was given */
186  *rcool = 0.;
187  *ainduc = 0.;
188  }
189 
190  ASSERT( bnfun_v >= 0. );
191  ASSERT( thermal.HeatNet>= 0. );
192  return( bnfun_v );
193 }
194 
195 /*GammaPrtShells for the element nelem and ion, print total photo rate, subshells,
196  * and call GamaPrt for important subshells */
197 void GammaPrtShells( long nelem , long ion )
198 {
199  double sum;
200  long int ns;
201 
202  DEBUG_ENTRY( "GammaPrtShells()" );
203 
204  fprintf(ioQQQ," GammaPrtShells nz\t%.2f \t%.2li %.2li ",
205  fnzone ,
206  nelem,
207  ion
208  );
209  sum = 0;
210  for( ns=0; ns < Heavy.nsShells[nelem][ion]; ++ns )
211  {
212  sum += ionbal.PhotoRate_Shell[nelem][ion][ns][0];
213  }
214  fprintf(ioQQQ,"\ttot\t%.2e", sum);
215  /* loop over all shells for this ion */
216  for( ns=0; ns < Heavy.nsShells[nelem][ion]; ++ns )
217  {
218  fprintf(ioQQQ,"\t%.2e", ionbal.PhotoRate_Shell[nelem][ion][ns][0] );
219 # if 0
220  /* always reevaluate the outer shell, and all shells if lgRedoStatic is set */
221  if( (ns==(Heavy.nsShells[nelem][ion]-1) || opac.lgRedoStatic) )
222  {
223  /* option to redo the rates only on occasion */
224  iplow = opac.ipElement[nelem][ion][ns][0];
225  iphi = opac.ipElement[nelem][ion][ns][1];
226  ipop = opac.ipElement[nelem][ion][ns][2];
227 
228  /* compute the photoionization rate, ionbal.lgPhotoIoniz_On is 1, set 0
229  * with "no photoionization" command */
230  ionbal.PhotoRate_Shell[nelem][ion][ns][0] =
231  GammaK(iplow,iphi,
232  ipop,t_yield::Inst().elec_eject_frac(nelem,ion,ns,0))*ionbal.lgPhotoIoniz_On;
233  }
234 # endif
235  }
236  fprintf(ioQQQ,"\n");
237  return;
238 }
239 
240 /**********************************************************************
241  *
242  * GammaPrt special version of gamma function to print strong contributors
243  * this only prints info - does not update heating rates properly since
244  * this is only a diagnostic
245  *
246  **********************************************************************/
247 
248 void GammaPrt(
249  /* first three par are identical to GammaK */
250  long int ipLoEnr,
251  long int ipHiEnr,
252  long int ipOpac,
253  /* io unit we will write to */
254  FILE * ioFILE,
255  /* total photo rate from previous call to GammaK */
256  double total,
257  /* we will print contributors that are more than this rate */
258  double threshold)
259 {
260  long int i,
261  j,
262  k;
263  double flxcor,
264  phisig;
265 
266  DEBUG_ENTRY( "GammaPrt()" );
267 
268  /* special special version of GAMFUN to punch step-by-step results */
269  /* >>chng 99 apr 02, test for lower greater than higher (when shell
270  * does not exist) added. This caused incorrect photo rates for
271  * non-existant shells */
272  if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr )
273  {
274  return;
275  }
276 
277  fprintf( ioFILE, " GammaPrt %.2f from ",fnzone);
278  fprintf( ioFILE,PrintEfmt("%9.2e",rfield.anu[ipLoEnr-1]));
279  fprintf( ioFILE, " to ");
280  fprintf( ioFILE,PrintEfmt("%9.2e",rfield.anu[ipHiEnr-1]));
281  fprintf( ioFILE, "R rates >");
282  fprintf( ioFILE,PrintEfmt("%9.2e",threshold));
283  fprintf( ioFILE, " of total=");
284  fprintf( ioFILE,PrintEfmt("%9.2e",total));
285  fprintf( ioFILE, " (frac inc, otslin, otscon, ConInterOut, outlin ConOTS_local_OTS_rate ) chL, C\n");
286 
287  if( threshold <= 0. || total <= 0. )
288  {
289  return;
290  }
291 
292  k = ipLoEnr;
293  j = MIN2(ipHiEnr,rfield.nflux);
294 
295  /* do theshold special, do not pick up otscon */
296  i = k-1;
297 
298  /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */
299  phisig = (rfield.flux[i] + rfield.otslin[i]+ rfield.ConInterOut[i]*rfield.lgOutOnly)*
300  opac.OpacStack[i-k+ipOpac];
301  if( phisig > threshold || phisig < 0.)
302  {
303  flxcor = rfield.flux[i] + rfield.otslin[i] + rfield.ConInterOut[i]*rfield.lgOutOnly;
304 
305  /* this really is array index on C scale */
306  fprintf( ioFILE, "[%5ld]" , i );
307  fprintf( ioFILE, PrintEfmt("%9.2e",rfield.anu[i]));
308  fprintf( ioFILE, PrintEfmt("%9.2e",phisig/total));
309  fprintf( ioFILE, "%5.2f%5.2f%5.2f%5.2f%5.2f%5.2f %4.4s %4.4s %.2e \n",
310  rfield.flux[i]/SDIV(flxcor),
311  rfield.otslin[i]/SDIV(flxcor),
312  /* otscon will appear below, but is not counted here, so do not print it (deceiving) */
313  0./SDIV(flxcor),
314  rfield.ConInterOut[i]*rfield.lgOutOnly/SDIV(flxcor),
315  (rfield.outlin[i]+rfield.outlin_noplot[i])/SDIV(flxcor),
316  rfield.ConOTS_local_OTS_rate[i]/SDIV(flxcor),
317  rfield.chLineLabel[i],
318  rfield.chContLabel[i],
319  opac.OpacStack[i-k+ipOpac]);
320  }
321 
322  for( i=k; i < j; i++ )
323  {
324  phisig = rfield.SummedCon[i]*opac.OpacStack[i-k+ipOpac];
325  if( phisig > threshold || phisig < 0.)
326  {
327  /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */
328  flxcor = rfield.flux[i] + rfield.otslin[i] + rfield.otscon[i] +
330 
331  fprintf( ioFILE, "%5ld", i );
332  fprintf(ioFILE,PrintEfmt("%9.2e",rfield.anu[i]));
333  fprintf(ioFILE,PrintEfmt("%9.2e",phisig/total));
334  fprintf( ioFILE, "%5.2f%5.2f%5.2f%5.2f%5.2f%5.2f %4.4s %4.4s %.2e \n",
335  rfield.flux[i]/SDIV(flxcor),
336  rfield.otslin[i]/SDIV(flxcor),
337  rfield.otscon[i]/SDIV(flxcor),
338  rfield.ConInterOut[i]*rfield.lgOutOnly/SDIV(flxcor),
339  (rfield.outlin[i] + rfield.outlin_noplot[i])/SDIV(flxcor),
340  rfield.ConOTS_local_OTS_rate[i]/SDIV(flxcor),
341  rfield.chLineLabel[i],
342  rfield.chContLabel[i],
343  opac.OpacStack[i-k+ipOpac]);
344  }
345  }
346  return;
347 }
348 
349 /*GammaK evaluate photoionization rate for single shell */
350 
351 /* this routine is a major pace setter for the code
352  * carefully check anything put into the following do loop */
353 
354 double GammaK(
355  /* ipLoEnr and ipHiEnr are pointers within frequency array to upper and
356  * lower bounds of evaluation */
357  long int ipLoEnr,
358  long int ipHiEnr,
359  /* ipOpac is offset to map onto OPSV*/
360  long int ipOpac,
361  /* yield1 is fraction of ionizations that emit 1 electron only,
362  * only used for energy balance */
363  double yield1)
364 {
365  long int i,
366  ilo,
367  ipOffset,
368  iup,
369  limit;
370  double GamHi,
371  eauger;
372 
373  double gamk_v ,
374  phisig;
375 
376  DEBUG_ENTRY( "GammaK()" );
377 
378  /* evaluate photoioinzation rate and photoelectric heating
379  * returns photoionization rate as GAMK, heating is H */
380 
381  /* >>chng 99 apr 02, test for lower greater than higher (when shell
382  * does not exist) added. This caused incorrect photo rates for
383  * non-existant shells */
384  if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr)
385  {
386  gamk_v = 0.;
387  thermal.HeatNet = 0.;
388  thermal.HeatHiEnr = 0.;
389  thermal.HeatLowEnr = 0.;
390  return( gamk_v );
391  }
392 
393  iup = MIN2(ipHiEnr,rfield.nflux);
394 
395  /* anu(i) is threshold, assume each auger electron has this energy
396  * less threshold energy, IP lost in initial photoionizaiton
397  * yield1 is the fraction of photos that emit 1 electron */
398  eauger = rfield.anu[ipLoEnr-1]*yield1;
399 
400  /* low energies where secondary ionization cannot occur
401  * will do threshold point, ipLoEnr, later */
402  gamk_v = 0.;
403  thermal.HeatNet = 0.;
404 
405  /* set up total offset for pointer addition not in loop */
406  ipOffset = ipOpac - ipLoEnr;
407 
408  /* >>>chng 99 apr 16, this had followed the loop below, moved here to
409  * be like rest of gamma functions */
410  /* first do the threshold point
411  * do not include otscon, which may contain diffuse continuum */
412  /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */
413  phisig = (rfield.flux[ipLoEnr-1] + rfield.otslin[ipLoEnr-1]+
414  rfield.ConInterOut[ipLoEnr-1]*rfield.lgOutOnly) * opac.OpacStack[ipOpac-1];
415  gamk_v += phisig;
416  thermal.HeatNet += phisig*rfield.anu[ipLoEnr-1];
417 
418  /* this loop only executed for energies than cannot sec ioniz */
419  limit = MIN2(iup,secondaries.ipSecIon-1);
420  for( i=ipLoEnr; i < limit; i++ )
421  {
422  phisig = rfield.SummedCon[i]*opac.OpacStack[i+ipOffset];
423  gamk_v += phisig;
424  thermal.HeatNet += phisig*rfield.anu[i];
425  }
426 
427  /* correct heating for work function */
428  /* >>chng 02 apr 10, from first to sec, due to neg heat in blister.in */
429  /*thermal.HeatNet += -gamk_v*eauger;*/
430  ASSERT( thermal.HeatNet >= 0. );
431  thermal.HeatNet = MAX2(0.,thermal.HeatNet - gamk_v*eauger);
432 
433  /* this is the total low-energy heating, in ryd, that cannot secondary ionize */
435 
436  /* now do part where secondary ionization possible */
437  thermal.HeatHiEnr = 0.;
438  GamHi = 0.;
439  /* make sure we don't count twice when ipSecIon = ipLoEnr */
440  ilo = MAX2(ipLoEnr+1,secondaries.ipSecIon);
441  for( i=ilo-1; i < iup; i++ )
442  {
443  /* produce of flux and cross section */
444  phisig = rfield.SummedCon[i]*opac.OpacStack[i+ipOffset];
445  GamHi += phisig;
446  thermal.HeatHiEnr += phisig*rfield.anu[i];
447  }
448 
449  /* add on the high energy part */
450  gamk_v += GamHi;
451  /* correct for work function */
452  thermal.HeatHiEnr -= GamHi*eauger;
453 
454  /* net heating include high energy heat, with correction for
455  * secondary ionization */
457 
458  /* final product is heating in erg */
462 
463  ASSERT( gamk_v >= 0. );
464  ASSERT( thermal.HeatNet>= 0. );
465  return( gamk_v );
466 }
467 
468 #if 0
469 //this is never used
470 /*GammaPL evaluate photoionization rate for power law photoionization cross section */
471 
472 double GammaPL(
473  /* this is prin quantum number for level */
474  long int n,
475  /* nelem = 0 for H */
476  long int nelem)
477 {
478  long int i,
479  ilo,
480  iup,
481  limit,
482  ipLoEnr,
483  ipHiEnr;
484  realnum GamHi,
485  csThresh,
486  gamPL_v,
487  phisig;
488 
489  DEBUG_ENTRY( "GammaPL()" );
490 
491  /* evaluate photoioinzation rate and photoelectric heating
492  * returns photoionization rate as GammaPL, heating is H
493  * n and nelem and principal quantum number and charge
494  */
495 
496  /* validate the input */
497  ASSERT( nelem >= 0);
498  ASSERT( nelem < LIMELM);
499  ASSERT( n >= 0);
500 
501  /* get pointer to photo threshold */
502  ipLoEnr = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][n];
503  ASSERT( ipLoEnr>0 );
504 
505  if( ipLoEnr >= rfield.nflux )
506  {
507  gamPL_v = 0.;
508  thermal.HeatNet = 0.;
509  thermal.HeatHiEnr = 0.;
510  thermal.HeatLowEnr = 0.;
511  return( gamPL_v );
512  }
513 
514  ipHiEnr = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s];
515  iup = MIN2(ipHiEnr,rfield.nflux);
516 
517  csThresh = t_ADfA::Inst().sth(n)*rfield.anu3[ipLoEnr-1]/ POW2(nelem+1.f);
518 
519  /* low energies where secondary ionization cannot occur
520  * will do threshold point, ipLoEnr, later */
521  gamPL_v = 0.;
522  thermal.HeatNet = 0.;
523 
524  /* first do the threshold point
525  * do not include otscon, which may contain diffuse continuum */
526  /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */
527  phisig = (rfield.flux[ipLoEnr-1] + rfield.otslin[ipLoEnr-1] +
528  rfield.ConInterOut[ipLoEnr-1]*rfield.lgOutOnly ) /rfield.anu3[ipLoEnr-1];
529  gamPL_v += phisig;
530 
531  /* this loop only executed for energies than cannot sec ioniz */
532  limit = MIN2(iup,secondaries.ipSecIon-1);
533  for( i=ipLoEnr; i < limit; i++ )
534  {
535  phisig = rfield.SummedCon[i]/rfield.anu3[i];
536  gamPL_v += phisig;
537  thermal.HeatNet += phisig*rfield.anu[i];
538  }
539 
540  /* convert photo rate into proper units with CS at threshold */
541  gamPL_v *= csThresh;
542 
543  /* correct heating for CS at threshold and work function */
544  thermal.HeatNet *= csThresh;
545  thermal.HeatNet += -gamPL_v*rfield.anu[ipLoEnr-1];
547 
548  /* now do part where secondary ionization possible */
549  thermal.HeatHiEnr = 0.;
550  GamHi = 0.;
551 
552  /* make sure we don't count twice when ipSecIon = ipLoEnr */
553  ilo = MAX2(ipLoEnr+1,secondaries.ipSecIon);
554  for( i=ilo-1; i < iup; i++ )
555  {
556  /* produce of flux and cross section */
557  phisig = rfield.SummedCon[i]/rfield.anu3[i];
558  GamHi += phisig;
559  thermal.HeatHiEnr += phisig*rfield.anu[i];
560  }
561 
562  /* add on the high energy part */
563  gamPL_v += GamHi*csThresh;
564 
565  /* correct heating for work function and convert to true cross section */
566  thermal.HeatHiEnr = (thermal.HeatHiEnr - GamHi*rfield.anu[ipLoEnr-1])*csThresh;
567 
568  /* net heating includes high energy heating corrected for secondary ionizations */
570 
571  /* final product is heating in erg */
575 
576  ASSERT( gamPL_v>= 0.);
577  return( gamPL_v );
578 }
579 #endif
580 
581 /*GammaBnPL evaluate hydrogenic photoionization rate for single level with induced recomb */
582 
583 double GammaBnPL(long int n,
584  long int nelem, /* 0 for H, etc */
585  double *ainduc,
586  double *rcool)
587 {
588  long int i,
589  ilo,
590  iup,
591  limit,
592  ipLoEnr,
593  ipHiEnr;
594  double GamHi,
595  bnfunPL_v,
596  csThresh,
597  emin,
598  g,
599  phisig,
600  RateInducRec,
601  RateInducRecCool,
602  prod,
603  thresh;
604 
605  DEBUG_ENTRY( "GammaBnPL()" );
606 
607  /* validate the incoming charge */
608  ASSERT( nelem >= 0);
609  ASSERT( nelem < LIMELM );
610 
611  /* special version of GAMFUN for arbitrary threshold, and induc rec
612  * compute ainduc, rate of inducted recombinations, rcool, induc rec cool */
613 
614  /* these are the upper and lower energy bounds, which we know from iso-sequence */
615  ipLoEnr = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][n];
616  ipHiEnr = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
617 
618  ilo = ipLoEnr;
619  iup = MIN2(ipHiEnr,rfield.nflux);
620  if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr )
621  {
622  bnfunPL_v = 0.;
623  thermal.HeatNet = 0.;
624  thermal.HeatHiEnr = 0.;
625  thermal.HeatLowEnr = 0.;
626  *ainduc = 0.;
627  *rcool = 0.;
628  return( bnfunPL_v );
629  }
630 
631  /* >>chng 96 oct 25, first photo elec energy may have been negative
632  * possible for first any point to be below threshold due to
633  * finite resolution of continuum mesh */
634  thresh = iso.xIsoLevNIonRyd[ipH_LIKE][nelem][n];
635  emin = MAX2(thresh,rfield.anu[ipLoEnr-1]);
636  /* >>chng 99 jun 08 heating systematically too small, change to correct
637  * threshold, and protect first cell */
638  emin = thresh;
639 
640  thermal.HeatNet = 0.;
641  g = 0.;
642  RateInducRec = 0.;
643  RateInducRecCool = 0.;
644 
645  /*csThresh = t_ADfA::Inst().sth(n)*POW3(thresh) / POW2(nelem+1.f);*/
646  csThresh = t_ADfA::Inst().sth(n)*rfield.anu3[ipLoEnr-1] / POW2(nelem+1.f);
647 
648  /* do first point without otscon, which may have diffuse cont */
649  /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */
650  phisig = (rfield.flux[ipLoEnr-1] + rfield.otslin[ipLoEnr-1]+
651  rfield.ConInterOut[ipLoEnr-1]*rfield.lgOutOnly)/
652  rfield.anu3[ipLoEnr-1];
653 
654  g += phisig;
655  thermal.HeatNet += phisig*rfield.anu[ipLoEnr-1];
656 
657  /* integral part of induced rec rate */
658  prod = phisig*rfield.ContBoltz[ipLoEnr-1];
659  RateInducRec += prod;
660  /* induced rec cooling */
661  RateInducRecCool += prod*(rfield.anu[ipLoEnr-1] - emin);
662 
663  limit = MIN2(iup,secondaries.ipSecIon-1);
664  for( i=ipLoEnr; i < limit; i++ )
665  {
666  phisig = rfield.SummedCon[i]/rfield.anu3[i];
667 
668  g += phisig;
669  thermal.HeatNet += phisig*rfield.anu[i];
670 
671  /* integral part of induced rec rate */
672  prod = phisig*rfield.ContBoltz[i];
673  RateInducRec += prod;
674  /* incuded rec cooling */
675  RateInducRecCool += prod*(rfield.anu[i] - emin);
676  }
677 
678  /* heating in Rydbergs, no secondary ionization */
679  thermal.HeatNet = MAX2(0.,(thermal.HeatNet - g*thresh)*csThresh);
681 
682  /* now do part where secondary ionization is possible */
683  thermal.HeatHiEnr = 0.;
684  GamHi = 0.;
685  /* make sure we don't count twice when ipSecIon = ipLoEnr */
686  ilo = MAX2(ipLoEnr+1,secondaries.ipSecIon);
687  for( i=ilo-1; i < iup; i++ )
688  {
689  phisig = rfield.SummedCon[i]/rfield.anu3[i];
690  GamHi += phisig;
691  thermal.HeatHiEnr += phisig*rfield.anu[i];
692 
693  /* integral part of induced rec rate */
694  prod = phisig*rfield.ContBoltz[i];
695  RateInducRec += prod;
696  /* incuded rec cooling */
697  RateInducRecCool += prod*(rfield.anu[i] - emin);
698  }
699  RateInducRec *= csThresh;
700  RateInducRecCool *= csThresh;
701 
702  /* this is total photo rate, low and high energy parts */
703  bnfunPL_v = (g + GamHi)*csThresh;
704 
705  /* heating in Rydbergs, uncorrected for secondary ionization */
706  thermal.HeatHiEnr = (thermal.HeatHiEnr - GamHi*thresh)*csThresh;
707 
708  /* net heating includes high energy heat corrected for secondary ionization */
710 
711  /* final product is heating in erg */
715 
716  /* this is an option to turn off induced processes */
717  if( rfield.lgInducProcess )
718  {
719  *rcool = RateInducRecCool*EN1RYD;
720  *ainduc = RateInducRec;
721  }
722  else
723  {
724  /* the "no induced" command was given */
725  *rcool = 0.;
726  *ainduc = 0.;
727  }
728 
729  ASSERT( bnfunPL_v>= 0.);
730  return( bnfunPL_v );
731 }
732 
733 /* GammaPrtRate will print resulting rates for ion and element */
735  /* io unit we will write to */
736  FILE * ioFILE,
737  /* stage of ionization on C scale, 0 for atom */
738  long int ion ,
739  /* 0 for H, etc */
740  long int nelem,
741  /* true - then print photo sources for each shell */
742  bool lgPRT )
743 {
744  long int nshell , ns;
745 
746  DEBUG_ENTRY( "GammaPrtRate()" );
747 
748  /* number of shells for this species */
749  nshell = Heavy.nsShells[nelem][ion];
750 
751  /* now print subshell photo rate */
752  fprintf(ioFILE , "GammaPrtRate: %li %li",ion , nelem );
753  for( ns=nshell-1; ns>=0; --ns )
754  {
755  fprintf(ioFILE , " %.2e" , ionbal.PhotoRate_Shell[nelem][ion][ns][0]);
756 
757  /* option to print individual contributors to each shell */
758  if( lgPRT )
759  {
760  fprintf(ioFILE , "\n");
761  GammaPrt(
762  /* first three par are identical to GammaK */
763  opac.ipElement[nelem][ion][ns][0],
764  opac.ipElement[nelem][ion][ns][1],
765  opac.ipElement[nelem][ion][ns][2],
766  /* io unit we will write to */
767  ioFILE,
768  /* total photo rate from previous call to GammaK */
769  ionbal.PhotoRate_Shell[nelem][ion][ns][0],
770  /* we will print contributors that are more than this rate */
771  ionbal.PhotoRate_Shell[nelem][ion][ns][0]*0.05);
772  }
773  }
774  fprintf(ioFILE , "\n");
775  return;
776 }

Generated for cloudy by doxygen 1.8.1.2