cloudy trunk
Loading...
Searching...
No Matches
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
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}
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
sys_float sexp(sys_float x)
Definition service.cpp:914
const int ipRecEsc
Definition cddefines.h:279
const int LIMELM
Definition cddefines.h:258
const int NISO
Definition cddefines.h:261
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
const int ipRecRad
Definition cddefines.h:283
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
LinSv * LineSv
Definition cdinit.cpp:70
static t_PredCont & Inst()
Definition cddefines.h:175
size_t size() const
Definition predcont.h:28
void set_offset(long offset)
Definition predcont.h:32
t_continuum continuum
Definition continuum.cpp:5
t_CoolHeavy CoolHeavy
Definition coolheavy.cpp:5
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
t_elementnames elementnames
t_geometry geometry
Definition geometry.cpp:5
t_Heavy Heavy
Definition heavy.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH3p
Definition iso.h:31
const int ipH3s
Definition iso.h:30
const int ipH2p
Definition iso.h:29
const int ipH3d
Definition iso.h:32
const int ipH2s
Definition iso.h:28
const int ipH_LIKE
Definition iso.h:62
t_LineSave LineSave
Definition lines.cpp:5
long int StuffComment(const char *chComment)
void linadd(double xInten, realnum wavelength, const char *chLab, char chInfo, const char *chComment)
void lindst(double xInten, realnum wavelength, const char *chLab, long int ipnt, char chInfo, bool lgOutToo, const char *chComment)
double RefIndex(double EnergyWN)
double emergent_line(double emissivity_in, double emissivity_out, long int ipCont)
t_opac opac
Definition opacity.cpp:5
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double EN1RYD
Definition physconst.h:179
UNUSED const double RYDLAM
Definition physconst.h:176
void sprt_wl(char *chString, realnum wl)
Definition prt.cpp:25
t_prt prt
Definition prt.cpp:10
void lines_continuum(void)
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
void chIonLbl(char *chIonLbl_v, const TransitionProxy &t)