cloudy  trunk
 All Data Structures 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-2008 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 "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 "lines_service.h"
17 #include "radius.h"
18 #include "continuum.h"
19 #include "lines.h"
20 
21 void lines_continuum(void)
22 {
23 
24  double f1,
25  f2 ,
26  bac ,
27  flow;
28  long i,nBand;
29 
30  DEBUG_ENTRY( "lines_continuum()" );
31 
32  /* code has all local emissivities zeroed out with cryptic comment about being
33  * situation dependent. Why? this is option to turn back on */
34  const bool KILL_CONT = false;
35 
36  i = StuffComment( "continua" );
37  linadd( 0., (realnum)i , "####", 'i',
38  " start continua");
39 
40  /***********************************************************************
41  * stuff in Bac ratio - continuum above the Balmer Jump
42  * this is trick, zeroing out saved continuum integrated so far,
43  * and adding the current version, so that the line array gives the
44  * value in the final continuum
45  *
46  * reflected continuum is different from others since relative to inner
47  * radius, others for for this radius
48  *************************************************************************/
49 
51  /***************************************************************************
52  * "Bac " , 3646, this is residual continuum at peak of Balmer Jump
53  * flux below - flux above
54  ***************************************************************************/
55  /* >>chng 00 dec 02, remove opac.tmn */
56  /* >>chng 00 dec 19, remove / radius.GeoDil */
57  /* extrapolated continuum above head */
58  /* >>chng 01 jul 13, from ConInterOut to ConEmitOut */
61  rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]/*/radius.GeoDil /opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]*/;
62 
63  /* extrapolated continuum below head */
64  /* >>chng 00 dec 19, remove / radius.GeoDil */
67  rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]/*/radius.GeoDil /opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]*/;
68 
69  /* convert to nuFnu units */
70  f1 = f1*0.250*0.250*EN1RYD*radius.r1r0sq;
71  f2 = f2*0.250*0.250*EN1RYD*radius.r1r0sq;
72  bac = (f1 - f2);
73 
74  /* memory not allocated until ipass >= 0
75  * clear summed intrinsic and emergent intensity of following
76  * entry - following call to linadd will enter the total and
77  * keep entering the total but is done for each zone hence need to
78  * keep resetting to zero*/
79  if( LineSave.ipass > 0 )
80  {
81  LineSv[LineSave.nsum].sumlin[0] = 0.;
82  LineSv[LineSave.nsum].sumlin[1] = 0.;
83  }
84 
85  linadd(MAX2(0.,bac)/radius.dVeff,3646,"Bac ",'i',
86  "residual flux at head of Balmer continuum, nuFnu ");
87  /* >>chng 03 feb 06, set to zero */
88  /* emslin saves the per unit vol emissivity of a line, which is normally
89  * what goes into linadd. We zero this unit emissivity which was set
90  * FOR THE PREVIOUS LINE since it is so situation dependent */
91  if( KILL_CONT && LineSave.ipass > 0 )
92  {
93  LineSv[LineSave.nsum-1].emslin[0] = 0.;
94  LineSv[LineSave.nsum-1].emslin[1] = 0.;
95  }
96 
97  /* memory not allocated until ipass >= 0 */
98  if( LineSave.ipass > 0 )
99  {
100  LineSv[LineSave.nsum].sumlin[0] = 0.;
101  LineSv[LineSave.nsum].sumlin[1] = 0.;
102  }
103 
104  linadd(f1/radius.dVeff,3645,"nFnu",'i',
105  "total flux above head of Balmer continuum, nuFnu ");
106  /* >>chng 03 feb 06, set to zero */
107  /* emslin saves the per unit vol emissivity of a line, which is normally
108  * what goes into linadd. We zero this unit emissivity which was set
109  * FOR THE PREVIOUS LINE since it is so situation dependent */
110  if( KILL_CONT && LineSave.ipass > 0 )
111  {
112  LineSv[LineSave.nsum-1].emslin[0] = 0.;
113  LineSv[LineSave.nsum-1].emslin[1] = 0.;
114  }
115 
116  /* memory not allocated until ipass >= 0 */
117  if( LineSave.ipass > 0 )
118  {
119  LineSv[LineSave.nsum].sumlin[0] = 0.;
120  LineSv[LineSave.nsum].sumlin[1] = 0.;
121  }
122 
123  linadd(f2/radius.dVeff,3647,"nFnu",'i',
124  "total flux above head of Balmer continuum, nuFnu ");
125  /* >>chng 03 feb 06, set to zero */
126  /* emslin saves the per unit vol emissivity of a line, which is normally
127  * what goes into linadd. We zero this unit emissivity which was set
128  * FOR THE PREVIOUS LINE since it is so situation dependent */
129  if( KILL_CONT && LineSave.ipass > 0 )
130  {
131  LineSv[LineSave.nsum-1].emslin[0] = 0.;
132  LineSv[LineSave.nsum-1].emslin[1] = 0.;
133  }
134 
135  /******************************************************************************
136  * "cout" , 3646, this is outward residual continuum at peak of Balmer Jump *
137  * equal to total in spherical geometry, half in opt thin open geometry *
138  ******************************************************************************/
139  /* >>chng 00 dec 02, remove opac.tmn */
140  /* >>chng 00 dec 19, remove / radius.GeoDil */
142  rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]/*/radius.GeoDil /opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]*/;
143 
144  /* >>chng 00 dec 19, remove / radius.GeoDil */
146  rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]/*/radius.GeoDil /opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]*/;
147 
148  /* net Balmer jump */
149  bac = (f1 - f2)*0.250*0.250*EN1RYD*radius.r1r0sq;
150 
151  /* memory not allocated until ipass >= 0 */
152  if( LineSave.ipass > 0 )
153  {
154  LineSv[LineSave.nsum].sumlin[0] = 0.;
155  LineSv[LineSave.nsum].sumlin[1] = 0.;
156  }
157 
158  linadd(MAX2(0.,bac)/radius.dVeff,3646,"cout",'i',
159  "residual flux in Balmer continuum, nuFnu ");
160  /* >>chng 03 feb 06, set to zero */
161  /* emslin saves the per unit vol emissivity of a line, which is normally
162  * what goes into linadd. We zero this unit emissivity which was set
163  * FOR THE PREVIOUS LINE since it is so situation dependent */
164  if( KILL_CONT && LineSave.ipass > 0 )
165  {
166  LineSv[LineSave.nsum-1].emslin[0] = 0.;
167  LineSv[LineSave.nsum-1].emslin[1] = 0.;
168  }
169 
170  /*********************************************************************
171  * "cref" , 3646, this is reflected continuum at peak of Balmer Jump*
172  * equal to zero in spherical geometry, half of total in op thin opn *
173  *********************************************************************/
174  /* >>chng 00 dec 02, remove opac.tmn */
175  /* >>chng 00 dec 19, remove / radius.GeoDil */
178 
181 
182  /* net Balmer jump */
183  bac = (f1 - f2)*0.250*0.250*EN1RYD;
184 
185  /* memory not allocated until ipass >= 0 */
186  if( LineSave.ipass > 0 )
187  {
188  LineSv[LineSave.nsum].sumlin[0] = 0.;
189  LineSv[LineSave.nsum].sumlin[1] = 0.;
190  }
191 
192  linadd(MAX2(0.,bac)/radius.dVeff,3646,"cref",'i',
193  "residual flux in Balmer continuum, nuFnu ");
194  /* >>chng 03 feb 06, set to zero */
195  /* emslin saves the per unit vol emissivity of a line, which is normally
196  * what goes into linadd. We zero this unit emissivity which was set
197  * FOR THE PREVIOUS LINE since it is so situation dependent */
198  if( KILL_CONT && LineSave.ipass > 0 )
199  {
200  LineSv[LineSave.nsum-1].emslin[0] = 0.;
201  LineSv[LineSave.nsum-1].emslin[1] = 0.;
202  }
203 
204  /*********************************************************************
205  * "thin" , 3646, tot optically thin continuum at peak of Balmer Jump*/
206  if( nzone > 0 )
207  {
208  /* rfield.ConEmitLocal is not defined initially, only evaluate when into model */
213  bac = (f1 - f2)*0.250*0.250*EN1RYD;
214  }
215  else
216  {
217  f1 = 0.;
218  f2 = 0.;
219  bac = 0.;
220  }
221 
222  linadd(MAX2(0.,bac),3646,"thin",'i',
223  "residual flux in Balmer continuum, nuFnu ");
224 
226  linadd(continuum.cn4861,4860,"Inci",'i',
227  "incident continuum nu*f_nu at H-beta, at illuminated face of cloud ");
228  /* >>chng 07 jun 13, at this point nsum is the number of lines in the stack
229  * so nsum-1 is the entry created by the above call to linadd
230  * we want to set the emergent incident continuum to the intrinsic
231  * incident continuum - there are reports of the incident continuum
232  * striking the cloud and the emergent - incident distinction does
233  * not apply - code had left emergent intensity at zero, bug reported
234  * by K Korista on discussion group 2007 jun 14 */
235  if( LineSave.ipass > 0 )
236  {
238  }
239 
241  linadd(continuum.cn1216,1215,"Inci",'i',
242  "incident continuum nu*f_nu near Ly-alpha, at illuminated face of cloud");
243  /* see comment concerning parallel code immediately above */
244  if( LineSave.ipass > 0 )
245  {
247  }
248 
251 
252  if( LineSave.ipass > 0 )
253  {
254  continuum.cn4861 = 0.;
255  continuum.cn1216 = 0.;
256  }
257 
261  dense.eden*dense.xIonDense[ipHYDROGEN][1]* 5.45e-12;
262  linadd(flow,0,"Ba C",'i',
263  "integrated Balmer continuum emission");
264 
266  {
273  dense.eden*dense.xIonDense[ipHYDROGEN][1]*3.53e-12;
274  }
275  else
276  {
279  dense.eden*dense.xIonDense[ipHYDROGEN][1]*3.53e-12;
280  }
281  linadd(flow,0,"PA C",'i',
282  "Paschen continuum emission ");
283 
284  /* >>chng 05 oct 16, add this output option, these are a
285  * series of continuum bands defined in the file bands_continuum.dat
286  * this makes it possible to enter any integrated total emission
287  * into the emission-line stack */
288  for( nBand=0; nBand<continuum.nContBand; ++nBand )
289  {
290  double total = 0.;
291  /* find total emission over band - units will be erg cm-2 s-1 */
292  for( i=continuum.ipContBandLow[nBand]; i<continuum.ipContBandHi[nBand]; ++i )
293  {
294  double xIntenOut =
295  /* the attenuated incident continuum */
296  rfield.flux[i-1] +
297 
298  /* the outward emitted continuum radiation field
299  * >>chng 06 aug 07, add geomery.covgeo here */
300  (rfield.ConEmitOut[i-1] +
301 
302  /* outward emitted lines */
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 =
308  ((double)rfield.ConEmitReflec[i-1]/SDIV((double)opac.E2TauAbsFace[i-1]) +
309  /* outward emitted lines */
310  rfield.reflin[i-1] )*geometry.covgeo;
311 
312  /* the fraction of this that gets out */
313  total += rfield.anu[i-1] *
314  emergent_line( xIntenIn , xIntenOut , i )
315  /* >>chng 06 aug 07, add opac.tmn here */
316  / SDIV(opac.tmn[i-1]);
317  }
318  /* we will call lindst with an energy half way between the two
319  * ends of the band. This will make an extinction correction that
320  * has already been applied above by multiplying by emergent_line.
321  * Find this factor and remove it before the call */
322  double corr = emergent_line( 0.5 , 0.5 ,
323  (continuum.ipContBandLow[nBand]+continuum.ipContBandHi[nBand])/2 );
324  if( corr < SMALLFLOAT )
325  total = 0.;
326  else
327  total /= corr;
328 
329  /* convert to physical units */
330  total *= EN1RYD*radius.r1r0sq/radius.dVeff;
331  /* memory not allocated until ipass >= 0 */
332  if( LineSave.ipass > 0 )
333  {
334  LineSv[LineSave.nsum].sumlin[0] = 0.;
335  LineSv[LineSave.nsum].sumlin[1] = 0.;
336  }
337 
338  lindst( total ,
339  continuum.ContBandWavelength[nBand] ,
341  (continuum.ipContBandLow[nBand]+continuum.ipContBandHi[nBand])/2 ,
342  'i' , false,
343  "continuum bands defined in bands_continuum.dat");
344  }
345 
346  linadd(MAX2(0.,CoolHeavy.brems_cool_net),0,"HFFc",'c',
347  "net free-free cooling, ALL species, free-free heating subtracted, so nearly cancels with cooling in LTE ");
348 
349  linadd(MAX2(0.,-CoolHeavy.brems_cool_net),0,"HFFh",'h',
350  "net free-free heating, nearly cancels with cooling in LTE ");
351 
352  linadd(CoolHeavy.brems_cool_h,0,"H FF",'i',
353  " H brems (free-free) cooling ");
354 
355  linadd(CoolHeavy.brems_heat_total,0,"FF H",'i',
356  "total free-free heating ");
357 
358  linadd(CoolHeavy.brems_cool_he,0,"HeFF",'i',
359  "He brems emission ");
360 
361  linadd(CoolHeavy.herec,0,"HeFB",'c',
362  "He recombination cooling ");
363 
364  linadd(CoolHeavy.heavfb,0,"MeFB",'c',
365  "heavy element recombination cooling ");
366 
367  linadd(CoolHeavy.brems_cool_metals,0,"MeFF",'i',
368  "heavy elements (metals) brems cooling, heat not subtracted ");
369 
371  "total brems emission - total cooling but not minus heating ");
372 
374  "part of H brems, in x-ray beyond 0.5KeV ");
375 
376  linadd(CoolHeavy.eebrm,0,"eeff",'c',
377  "electron - electron brems ");
378 
379  /* predict emitted continuum at series of continuum points */
380  /* variables are located in predcont.h,
381  * NPREDCONT number of continuum points,
382  * next two arrays are defined in zerologic.c
383  * EnrPredCont - energy (Ryd) where we want to predict the continuum,
384  * these are set in zerologic.c,
385  * ipPreCont - pointers to energy within continuum array,
386  * these are set in ContCreatePointers
387  *
388  * the array of continuum energies where this will be done is
389  * stored as EnrPredCont which lives in predcont.h
390  * and is initialized in zerologic.c
391  *
392  * the entry nFnu will only be printed if the command
393  * print diffuse continuum
394  * is entered -
395  *
396  * this code should be kept parallel with that in dopunch, where
397  * punch continuum is produced, since two must agree */
398  for( i=0; i < NPREDCONT; i++ )
399  {
400  double SourceTransmitted , Cont_nInu;
401  double SourceReflected,
402  DiffuseOutward,
403  DiffuseInward;
404  double renorm;
405 
406  /* put wavelength in Angstroms into dummy structure, so that we can use iWavLen
407  * to get a proper wavelength with units, continuum energies are stored in EnrPredCont */
409  /*lambda = iWavLen(&TauDummy , &chUnits , &chShift );*/
410 
411  /* >>chng 00 dec 02, there were three occurrences of /opac.tmn which had the
412  * effect of raising the summed continuum by the local opacity correction factor.
413  * in the case of the Lyman continuum this raised the reported value by orders
414  * of magnitude. There have been commented out in the following for now. */
415  /* reflected total continuum (diff+incident emitted inward direction) */
416 
417  /* >>chng 00 dec 08, implement the "set nFnu [SOURCE_REFLECTED] ... command, PvH */
418  /* >>chng 00 dec 19, remove / radius.GeoDil */
419  renorm = rfield.anu2[ipPredCont[i]]*EN1RYD/rfield.widflx[ipPredCont[i]]/*/radius.GeoDil*/;
420 
421  /* this is the reflected diffuse continuum */
422  if( prt.lgDiffuseInward )
423  {
424  DiffuseInward = rfield.ConEmitReflec[ipPredCont[i]]*renorm;
425  /* /opac.tmn[ipPredCont[i]]*/
426  }
427  else
428  {
429  DiffuseInward = 0.;
430  }
431 
432  /* the outward diffuse continuum */
433  if( prt.lgDiffuseOutward )
434  {
435  DiffuseOutward = rfield.ConEmitOut[ipPredCont[i]]*renorm*radius.r1r0sq;
436  }
437  else
438  {
439  DiffuseOutward = 0.;
440  }
441 
442  /* reflected part of INCIDENT continuum (only incident, not diffuse, which was above) */
443  if( prt.lgSourceReflected )
444  {
445  SourceReflected = rfield.ConRefIncid[ipPredCont[i]]*renorm;
446  /* /opac.tmn[ipPredCont[i]]*/
447  }
448  else
449  {
450  SourceReflected = 0.;
451  }
452 
453  /* the attenuated incident continuum */
455  {
456  SourceTransmitted = rfield.flux[ipPredCont[i]]*renorm*radius.r1r0sq;
457  }
458  else
459  {
460  SourceTransmitted = 0.;
461  }
462 
463  /* memory has not been allocated until ipass >= 0, so must not access this element,
464  * this element will be used to save the following quantity */
465  if( LineSave.ipass > 0 )
466  {
467  LineSv[LineSave.nsum].sumlin[0] = 0.;
468  LineSv[LineSave.nsum].sumlin[1] = 0.;
469  }
470 
471  linadd((DiffuseInward+SourceReflected+DiffuseOutward+SourceTransmitted)/radius.dVeff,
472  TauDummy.WLAng,"nFnu",'i',
473  "total continuum at selected energy points " );
474 
475  /* emslin saves the per unit vol emissivity of a line, which is normally
476  * what goes into linadd. We zero this unit emissivity which was set
477  * FOR THE PREVIOUS LINE since it is so situation dependent */
478  if( KILL_CONT && LineSave.ipass > 0 )
479  {
480  LineSv[LineSave.nsum-1].emslin[0] = 0.;
481  LineSv[LineSave.nsum-1].emslin[1] = 0.;
482  }
483 
484  /* this is the normal set to zero to trick the NEXT line into going in properly */
485  if( LineSave.ipass > 0 )
486  {
487  LineSv[LineSave.nsum].sumlin[0] = 0.;
488  LineSv[LineSave.nsum].sumlin[1] = 0.;
489  }
490 
491  /* the nsum-1 -- emslin and nsum -- sumlin is not a bug, look above - they do
492  * different things to different saves */
493  Cont_nInu = rfield.flux[ipPredCont[i]]*renorm*radius.r1r0sq +
494  rfield.ConRefIncid[ipPredCont[i]]*renorm;
495 
496 # if 0
497  /* this code can be used to create assert statements for the continuum shape */
498  if( !i )
499  fprintf(ioQQQ,"\n");
500  char chWL[1000];
501  sprt_wl( chWL , TauDummy.WLAng );
502  fprintf( ioQQQ,"assert line luminosity \"nInu\" %s %.3f\n",
503  chWL ,
504  log10(SDIV(Cont_nInu/radius.dVeff)) + radius.Conv2PrtInten );
505 # endif
506 
507  linadd( Cont_nInu/radius.dVeff,TauDummy.WLAng,"nInu",'i',
508  "transmitted and reflected incident continuum at selected energy points " );
509 
510  /* emslin saves the per unit volume emissivity of a line, which is normally
511  * what goes into linadd. We zero this unit emissivity since it is so situation dependent */
512  if( KILL_CONT && LineSave.ipass > 0 )
513  {
514  LineSv[LineSave.nsum-1].emslin[0] = 0.;
515  LineSv[LineSave.nsum-1].emslin[1] = 0.;
516  }
517 
518  /* memory has not been allocated until ipass >= 0 */
519  if( LineSave.ipass > 0 )
520  {
521  LineSv[LineSave.nsum].sumlin[0] = 0.;
522  LineSv[LineSave.nsum].sumlin[1] = 0.;
523  }
524 
525  linadd( (DiffuseInward+SourceReflected)/radius.dVeff,TauDummy.WLAng,"InwT",'i',
526  "total reflected continuum, total inward emission plus reflected (XXdiffuseXX) total continuum ");
527 
528  if( KILL_CONT && LineSave.ipass > 0 )
529  {
530  LineSv[LineSave.nsum-1].emslin[0] = 0.;
531  LineSv[LineSave.nsum-1].emslin[1] = 0.;
532  }
533 
534  /* memory has not been allocated until ipass >= 0 */
535  if( LineSave.ipass > 0 )
536  {
537  LineSv[LineSave.nsum].sumlin[0] = 0.;
538  LineSv[LineSave.nsum].sumlin[1] = 0.;
539  }
540 
541  linadd(SourceReflected/radius.dVeff,TauDummy.WLAng,"InwC",'i',
542  "reflected incident continuum (only incident) ");
543 
544  if( KILL_CONT && LineSave.ipass > 0 )
545  {
546  LineSv[LineSave.nsum-1].emslin[0] = 0.;
547  LineSv[LineSave.nsum-1].emslin[1] = 0.;
548  }
549  }
550  return;
551 }

Generated for cloudy by doxygen 1.8.1.2