cloudy trunk
Loading...
Searching...
No Matches
prt_lines.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 main routine to put emission line intensities into line stack,
4 * calls lineset1, 2, 3, 4 */
5/*Drive_cdLine do the drive cdLine command */
6/*FindStrongestLineLabels find strongest lines contributing to point in continuum energy mesh, output in some save commands */
7#include "cddefines.h"
8#include "physconst.h"
9#include "taulines.h"
10#include "thermal.h"
11#include "yield.h"
12#include "ipoint.h"
13#include "ionbal.h"
14#include "cddrive.h"
15#include "iterations.h"
16#include "trace.h"
17#include "dense.h"
18#include "prt.h"
19#include "rt.h"
20#include "coolheavy.h"
21#include "rfield.h"
22#include "phycon.h"
23#include "elementnames.h"
24#include "iso.h"
25#include "hyperfine.h"
26#include "hydrogenic.h"
27#include "lines_service.h"
28#include "atmdat.h"
29#include "lines.h"
30#include "radius.h"
31
32STATIC void Drive_cdLine( void );
33
34void lines(void)
35{
36 char chLabel[5];
37 long int i,
38 ipnt,
39 nelem;
40 double BigstExtra,
41 ExtraCool,
42 f2, sum;
43
44 DEBUG_ENTRY( "lines()" );
45
46 /* LineSave.ipass
47 * -1 - space not yet allocated - just count number of lines entered into stack
48 * 0 - first call with space allocated - must create labels and add in wavelengths
49 * +1 - later calls - just add intensity
50 */
51
52 /* major routines used here:
53 *
54 * PutLine( tarray )
55 * this uses information in tarray to give various
56 * contributions to lines, and their intensities
57 *
58 * PutExtra( extra )
59 * extra is some extra intensity to add to the line
60 * it will go into the totl contribution put out by PutLine,
61 * and this contribution should be indicated by independent
62 * call to linadd
63 * */
64
65 if( trace.lgTrace )
66 {
67 fprintf( ioQQQ, " lines called\n" );
68 }
69
70 /* the drive cdline command, which checks that all lines can be pulled out by cdLine */
71 if( trace.lgDrv_cdLine && LineSave.ipass > 0 )
73
74 /* total luminosity radiated by this model - will be compared with energy in incident
75 * continuum when calculation is complete */
76 thermal.power += thermal.htot*radius.dVeffAper;
77
78 /* remember the total free-free heating */
79 fixit(); // get rid of brems_heat_total entirely (only net heating is added into stack, so use net here to avoid nonsense ratios elsewhere)
80 //thermal.FreeFreeTotHeat += CoolHeavy.brems_heat_total*radius.dVeffAper;
81 thermal.FreeFreeTotHeat += thermal.heating[0][11]*radius.dVeffAper;
82
83 /* total Compton heating - cooling */
84 rfield.comtot += rfield.cmheat*radius.dVeffAper;
85 thermal.totcol += thermal.ctot*radius.dVeffAper;
86
87 /* up up induced recombination cooling */
88 for( nelem=0; nelem<LIMELM; ++nelem )
89 {
90 hydro.cintot += iso_sp[ipH_LIKE][nelem].RecomInducCool_Rate*radius.dVeffAper;
91 }
92
93 /* nsum is line pointer within large stack of line intensities */
94 LineSave.nsum = 0;
95 LineSave.nComment = 0;
96
97 /* this is used by lindst to proportion inward and outward. should be 50% for
98 * optically thin line. putline sets this to actual value for particular line
99 * and calls lindst then rests to 50% */
100 rt.fracin = 0.5;
101
102 /* last arg in call to lindst and linadd is info on line
103 * info is char variable indicating type of line this is
104 * 'c' cooling
105 * 'h' heating
106 * 'i' information only
107 * 'r' recombination line
108 *
109 * all components of lines are entered into the main line stack here
110 * when printing, filters exist to not print Inwd component */
111
112 /* initialize routine that can generate pointers for forbidden lines,
113 * these are lines that are not transferred otherwise,
114 * in following routines there will be pairs of calls, first to
115 * PntForLine to get pointer, then lindst to add to stack */
116 PntForLine(0.,"FILL",&i);
117
118 /* evaluate rec coefficient for rec lines of C, N, O
119 * some will be used in LineSet2 and then zeroed out,
120 * others left alone and used below */
121 t_ADfA::Inst().rec_lines(phycon.te,LineSave.RecCoefCNO);
122
123 /* initialize ExtraInten with a zero */
124 PutExtra(0.);
125
126 /* put in something impossible in element 0 of line stack */
127 linadd(0.f,0,"zero",'i' , "null placeholder");
128
129 /* this is compared with true volume in final. The number can't
130 * actually be unity since this would overflow on a 32 bit machine */
131 /* integrate the volume as a sanity check */
132 linadd( 1.e-10 , 1 , "Unit" , 'i' , "unit integration placeholder");
133 static long int ipOneAng=-1;
134 if( LineSave.ipass<0 )
135 ipOneAng = ipoint( RYDLAM );
136 lindst( 1.e-10 , 1. , "UntD" , ipOneAng , 'i' , false,"unit integration placeholder");
137
138 /* initial set of general properties */
140
141 /* do all continua */
143
144 /* information about grains */
145 lines_grains();
146
147 /* update all satellite lines */
148 for( long nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
150
151 /* do all hydrogenic ions */
152 lines_hydro();
153
154 /* enter He-iso sequence lines */
155 lines_helium();
156
157#if 0
158 /* This is Ryan's code for dumping lots of Helium lines according to
159 * quantum number rather than wavelength, principally for comparison with Rob
160 * Bauman's code. */
161 if( iteration > 1 )
162 {
163 fprintf( ioQQQ,"ipHi\tipLo\tnu\tlu\tsu\tnl\tll\tsl\tWL\tintens\n" );
164 for( long ipHi=5; ipHi<= iso_sp[ipHE_LIKE][ipHELIUM].numLevels_local - iso_sp[ipHE_LIKE][ipHELIUM].nCollapsed_local; ipHi++ )
165 {
166 for( long ipLo=0; ipLo<ipHi; ipLo++ )
167 {
168 if( iso_sp[ipHE_LIKE][ipHELIUM].trans(ipHi,ipLo).ipCont() > 0 )
169 {
170 double relint, absint;
171
172 if( cdLine("He 1",
173 iso_sp[ipHE_LIKE][ipHELIUM].trans(ipHi,ipLo).WLAng(),
174 &relint, &absint ) )
175 {
176 //iso_sp[ipHE_LIKE][ipHELIUM].trans(ipHi,ipLo).Hi()->chLabel
177
178 if( iso_sp[ipHE_LIKE][ipHELIUM].trans(ipHi,ipLo).WLAng() < 1.1E4 &&
179 iso_sp[ipHE_LIKE][ipHELIUM].trans(ipHi,ipLo).WLAng() > 3.59E3 &&
180 ipLo!=3 && ipLo!=4 && relint >= 0.0009 )
181 {
182 fprintf( ioQQQ,"%li\t%li\t%li\t%li\t%li\t%li\t%li\t%li\t%e\t%e\n",
183 ipHi,
184 ipLo,
185 iso_sp[ipHE_LIKE][ipHELIUM].st[ipHi].n(),
186 iso_sp[ipHE_LIKE][ipHELIUM].st[ipHi].l(),
187 iso_sp[ipHE_LIKE][ipHELIUM].st[ipHi].S(),
188 iso_sp[ipHE_LIKE][ipHELIUM].st[ipLo].n(),
189 iso_sp[ipHE_LIKE][ipHELIUM].st[ipLo].l(),
190 iso_sp[ipHE_LIKE][ipHELIUM].st[ipLo].S(),
191 iso_sp[ipHE_LIKE][ipHELIUM].trans(ipHi,ipLo).WLAng(),
192 relint );
193 }
194 }
195 }
196 }
197 }
198 }
199#endif
200
201 // these must come before the old level 1 or 2 lines. we now have the option
202 // to use the external database by default, or turn it off (set chianti off)
203 // and fall back to the old line treatment. When database is off those lines
204 // do not exist. But when database is on the level 1 lines are still evaluated
205 // but with the ion abundance set to zero. If the level 1 lines were entered
206 // into the stack then searches for the line with cdLine would turn up the level 1
207 // line, which would be zero.
208 // The long term goal is to have all lines be external databases & rm level 1&2 liens
209
210 /* external database lines */
211 i = StuffComment( "database lines" );
212 linadd( 0., (realnum)i , "####", 'i' ,
213 "database lines ");
214 for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
215 {
216 for( EmissionList::iterator em=dBaseTrans[ipSpecies].Emis().begin();
217 em != dBaseTrans[ipSpecies].Emis().end(); ++em)
218 {
219 /* \todo 2 say which database in the comment */
220 if( (*em).Tran().ipCont() > 0)
221 {
222 PutLine((*em).Tran(), "lines from third-party databases", (*(*em).Tran().Hi()).chLabel());
223 }
224 }
225 }
226
227 /* do heavies, lithium through neon */
229
230 /* do heavies, sodium through argon */
232
233 /* do heavies, potassium through zinc */
235
236 /* add up line intensities for certain set of lines */
237 sum = PrtLineSum();
238 /* zero out the location that will receive this information,
239 * remember that memory not allocated until ipass >= 0 */
240 if( LineSave.ipass > 0 )
241 {
242 LineSv[LineSave.nsum].SumLine[0] = 0.;
243 LineSv[LineSave.nsum].SumLine[1] = 0.;
244 }
245 /* optional sum of certain emission lines, set with "print sum" */
246 linadd(sum/radius.dVeffAper,0,"Stoy",'i' ,
247 "Stoy method energy sum ");
248
249 /* next come some recombination lines */
250 i = StuffComment( "recombination" );
251 linadd( 0., (realnum)i , "####", 'i' ,
252 "recombination lines");
253
254 /***********************************************************************
255 * large number of C, N, and O recombination lines *
256 *************************************************************************/
257
258 for( i=0; i < 471; i++ )
259 {
260 /* generate label for the line if ipass is -1 or 0 - saved in arrays
261 * so no need to do this during production */
262 if( LineSave.ipass <= 0 )
263 {
264 /* generate label for the line */
265 strcpy( chLabel, elementnames.chElementSym[(long)(LineSave.RecCoefCNO[0][i])-1] );
266 strcat( chLabel, elementnames.chIonStage[(long)(LineSave.RecCoefCNO[0][i]-
267 LineSave.RecCoefCNO[1][i]+1.01)-1] );
268 }
269 else
270 chLabel[0] = ' ';
271
272 /* number of rec per unit vol
273 * do not predict these simple reccombination intensities at high densities
274 * since lines neglect collisional deexciation and line optical depth effects.
275 * They were not intended for high densities or column densities.
276 * As a result they become unphysically bright at high densities and
277 * violate the black body limit. There would be major
278 * energy conservation problems if they were added in the outward beam in
279 * dense simulations.
280 * */
281 if( dense.eden < 1e8 )
282 {
283 nelem = (long)LineSave.RecCoefCNO[0][i]-1;
284 long int ion = (long)(LineSave.RecCoefCNO[0][i]-LineSave.RecCoefCNO[1][i]+2)-1;
285 f2 = LineSave.RecCoefCNO[3][i]*dense.eden*
286 dense.xIonDense[nelem][ion];
287
288 /* convert to intensity */
289 f2 = f2*1.99e-8/LineSave.RecCoefCNO[2][i];
290 }
291 else
292 {
293 f2 = 0.;
294 }
295 /* stuff it into the stack */
296 PntForLine(LineSave.RecCoefCNO[2][i],chLabel,&ipnt);
297 lindst(f2,LineSave.RecCoefCNO[2][i],chLabel,ipnt, 'r',true ,
298 "recombination line");
299 }
300
301 /* next come the atom_level2 lines */
302 i = StuffComment( "level2 lines" );
303 linadd( 0., (realnum)i , "####", 'i' ,
304 "level2 lines");
305
306 /* add in all the other level 2 wind lines
307 * Dima's 6k lines */
308 ExtraCool = 0.;
309 BigstExtra = 0.;
310 for( i=0; i < nWindLine; i++ )
311 {
312 if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
313 {
314 PutLine(TauLine2[i],"level 2 line");
315 if( TauLine2[i].Coll().cool() > BigstExtra )
316 {
317 BigstExtra = TauLine2[i].Coll().cool();
318 thermal.ipMaxExtra = i+1;
319 }
320 ExtraCool += TauLine2[i].Coll().cool();
321 }
322 }
323 /* keep track of how important this is */
324 thermal.GBarMax = MAX2(thermal.GBarMax,(realnum)(ExtraCool/thermal.ctot));
325
326 /* next come the hyperfine structure lines */
327 i = StuffComment( "hyperfine structure" );
328 linadd( 0., (realnum)i , "####", 'i' ,
329 "hyperfine structure lines ");
330
331 /* this is total cooling due to all HF lines */
332 linadd( hyperfine.cooling_total, 0., "hfin", 'i' ,
333 "total cooling all hyperfine structure lines ");
334
335 /* remember largest local cooling for possible printout in comments */
336 hyperfine.cooling_max = (realnum)MAX2(hyperfine.cooling_max,hyperfine.cooling_total/thermal.ctot);
337
338 /* the hyperfine lines */
339 for( i=0; i < nHFLines; i++ )
340 {
341 PutLine(HFLines[i],
342 "hyperfine structure line");
343 }
344
345 /* next come the inner shell fluorescent lines */
346 i = StuffComment( "inner shell" );
347 linadd( 0., (realnum)i , "####", 'i' ,
348 "inner shell lines");
349
350 /* the group of inner shell fluorescent lines */
351 for( i=0; i < t_yield::Inst().nlines(); ++i )
352 {
353 double xInten =
354 /* density of parent ion, cm-3 */
355 dense.xIonDense[t_yield::Inst().nelem(i)][t_yield::Inst().ion(i)] *
356 /* photo rate per atom per second, s-1 */
357 ionbal.PhotoRate_Shell[t_yield::Inst().nelem(i)][t_yield::Inst().ion(i)][t_yield::Inst().nshell(i)][0]*
358 /* fluor yield - dimensionless */
359 t_yield::Inst().yield(i) *
360 /* photon energy - ryd, converted into ergs */
362
363 /* create label if initializing line stack */
364 if( LineSave.ipass == 0 )
365 {
366 /* only generate the line label if it is going to be used */
367 strcpy( chLabel , elementnames.chElementSym[t_yield::Inst().nelem(i)] );
368 strcat( chLabel , elementnames.chIonStage[t_yield::Inst().ion_emit(i)] );
369# if 0
370 /* only print yields for atoms */
371 if( t_yield::Inst().ion(i) == 0 && t_yield::Inst().nelem()(i) == ipIRON )
372 fprintf(ioQQQ,"DEBUGyeild\t%s\t%.3f\t%.3e\n",
373 /* line designation, energy in eV, yield */
374 chLabel , t_yield::Inst().energy()(i)*EVRYD, t_yield::Inst().yield(i) );
375# endif
376 }
377
378 /* the group of inner shell fluorescent lines */
379 lindst(
380 /* intensity of line */
381 xInten,
382 /* wavelength of line in Angstroms */
383 (realnum)RYDLAM / t_yield::Inst().energy(i),
384 /* label */
385 chLabel ,
386 /* continuum array offset for line as set in ipoint */
387 t_yield::Inst().ipoint(i),
388 /* type of line - count as a recombination line */
389 'r',
390 /* include line in continuum? */
391 true ,
392 "inner shell line");
393 }
394
395 /* >>chng 06 jan 03, confirm that number of lines never changes once we have
396 * created the labels */
397 {
398 static long nLineSave=-1 , ndLineSave=-1;
399 if( LineSave.ipass == 0 )
400 {
401 nLineSave = LineSave.nsum;
402 ndLineSave = LineSave.nsum;
403 }
404 else if( LineSave.ipass > 0 )
405 {
406 /* this can't happen */
407 if( nLineSave<= 0 || ndLineSave < 0 )
409
410 /* now make sure that we have the same number of lines as we had previously
411 * created labels. This would not pass if we did not add in exactly the same
412 * number of lines on each pass */
413 if( nLineSave != LineSave.nsum )
414 {
415 fprintf( ioQQQ, "DISASTER number of lines in LineSave.nsum changed between pass 0 and 1 - this is impossible\n" );
416 fprintf( ioQQQ, "DISASTER LineSave.nsum is %li and nLineSave is %li\n",
417 LineSave.nsum ,
418 nLineSave);
419 ShowMe();
421 }
422 if( ndLineSave != LineSave.nsum )
423 {
424 fprintf( ioQQQ, "DISASTER number of lines in LineSave.nsum changed between pass 0 and 1 - this is impossible\n" );
425 fprintf( ioQQQ, "DISASTER LineSave.nsum is %li and ndLineSave is %li\n",
426 LineSave.nsum ,
427 ndLineSave);
428 ShowMe();
430 }
431 }
432 }
433
434 /* now do all molecules - do last since so many H2 lines */
436
437 if( trace.lgTrace )
438 {
439 fprintf( ioQQQ, " lines returns\n" );
440 }
441 return;
442}
443
444/*Drive_cdLine do the drive cdLine command */
446{
447 long int j;
448 bool lgMustPrintHeader = true;
449 double absval , rel;
450
451 DEBUG_ENTRY( "Drive_cdLine()" );
452
453 for( j=1; j < LineSave.nsum; j++ )
454 {
455 if( cdLine( LineSv[j].chALab , LineSv[j].wavelength , &absval , &rel ) <= 0 )
456 {
457 /* print header if first miss */
459 {
460 fprintf(ioQQQ,"n\tlab\twl\n");
461 lgMustPrintHeader = false;
462 }
463
464 fprintf(ioQQQ,"%li\t%s\t%f\n", j, LineSv[j].chALab , LineSv[j].wavelength );
465 }
466 }
467 fprintf( ioQQQ, " Thanks for checking on the cdLine routine!\n" );
469}
470
FILE * ioQQQ
Definition cddefines.cpp:7
long int iteration
Definition cddefines.cpp:16
#define STATIC
Definition cddefines.h:97
const int ipIRON
Definition cddefines.h:330
const int LIMELM
Definition cddefines.h:258
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int NISO
Definition cddefines.h:261
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
NORETURN void TotalInsanity(void)
Definition service.cpp:886
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void ShowMe(void)
Definition service.cpp:181
void fixit(void)
Definition service.cpp:991
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition cddrive.cpp:1228
long nWindLine
Definition cdinit.cpp:19
LinSv * LineSv
Definition cdinit.cpp:70
EmissionProxy::iterator iterator
Definition emission.h:317
static t_ADfA & Inst()
Definition cddefines.h:175
void rec_lines(double t, realnum r[][471])
realnum yield(long n) const
Definition yield.h:65
realnum energy(long n) const
Definition yield.h:63
int nlines() const
Definition yield.h:76
int nshell(long n) const
Definition yield.h:61
int nelem(long n) const
Definition yield.h:59
int ion(long n) const
Definition yield.h:60
long ipoint(double energy_ryd)
t_dense dense
Definition dense.cpp:24
t_elementnames elementnames
t_hydro hydro
Definition hydrogenic.cpp:5
t_hyperfine hyperfine
Definition hyperfine.cpp:5
t_ionbal ionbal
Definition ionbal.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipHE_LIKE
Definition iso.h:63
void iso_satellite_update(long nelem)
const int ipH_LIKE
Definition iso.h:62
t_LineSave LineSave
Definition lines.cpp:5
long int StuffComment(const char *chComment)
void lines_molecules(void)
void lines_helium(void)
void lines_lv1_li_ne(void)
void lines_lv1_k_zn(void)
void lines_general(void)
void lines_continuum(void)
void lines_lv1_na_ar(void)
void lines_hydro(void)
void lines_grains(void)
void linadd(double xInten, realnum wavelength, const char *chLab, char chInfo, const char *chComment)
void PntForLine(double wavelength, const char *chLabel, long int *ipnt)
void lindst(double xInten, realnum wavelength, const char *chLab, long int ipnt, char chInfo, bool lgOutToo, const char *chComment)
static realnum * wavelength
#define S(I_, J_)
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double EN1RYD
Definition physconst.h:179
UNUSED const double EVRYD
Definition physconst.h:189
UNUSED const double RYDLAM
Definition physconst.h:176
double PrtLineSum(void)
void lines(void)
Definition prt_lines.cpp:34
STATIC void Drive_cdLine(void)
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
t_rt rt
Definition rt.cpp:5
static bool lgMustPrintHeader
long int nSpecies
Definition taulines.cpp:21
TransitionList TauLine2("TauLine2", &AnonStates)
vector< TransitionList > dBaseTrans
Definition taulines.cpp:17
TransitionList HFLines("HFLines", &AnonStates)
long int nHFLines
Definition taulines.cpp:31
t_thermal thermal
Definition thermal.cpp:5
t_trace trace
Definition trace.cpp:5
void PutExtra(double Extra)
void PutLine(const TransitionProxy &t, const char *chComment, const char *chLabelTemp)