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