cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
prt_lines_continuum.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 /*lines_continuum put energetics, H, and He lines into line intensity stack */
4 #include "cddefines.h"
5 #include "taulines.h"
6 #include "iso.h"
7 #include "geometry.h"
8 #include "heavy.h"
9 #include "dense.h"
10 #include "prt.h"
11 #include "opacity.h"
12 #include "coolheavy.h"
13 #include "phycon.h"
14 #include "rfield.h"
15 #include "predcont.h"
16 #include "radius.h"
17 #include "continuum.h"
18 #include "lines.h"
19 #include "freebound.h"
20 #include "lines_service.h"
21 
22 void lines_continuum(void)
23 {
24 
25  double f1,
26  f2 ,
27  bac ,
28  flow;
29  long i,nBand;
30 
31  DEBUG_ENTRY( "lines_continuum()" );
32 
33  /* code has all local emissivities zeroed out with cryptic comment about being
34  * situation dependent. Why? this is option to turn back on */
35  const bool KILL_CONT = false;
36 
37  i = StuffComment( "continua" );
38  linadd( 0., (realnum)i , "####", 'i',
39  " start continua");
40 
41  /* these entries only work correctly if the APERTURE command is not in effect */
42  if( geometry.iEmissPower == 2 )
43  {
44  /***********************************************************************
45  * stuff in Bac ratio - continuum above the Balmer Jump
46  * this is trick, zeroing out saved continuum integrated so far,
47  * and adding the current version, so that the line array gives the
48  * value in the final continuum
49  *
50  * reflected continuum is different from others since relative to inner
51  * radius, others for for this radius
52  *************************************************************************/
53 
55  /***************************************************************************
56  * "Bac " , 3646, this is residual continuum at peak of Balmer Jump
57  * flux below - flux above
58  ***************************************************************************/
59  /* >>chng 00 dec 02, remove opac.tmn */
60  /* >>chng 00 dec 19, remove / radius.GeoDil */
61  /* extrapolated continuum above head */
62  /* >>chng 01 jul 13, from ConInterOut to ConEmitOut */
63  f1 = (rfield.ConEmitOut[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1] +
64  rfield.ConEmitReflec[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1]/radius.r1r0sq )/
65  rfield.widflx(iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1);
66 
67  /* extrapolated continuum below head */
68  /* >>chng 00 dec 19, remove / radius.GeoDil */
69  f2 = (rfield.ConEmitOut[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-2]+
70  rfield.ConEmitReflec[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-2]/radius.r1r0sq )/
71  rfield.widflx(iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-2);
72 
73  /* convert to nuFnu units */
74  f1 = f1*0.250*0.250*EN1RYD*radius.r1r0sq;
75  f2 = f2*0.250*0.250*EN1RYD*radius.r1r0sq;
76  bac = (f1 - f2);
77 
78  /* memory not allocated until ipass >= 0
79  * clear summed intrinsic and emergent intensity of following
80  * entry - following call to linadd will enter the total and
81  * keep entering the total but is done for each zone hence need to
82  * keep resetting to zero*/
83  if( LineSave.ipass > 0 )
84  {
85  LineSave.lines[LineSave.nsum].SumLineZero();
86  }
87 
88  linadd(MAX2(0.,bac)/radius.dVeffAper,3646,"Bac ",'i',
89  "residual flux at head of Balmer continuum, nuFnu ");
90  /* >>chng 03 feb 06, set to zero */
91  /* emslin saves the per unit vol emissivity of a line, which is normally
92  * what goes into linadd. We zero this unit emissivity which was set
93  * FOR THE PREVIOUS LINE since it is so situation dependent */
94  if( KILL_CONT && LineSave.ipass > 0 )
95  {
96  LineSave.lines[LineSave.nsum-1].emslinZero();
97  }
98 
99  /* memory not allocated until ipass >= 0 */
100  if( LineSave.ipass > 0 )
101  {
102  LineSave.lines[LineSave.nsum].SumLineZero();
103  }
104 
105  linadd(f1/radius.dVeffAper,3645,"nFnu",'i',
106  "total flux above head of Balmer continuum, nuFnu ");
107  /* >>chng 03 feb 06, set to zero */
108  /* emslin saves the per unit vol emissivity of a line, which is normally
109  * what goes into linadd. We zero this unit emissivity which was set
110  * FOR THE PREVIOUS LINE since it is so situation dependent */
111  if( KILL_CONT && LineSave.ipass > 0 )
112  {
113  LineSave.lines[LineSave.nsum-1].emslinZero();
114  }
115 
116  /* memory not allocated until ipass >= 0 */
117  if( LineSave.ipass > 0 )
118  {
119  LineSave.lines[LineSave.nsum].SumLineZero();
120  }
121 
122  linadd(f2/radius.dVeffAper,3647,"nFnu",'i',
123  "total flux above head of Balmer continuum, nuFnu ");
124  /* >>chng 03 feb 06, set to zero */
125  /* emslin saves the per unit vol emissivity of a line, which is normally
126  * what goes into linadd. We zero this unit emissivity which was set
127  * FOR THE PREVIOUS LINE since it is so situation dependent */
128  if( KILL_CONT && LineSave.ipass > 0 )
129  {
130  LineSave.lines[LineSave.nsum-1].emslinZero();
131  }
132 
133  /******************************************************************************
134  * "cout" , 3646, this is outward residual continuum at peak of Balmer Jump *
135  * equal to total in spherical geometry, half in opt thin open geometry *
136  ******************************************************************************/
137  /* >>chng 00 dec 02, remove opac.tmn */
138  /* >>chng 00 dec 19, remove / radius.GeoDil */
139  f1 = rfield.ConEmitOut[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1]/
140  rfield.widflx(iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1);
141 
142  /* >>chng 00 dec 19, remove / radius.GeoDil */
143  f2 = rfield.ConEmitOut[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-2]/
144  rfield.widflx(iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-2);
145 
146  /* net Balmer jump */
147  bac = (f1 - f2)*0.250*0.250*EN1RYD*radius.r1r0sq;
148 
149  /* memory not allocated until ipass >= 0 */
150  if( LineSave.ipass > 0 )
151  {
152  LineSave.lines[LineSave.nsum].SumLineZero();
153  }
154 
155  linadd(MAX2(0.,bac)/radius.dVeffAper,3646,"cout",'i',
156  "residual flux in Balmer continuum, nuFnu ");
157  /* >>chng 03 feb 06, set to zero */
158  /* emslin saves the per unit vol emissivity of a line, which is normally
159  * what goes into linadd. We zero this unit emissivity which was set
160  * FOR THE PREVIOUS LINE since it is so situation dependent */
161  if( KILL_CONT && LineSave.ipass > 0 )
162  {
163  LineSave.lines[LineSave.nsum-1].emslinZero();
164  }
165 
166  /*********************************************************************
167  * "cref" , 3646, this is reflected continuum at peak of Balmer Jump*
168  * equal to zero in spherical geometry, half of total in op thin opn *
169  *********************************************************************/
170  /* >>chng 00 dec 02, remove opac.tmn */
171  /* >>chng 00 dec 19, remove / radius.GeoDil */
172  f1 = rfield.ConEmitReflec[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1]/
173  rfield.widflx(iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1);
174 
175  f2 = rfield.ConEmitReflec[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-2]/
176  rfield.widflx(iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-2);
177 
178  /* net Balmer jump */
179  bac = (f1 - f2)*0.250*0.250*EN1RYD;
180 
181  /* memory not allocated until ipass >= 0 */
182  if( LineSave.ipass > 0 )
183  {
184  LineSave.lines[LineSave.nsum].SumLineZero();
185  }
186 
187  linadd(MAX2(0.,bac)/radius.dVeffAper,3646,"cref",'i',
188  "residual flux in Balmer continuum, nuFnu ");
189  /* >>chng 03 feb 06, set to zero */
190  /* emslin saves the per unit vol emissivity of a line, which is normally
191  * what goes into linadd. We zero this unit emissivity which was set
192  * FOR THE PREVIOUS LINE since it is so situation dependent */
193  if( KILL_CONT && LineSave.ipass > 0 )
194  {
195  LineSave.lines[LineSave.nsum-1].emslinZero();
196  }
197 
198  /*********************************************************************
199  * "thin" , 3646, tot optically thin continuum at peak of Balmer Jump*/
200  if( nzone > 0 )
201  {
202  /* rfield.ConEmitLocal is not defined initially, only evaluate when into model */
203  f1 = rfield.ConEmitLocal[nzone][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1]/
204  rfield.widflx(iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-1);
205  f2 = rfield.ConEmitLocal[nzone][iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-2]/
206  rfield.widflx(iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].ipIsoLevNIonCon-2);
207  }
208  else
209  {
210  f1 = 0.;
211  f2 = 0.;
212  }
213 
214  bac = (f1 - f2)*0.250*0.250*EN1RYD;
215 
216  linadd(MAX2(0.,bac),3646,"thin",'i',
217  "residual flux in Balmer continuum, nuFnu ");
218 
219  linadd(continuum.cn4861/radius.dVeffAper,4860,"Inci",'i',
220  "incident continuum nu*f_nu at H-beta, at illuminated face of cloud ");
221 
222  linadd(continuum.cn1216/radius.dVeffAper,1215,"Inci",'i',
223  "incident continuum nu*f_nu near Ly-alpha, at illuminated face of cloud");
224 
225  if( LineSave.ipass > 0 )
226  {
227  continuum.cn4861 = 0.;
228  continuum.cn1216 = 0.;
229  }
230  }
231 
232  flow = (iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2p].RadRecomb[ipRecRad] +
233  iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2s].RadRecomb[ipRecRad])*
234  iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2p].RadRecomb[ipRecEsc]*
235  dense.eden*dense.xIonDense[ipHYDROGEN][1]* 5.45e-12;
236  linadd(flow,0,"Ba C",'i',
237  "integrated Balmer continuum emission");
238 
239  if( iso_sp[ipH_LIKE][ipHYDROGEN].n_HighestResolved_max >= 3 )
240  {
241  flow = ( iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH3s].RadRecomb[ipRecRad]*
242  iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH3s].RadRecomb[ipRecEsc] +
243  iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH3p].RadRecomb[ipRecRad]*
244  iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH3p].RadRecomb[ipRecEsc] +
245  iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH3d].RadRecomb[ipRecRad]*
246  iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH3d].RadRecomb[ipRecEsc] ) *
247  dense.eden*dense.xIonDense[ipHYDROGEN][1]*3.53e-12;
248  }
249  else
250  {
251  flow = iso_sp[ipH_LIKE][ipHYDROGEN].fb[3].RadRecomb[ipRecRad]*
252  iso_sp[ipH_LIKE][ipHYDROGEN].fb[3].RadRecomb[ipRecEsc]*
253  dense.eden*dense.xIonDense[ipHYDROGEN][1]*3.53e-12;
254  }
255  linadd(flow,0,"PA C",'i',
256  "Paschen continuum emission ");
257 
258  /* these are a series of continuum bands defined in the file
259  * continuum_bands.ini - this makes it possible to enter any
260  * integrated total emission into the emission-line stack */
261  /* these entries only work correctly if the APERTURE command is not in effect */
262  if( geometry.iEmissPower == 2 )
263  {
264  for( nBand=0; nBand < continuum.nContBand; ++nBand )
265  {
266  double EmergentContinuum = 0.;
267  double DiffuseEmission = 0.;
268  if( LineSave.ipass > 0 )
269  {
270  /* find total emission over band - units will be erg cm-2 s-1 */
271  for( i=continuum.ipContBandLow[nBand]; i<=continuum.ipContBandHi[nBand]; ++i )
272  {
273  // correction for fraction of low or hi cell
274  // that lies within the band
275  double EdgeCorrection = 1.;
276  if( i==continuum.ipContBandLow[nBand] )
277  EdgeCorrection = continuum.BandEdgeCorrLow[nBand];
278  else if( i==continuum.ipContBandHi[nBand])
279  EdgeCorrection = continuum.BandEdgeCorrHi[nBand];
280 
281  double xIntenOut =
282  /* the attenuated incident continuum */
283  flux_correct_isotropic( 0, i-1 ) +
284 
285  // the outward emitted continuous radiation field
286  (rfield.ConEmitOut[0][i-1] +
287 
288  /* outward emitted lines */
289  rfield.outlin[0][i-1])*geometry.covgeo;
290  xIntenOut *= EdgeCorrection;
291 
292  /* div by opac.E2TauAbsFace[i] because ConEmitReflec has already
293  * been corrected for this - emergent_line will introduce a
294  * further correction so this will cancel the extra factor */
295  /* NB: Comparison to 1e-37 suppresses overflows when realnum
296  * is double (FLT_IS_DBL). */
297  double xIntenIn = 0.;
298  if( opac.E2TauAbsFace[i-1] > 1e-37 )
299  xIntenIn = (double)rfield.ConEmitReflec[0][i-1]/
300  (double)opac.E2TauAbsFace[i-1]*geometry.covgeo;
301  /* outward emitted lines */
302  xIntenIn += rfield.reflin[0][i-1]*geometry.covgeo;
303  xIntenIn *= EdgeCorrection;
304 
305  /* the fraction of this that gets out */
306  EmergentContinuum += rfield.anu(i-1) *
307  emergent_line( xIntenIn , xIntenOut , i )
308  / SDIV(opac.tmn[i-1]);
309 
310  // diffuse local emission
311  DiffuseEmission += (rfield.ConEmitLocal[nzone][i-1] +
313  EdgeCorrection;
314  }
315  }
316  /* we will call lindst with an energy half way between the two
317  * ends of the band. This will make an extinction correction that
318  * has already been applied above by multiplying by emergent_line.
319  * Find this factor and remove it before the call */
320  double corr = emergent_line( 0.5 , 0.5 ,
321  (continuum.ipContBandLow[nBand]+continuum.ipContBandHi[nBand])/2 );
322  /* NB: Comparison to 1e-37 suppresses overflows when realnum
323  * is double (FLT_IS_DBL). */
324  if( corr < 1e-37 )
325  EmergentContinuum = 0.;
326  else
327  EmergentContinuum /= corr;
328 
329  /* convert to physical units */
330  EmergentContinuum *= EN1RYD*radius.r1r0sq/radius.dVeffAper;
331  DiffuseEmission *= EN1RYD;
332  /* memory not allocated until ipass >= 0 */
333  if( LineSave.ipass > 0 )
334  {
335  LineSave.lines[LineSave.nsum].SumLineZero();
336  }
337 
338  lindst( EmergentContinuum ,
339  // the negative wavelength is a sentinel that the wavelength
340  // may be bogus - very often the "center" of the band, as defined
341  // by observers, is far to the blue end of the range. use the
342  // observed definition of the wavelength. This introduces an error
343  // since the wavelength is used to determine the transfer of the
344  // band against background opacities
345  -continuum.ContBandWavelength[nBand] ,
347  (continuum.ipContBandLow[nBand]+continuum.ipContBandHi[nBand])/2 ,
348  't' , false,
349  "continuum bands defined in continuum_bands.ini");
350 
351  // emissivity has no meaning for these bands - quantity is net
352  // transmitted radiation field
353  if( LineSave.ipass > 0 )
354  {
355  LineSave.lines[LineSave.nsum-1].emslinSet(0,DiffuseEmission);
356  LineSave.lines[LineSave.nsum-1].emslinThin();
357  }
358  }
359  }
360 
361  linadd(MAX2(0.,CoolHeavy.brems_cool_net),0,"HFFc",'c',
362  "net free-free cooling, ALL species, free-free heating subtracted, so nearly cancels with cooling in LTE ");
363 
364  linadd(MAX2(0.,-CoolHeavy.brems_cool_net),0,"HFFh",'h',
365  "net free-free heating, nearly cancels with cooling in LTE ");
366 
367  linadd(CoolHeavy.brems_cool_h,0,"H FF",'i',
368  " H brems (free-free) cooling ");
369 
370  linadd(CoolHeavy.brems_heat_total,0,"FF H",'i',
371  "total free-free heating ");
372 
373  linadd(CoolHeavy.brems_cool_he,0,"HeFF",'i',
374  "He brems emission ");
375 
376  linadd(CoolHeavy.heavfb,0,"MeFB",'c',
377  "heavy element recombination cooling ");
378 
379  linadd(CoolHeavy.brems_cool_metals,0,"MeFF",'i',
380  "heavy elements (metals) brems cooling, heat not subtracted ");
381 
383  "total brems emission - total cooling but not minus heating ");
384 
386  "part of H brems, in x-ray beyond 0.5KeV ");
387 
388  linadd(CoolHeavy.eebrm,0,"eeff",'c',
389  "electron - electron brems ");
390 
391  linadd(CoolHeavy.colmet,0,"Mion",'c',
392  " cooling due to collisional ionization of heavy elements" );
393 
394  /* predict emitted continuum at series of continuum points */
395  /* class is located in predcont.h,
396  * PredCont - contains vector of pair of Energy and ip where
397  * we want to predict the continuum,
398  *
399  * the entry nFnu will only be printed if the command
400  * print diffuse continuum
401  * is entered -
402  *
403  * this code should be kept parallel with that in dopunch, where
404  * save continuum is produced, since two must agree */
405 
406  t_PredCont& PredCont = t_PredCont::Inst();
407  if( LineSave.ipass == 0 )
408  PredCont.set_offset(LineSave.nsum);
409 
410  /* these entries only work correctly if the APERTURE command is not in effect */
411  if( geometry.iEmissPower == 2 )
412  {
413  for( i=0; i < long(PredCont.size()); i++ )
414  {
415  double SourceTransmitted , Cont_nInu;
416  double SourceReflected, DiffuseOutward, DiffuseInward;
417  double renorm;
418 
419  /* put wavelength in Angstroms into dummy structure, so that we can use iWavLen
420  * to get a proper wavelength with units, continuum energies are stored in PredCont */
421  (*TauDummy).WLAng() = (realnum)PredCont[i].Angstrom();
422  /*lambda = iWavLen(TauDummy , &chUnits , &chShift );*/
423 
424  /* >>chng 00 dec 02, there were three occurrences of /opac.tmn which had the
425  * effect of raising the summed continuum by the local opacity correction factor.
426  * in the case of the Lyman continuum this raised the reported value by orders
427  * of magnitude. There have been commented out in the following for now. */
428  /* reflected total continuum (diff+incident emitted inward direction) */
429 
430  /* >>chng 00 dec 08, implement the "set nFnu [SOURCE_REFLECTED] ... command, PvH */
431  /* >>chng 00 dec 19, remove / radius.GeoDil */
432  renorm = rfield.anu2(PredCont[i].ip_C())*EN1RYD/rfield.widflx(PredCont[i].ip_C());
433 
434  /* this is the reflected diffuse continuum */
435  if( prt.lgDiffuseInward )
436  {
437  DiffuseInward = rfield.ConEmitReflec[0][PredCont[i].ip_C()]*renorm;
438  }
439  else
440  {
441  DiffuseInward = 0.;
442  }
443 
444  /* the outward diffuse continuum */
445  if( prt.lgDiffuseOutward )
446  {
447  DiffuseOutward = rfield.ConEmitOut[0][PredCont[i].ip_C()]*renorm*radius.r1r0sq;
448  }
449  else
450  {
451  DiffuseOutward = 0.;
452  }
453 
454  /* reflected part of INCIDENT continuum (only incident, not diffuse, which was above) */
455  if( prt.lgSourceReflected )
456  {
457  SourceReflected = rfield.ConRefIncid[0][PredCont[i].ip_C()]*renorm;
458  }
459  else
460  {
461  SourceReflected = 0.;
462  }
463 
464  /* the attenuated incident continuum */
466  {
467  SourceTransmitted = rfield.flux[0][PredCont[i].ip_C()]*renorm*radius.r1r0sq;
468  }
469  else
470  {
471  SourceTransmitted = 0.;
472  }
473 
474  /* memory has not been allocated until ipass >= 0, so must not access this element,
475  * this element will be used to save the following quantity */
476  if( LineSave.ipass > 0 )
477  {
478  LineSave.lines[LineSave.nsum].SumLineZero();
479  }
480 
481  linadd((DiffuseInward+SourceReflected+DiffuseOutward+SourceTransmitted)/radius.dVeffAper,
482  (*TauDummy).WLAng(),"nFnu",'i',
483  "total continuum at selected energy points " );
484 
485  /* emslin saves the per unit vol emissivity of a line, which is normally
486  * what goes into linadd. We zero this unit emissivity which was set
487  * FOR THE PREVIOUS LINE since it is so situation dependent */
488  if( KILL_CONT && LineSave.ipass > 0 )
489  {
490  LineSave.lines[LineSave.nsum-1].emslinZero();
491  }
492 
493  /* this is the normal set to zero to trick the NEXT line into going in properly */
494  if( LineSave.ipass > 0 )
495  {
496  LineSave.lines[LineSave.nsum].SumLineZero();
497  }
498 
499  /* the nsum-1 -- emslin and nsum -- SumLine is not a bug, look above - they do
500  * different things to different saves */
501  Cont_nInu = rfield.flux[0][PredCont[i].ip_C()]*renorm*radius.r1r0sq +
502  rfield.ConRefIncid[0][PredCont[i].ip_C()]*renorm;
503 
504 # if 0
505  /* this code can be used to create assert statements for the continuum shape */
506  if( !i )
507  fprintf(ioQQQ,"\n");
508  char chWL[1000];
509  sprt_wl( chWL , (*TauDummy).WLAng() );
510  fprintf( ioQQQ,"assert line luminosity \"nInu\" %s %.3f\n",
511  chWL ,
512  log10(SDIV(Cont_nInu/radius.dVeffAper) * radius.Conv2PrtInten) );
513 # endif
514 
515  linadd( Cont_nInu/radius.dVeffAper,(*TauDummy).WLAng(),"nInu",'i',
516  "transmitted and reflected incident continuum at selected energy points " );
517 
518  /* emslin saves the per unit volume emissivity of a line, which is normally
519  * what goes into linadd. We zero this unit emissivity since it is so situation dependent */
520  if( KILL_CONT && LineSave.ipass > 0 )
521  {
522  LineSave.lines[LineSave.nsum-1].emslinZero();
523  }
524 
525  /* memory has not been allocated until ipass >= 0 */
526  if( LineSave.ipass > 0 )
527  {
528  LineSave.lines[LineSave.nsum].SumLineZero();
529  }
530 
531  linadd( (DiffuseInward+SourceReflected)/radius.dVeffAper,(*TauDummy).WLAng(),"InwT",'i',
532  "total reflected continuum, total inward emission plus reflected (diffuse) total continuum ");
533 
534  if( KILL_CONT && LineSave.ipass > 0 )
535  {
536  LineSave.lines[LineSave.nsum-1].emslinZero();
537  }
538 
539  /* memory has not been allocated until ipass >= 0 */
540  if( LineSave.ipass > 0 )
541  {
542  LineSave.lines[LineSave.nsum].SumLineZero();
543  }
544 
545  linadd(SourceReflected/radius.dVeffAper,(*TauDummy).WLAng(),"InwC",'i',
546  "reflected incident continuum (only incident) ");
547 
548  if( KILL_CONT && LineSave.ipass > 0 )
549  {
550  LineSave.lines[LineSave.nsum-1].emslinZero();
551  }
552  }
553  }
554 
555  i = StuffComment( "RRC" );
556  linadd( 0., (realnum)i , "####", 'i',"radiative recombination continua");
557 
558  // radiative recombination continua, RRC, for iso sequences
559  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
560  {
561  for( long nelem=ipISO; nelem < LIMELM; nelem++ )
562  {
563  if( nelem < 2 || dense.lgElmtOn[nelem] )
564  {
565  for( long n=0; n < iso_sp[ipISO][nelem].numLevels_max; n++ )
566  {
567  if( LineSave.ipass < 0 )
568  // this pass only counting lines
569  linadd(0.,0.,"dumy",'i',"radiative recombination continuum");
570  else if( LineSave.ipass == 0 )
571  {
572  // save wavelength and label
573  /* chIonLbl generates a null terminated 4 char string, of form "C 2"
574  * the result, chLable, is only used when ipass == 0, can be undefined otherwise */
575  realnum wl = (realnum)(RYDLAM / iso_sp[ipISO][nelem].fb[n].xIsoLevNIonRyd);
576  wl /= (realnum)RefIndex( 1e8/wl );
577  linadd( 0. , wl ,chIonLbl(iso_sp[ipISO][nelem].trans(1,0)).c_str(),'i',
578  "radiative recombination continuum");
579  }
580  else
581  {
582  // save intensity
583  linadd(iso_sp[ipISO][nelem].fb[n].RadRecCon,0,"dumy",'i',
584  "radiative recombination continuum");
585  }
586  }
587  }
588  }
589  }
590 
591  // RRC for non iso sequence ions
592 
593  /* add recombination continua for elements heavier than those done with iso seq */
594  for( long nelem=NISO; nelem < LIMELM; nelem++ )
595  {
596  /* do not include species with iso-sequence in following */
597  /* >>chng 03 sep 09, upper bound was wrong, did not include NISO */
598  for( long ion=0; ion < nelem-NISO+1; ion++ )
599  {
600  if( dense.lgElmtOn[nelem] )
601  {
602  if( LineSave.ipass < 0 )
603  // this pass only counting lines
604  linadd(0.,0.,"dumy",'i',"radiative recombination continuum");
605  else if( LineSave.ipass == 0 )
606  {
607  string chLabel = chIonLbl( nelem+1, ion+1 );
608  // the constant 0.9998787 has been removed elsewhere in the code but has been
609  // added here to prevent wavelengths from changing between c17.01 and later c17
610  // versions; this could break existing scripts; see also ticket #407
611  realnum wl = (realnum)(RYDLAM / (Heavy.Valence_IP_Ryd[nelem][ion]* 0.9998787));
612  wl /= (realnum)RefIndex( 1e8/wl );
613  linadd( 0. , wl ,chLabel.c_str(),'i',
614  "radiative recombination continuum");
615  }
616  else
617  {
618  // save intensity
619  linadd(Heavy.RadRecCon[nelem][ion],0,"dumy",'i',
620  "radiative recombination continuum");
621  }
622  }
623  }
624  }
625 
626  return;
627 }
double RadRecCon[LIMELM][LIMELM]
Definition: heavy.h:18
string chIonLbl(const TransitionProxy &t)
Definition: transition.cpp:233
double widflx(size_t i) const
Definition: mesh.h:156
t_opac opac
Definition: opacity.cpp:5
realnum ** flux
Definition: rfield.h:68
realnum * DiffuseLineEmission
Definition: rfield.h:193
t_Heavy Heavy
Definition: heavy.cpp:5
const int NISO
Definition: cddefines.h:311
long int * ipContBandLow
Definition: continuum.h:88
void lindst(double xEmiss, realnum wavelength, const char *chLab, long int ipnt, char chInfo, bool lgOutToo, const char *chComment)
realnum * BandEdgeCorrHi
Definition: continuum.h:91
long int iEmissPower
Definition: geometry.h:71
t_phycon phycon
Definition: phycon.cpp:6
double brems_cool_he
Definition: coolheavy.h:36
t_LineSave LineSave
Definition: lines.cpp:10
bool lgDiffuseInward
Definition: prt.h:209
realnum cn4861
Definition: continuum.h:75
sys_float sexp(sys_float x)
Definition: service.cpp:999
double RefIndex(double EnergyWN)
t_CoolHeavy CoolHeavy
Definition: coolheavy.cpp:5
realnum covgeo
Definition: geometry.h:45
realnum ** outlin
Definition: rfield.h:189
FILE * ioQQQ
Definition: cddefines.cpp:7
long int nzone
Definition: cddefines.cpp:14
double eebrm
Definition: coolheavy.h:18
vector< freeBound > fb
Definition: iso.h:481
void lines_continuum(void)
double anu(size_t i) const
Definition: mesh.h:120
vector< LinSv > lines
Definition: lines.h:132
void set_offset(long offset)
Definition: predcont.h:32
double brems_cool_net
Definition: coolheavy.h:36
t_dense dense
Definition: global.cpp:15
static T & Inst()
Definition: cddefines.h:209
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
bool lgSourceReflected
Definition: prt.h:207
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
realnum ** ConEmitLocal
Definition: rfield.h:139
t_geometry geometry
Definition: geometry.cpp:5
bool lgDiffuseOutward
Definition: prt.h:210
LinSv * linadd(double xEmiss, realnum wavelength, const char *chLab, char chInfo, const char *chComment)
t_continuum continuum
Definition: continuum.cpp:6
double brems_cool_h
Definition: coolheavy.h:36
t_rfield rfield
Definition: rfield.cpp:9
float realnum
Definition: cddefines.h:124
realnum * ContBandWavelength
Definition: continuum.h:87
const int ipRecRad
Definition: cddefines.h:333
double brems_heat_total
Definition: coolheavy.h:36
double anu2(size_t i) const
Definition: mesh.h:124
const int ipRecEsc
Definition: cddefines.h:329
const int ipH3s
Definition: iso.h:32
double heavfb
Definition: coolheavy.h:18
const int ipH3d
Definition: iso.h:34
t_radius radius
Definition: radius.cpp:5
long int nContBand
Definition: continuum.h:85
t_prt prt
Definition: prt.cpp:14
realnum ** reflin
Definition: rfield.h:196
realnum ** ConEmitOut
Definition: rfield.h:151
bool lgElmtOn[LIMELM]
Definition: dense.h:160
size_t size() const
Definition: predcont.h:28
double Conv2PrtInten
Definition: radius.h:153
double brems_cool_metals
Definition: coolheavy.h:36
const int ipH2p
Definition: iso.h:31
void sprt_wl(char *chString, realnum wl)
Definition: prt.cpp:56
const int ipH2s
Definition: iso.h:30
const int ipH_LIKE
Definition: iso.h:64
const int LIMELM
Definition: cddefines.h:308
realnum * BandEdgeCorrLow
Definition: continuum.h:91
double Valence_IP_Ryd[LIMELM][LIMELM]
Definition: heavy.h:24
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
realnum * E2TauAbsFace
Definition: opacity.h:137
double eden
Definition: dense.h:201
#define MAX2(a, b)
Definition: cddefines.h:824
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
double emergent_line(double emissivity_in, double emissivity_out, long int ipCont)
bool lgSourceTransmitted
Definition: prt.h:208
sys_float SDIV(sys_float x)
Definition: cddefines.h:1002
long int numLevels_max
Definition: iso.h:524
long int nsum
Definition: lines.h:87
realnum * tmn
Definition: opacity.h:149
realnum ** ConEmitReflec
Definition: rfield.h:145
double r1r0sq
Definition: radius.h:31
long int * ipContBandHi
Definition: continuum.h:88
double te
Definition: phycon.h:21
const int ipHYDROGEN
Definition: cddefines.h:349
double flux_correct_isotropic(const bool lgSaveIsotr, const int nEmType, const int iflux)
Definition: rfield.cpp:112
realnum ** ConRefIncid
Definition: rfield.h:157
const int ipH3p
Definition: iso.h:33
long int ipass
Definition: lines.h:96
double dVeffAper
Definition: radius.h:93
realnum cn1216
Definition: continuum.h:75
char ** chContBandLabels
Definition: continuum.h:86
long int StuffComment(const char *chComment)
Definition: prt_final.cpp:1938
double colmet
Definition: coolheavy.h:18