cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
prt_lines.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 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 #include "cddefines.h"
7 #include "physconst.h"
8 #include "taulines.h"
9 #include "thermal.h"
10 #include "yield.h"
11 #include "ionbal.h"
12 #include "cddrive.h"
13 #include "trace.h"
14 #include "dense.h"
15 #include "prt.h"
16 #include "rt.h"
17 #include "coolheavy.h"
18 #include "rfield.h"
19 #include "phycon.h"
20 #include "elementnames.h"
21 #include "iso.h"
22 #include "hyperfine.h"
23 #include "hydrogenic.h"
24 #include "lines_service.h"
25 #include "atmdat.h"
26 #include "lines.h"
27 #include "radius.h"
28 STATIC void Drive_cdLine( void );
29 
30 void lines(void)
31 {
32  char chLabel[5];
33  static bool lgRecOn;
34  long int i,
35  ipnt,
36  nelem;
37  double BigstExtra,
38  ExtraCool,
39  f2, sum;
40 
41  DEBUG_ENTRY( "lines()" );
42 
43  /* LineSave.ipass
44  * -1 - space not yet allocated - just count number of lines entered into stack
45  * 0 - first call with space allocated - must create labels and add in wavelengths
46  * +1 - later calls - just add intensity
47  */
48 
49  /* major routines used here:
50  *
51  * PutLine( tarray )
52  * this uses information in tarray to give various
53  * contributions to lines, and their intensities
54  *
55  * PutExtra( extra )
56  * extra is some extra intensity to add to the line
57  * it will go into the totl contribution put out by PutLine,
58  * and this contribution should be indicated by independent
59  * call to linadd
60  * */
61 
62  if( trace.lgTrace )
63  {
64  fprintf( ioQQQ, " lines called\n" );
65  }
66 
67  /* the drive cdline command, which checks that all lines can be pulled out by cdLine */
68  if( trace.lgDrv_cdLine && LineSave.ipass > 0 )
69  Drive_cdLine();
70 
71  /* total luminosity radiated by this model - will be compared with energy in incident
72  * continuum when calculation is complete */
74 
75  /* remember the total free-free heating */
77 
78  /* total Compton heating - cooling */
81 
82  /* up up induced recombination cooling */
83  for( nelem=0; nelem<LIMELM; ++nelem )
84  {
86  }
87 
88  /* nsum is line pointer within large stack of line intensities */
89  LineSave.nsum = 0;
90  LineSave.nComment = 0;
91 
92  /* this is used by lindst to proportion inward and outward. should be 50% for
93  * optically thin line. putline sets this to actual value for particular line
94  * and calls lindst then rests to 50% */
95  rt.fracin = 0.5;
96 
97  /* last arg in call to lindst and linadd is info on line
98  * info is char variable indicating type of line this is
99  * 'c' cooling
100  * 'h' heating
101  * 'i' information only
102  * 'r' recombination line
103  *
104  * all components of lines are entered into the main line stack here
105  * when printing, filters exist to not print Inwd component */
106 
107  /* initialize routine that can generate pointers for forbidden lines,
108  * these are lines that are not transferred otherwise,
109  * in following routines there will be pairs of calls, first to
110  * PntForLine to get pointer, then lindst to add to stack */
111  PntForLine(0.,"FILL",&i);
112 
113  /* evaluate rec coefficient for rec lines of C, N, O
114  * some will be used in LineSet2 and then zeroed out,
115  * others left alone and used below */
117 
118  /* initialize ExtraInten with a zero */
119  PutExtra(0.);
120 
121  /* put in something impossible in element 0 of line stack */
122  linadd(0.f,0,"zero",'i' , "null placeholder");
123 
124  /* this is compared with true volume in final. The number can't
125  * actually be unity since this would overflow on a 32 bit machine */
126  /* integrate the volume as a sanity check */
127  linadd( 1.e-10 , 1 , "Unit" , 'i' , "unit integration placeholder");
128 
129  /* initial set of general properties */
130  lines_general();
131 
132  /* do all continua */
133  lines_continuum();
134 
135  /* information about grains */
136  lines_grains();
137 
138  /* update all satellite lines */
140 
141  /* do all hydrogenic ions */
142  lines_hydro();
143 
144  /* enter He-iso sequence lines */
145  lines_helium();
146 
147 #if 0
148  /* This is Ryan's code for dumping lots of Helium lines according to
149  * quantum number rather than wavelength, principally for comparison with Rob
150  * Bauman's code. */
151  if( iteration > 1 )
152  {
153  fprintf( ioQQQ,"ipHi\tipLo\tnu\tlu\tsu\tnl\tll\tsl\tWL\tintens\n" );
154  for( long ipHi=5; ipHi<= iso.numLevels_local[ipHE_LIKE][ipHELIUM] - iso.nCollapsed_local[ipHE_LIKE][ipHELIUM]; ipHi++ )
155  {
156  for( long ipLo=0; ipLo<ipHi; ipLo++ )
157  {
158  if( Transitions[ipHE_LIKE][ipHELIUM][ipHi][ipLo].ipCont > 0 )
159  {
160  double relint, absint;
161 
162  if( cdLine("He 1",
163  Transitions[ipHE_LIKE][ipHELIUM][ipHi][ipLo].WLAng,
164  &relint, &absint ) )
165  {
166  //Transitions[ipHE_LIKE][ipHELIUM][ipHi][ipLo].Hi->chLabel
167 
168  if( Transitions[ipHE_LIKE][ipHELIUM][ipHi][ipLo].WLAng < 1.1E4 &&
169  Transitions[ipHE_LIKE][ipHELIUM][ipHi][ipLo].WLAng > 3.59E3 &&
170  ipLo!=3 && ipLo!=4 && relint >= 0.0009 )
171  {
172  fprintf( ioQQQ,"%li\t%li\t%li\t%li\t%li\t%li\t%li\t%li\t%e\t%e\n",
173  ipHi,
174  ipLo,
175  StatesElem[ipHE_LIKE][ipHELIUM][ipHi].n,
176  StatesElem[ipHE_LIKE][ipHELIUM][ipHi].l,
177  StatesElem[ipHE_LIKE][ipHELIUM][ipHi].S,
178  StatesElem[ipHE_LIKE][ipHELIUM][ipLo].n,
179  StatesElem[ipHE_LIKE][ipHELIUM][ipLo].l,
180  StatesElem[ipHE_LIKE][ipHELIUM][ipLo].S,
181  Transitions[ipHE_LIKE][ipHELIUM][ipHi][ipLo].WLAng,
182  relint );
183  }
184  }
185  }
186  }
187  }
188  }
189 #endif
190 
191  /* do heavies, lithium through neon */
192  lines_lv1_li_ne();
193 
194  /* do heavies, sodium through argon */
195  lines_lv1_na_ar();
196 
197  /* do heavies, potassium through zinc */
198  lines_lv1_k_zn();
199 
200  /* add up line intensities for certain set of lines */
201  sum = PrtLineSum( " SUM" );
202  /* zero out the location that will receive this information,
203  * remember that memory not allocated until ipass >= 0 */
204  if( LineSave.ipass > 0 )
206  /* optional sum of certain emission lines, set with "print sum" */
207  linadd(sum/radius.dVeff,0,"Stoy",'i' ,
208  "Stoy method energy sum ");
209 
210  /* next come some recombination lines */
211  i = StuffComment( "recombination" );
212  linadd( 0., (realnum)i , "####", 'i' ,
213  "recombination lines");
214 
215  /***********************************************************************
216  * large number of C, N, and O recombination lines *
217  *************************************************************************/
218 
219  if( LineSave.ipass <= 0 )
220  {
221  /* initialize to true, may be set false if density is very high below */
222  lgRecOn = true;
223  }
224 
225  for( i=0; i < 471; i++ )
226  {
227  /* generate label for the line */
228  strcpy( chLabel, elementnames.chElementSym[(long)(LineSave.RecCoefCNO[0][i])-1] );
229  strcat( chLabel, elementnames.chIonStage[(long)(LineSave.RecCoefCNO[0][i]-
230  LineSave.RecCoefCNO[1][i]+1.01)-1] );
231 
232  /* number of rec per unit vol
233  * >>chng 97 aug 28, do not predict rec intensities at high densities
234  * blr.in and oldblr.in go nuts if this is added in outward beam
235  * these recombination coefficients were computed in the low density limit -
236  * they were not intended for high densities */
237  if( dense.eden < 1e8 && lgRecOn )
238  {
239  f2 = LineSave.RecCoefCNO[3][i]*dense.eden*
240  dense.xIonDense[(long)(LineSave.RecCoefCNO[0][i])-1][(long)(LineSave.RecCoefCNO[0][i]-LineSave.RecCoefCNO[1][i]+2)-1];
241 
242  /* convert to intensity */
243  f2 = f2*1.99e-8/LineSave.RecCoefCNO[2][i];
244  }
245  else
246  {
247  lgRecOn = false;
248  f2 = 0.;
249  }
250  /* finally stuff it into the stack */
251  PntForLine(LineSave.RecCoefCNO[2][i],chLabel,&ipnt);
252  lindst(f2,LineSave.RecCoefCNO[2][i],chLabel,ipnt, 'r',true ,
253  "recombination line");
254  }
255 
256  /* next come the atom_level2 lines */
257  i = StuffComment( "level2 lines" );
258  linadd( 0., (realnum)i , "####", 'i' ,
259  "level2 lines");
260 
261  /* add in all the other level 2 wind lines
262  * Dima's 6k lines */
263  ExtraCool = 0.;
264  BigstExtra = 0.;
265  for( i=0; i < nWindLine; i++ )
266  {
267  if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
268  {
269  PutLine(&TauLine2[i],"level 2 line");
270  if( TauLine2[i].Coll.cool > BigstExtra )
271  {
272  BigstExtra = TauLine2[i].Coll.cool;
273  thermal.ipMaxExtra = i+1;
274  }
275  ExtraCool += TauLine2[i].Coll.cool;
276  }
277  }
278  /* keep track of how important this is */
280 
281  /* next come the hyperfine structure lines */
282  i = StuffComment( "hyperfine structure" );
283  linadd( 0., (realnum)i , "####", 'i' ,
284  "hyperfine structure lines ");
285 
286  /* this is total cooling due to all HF lines */
287  linadd( hyperfine.cooling_total, 0., "hfin", 'i' ,
288  "total cooling all hyperfine structure lines ");
289 
290  /* remember largest local cooling for possible printout in comments */
292 
293  /* the hyperfine lines */
294  for( i=0; i < nHFLines; i++ )
295  {
296  PutLine(&HFLines[i],
297  "hyperfine structure line");
298  }
299 
300  /* Databases: Atoms & Molecules: Added Oct 26 07*/
301  i = StuffComment( "database lines" );
302  linadd( 0., (realnum)i , "####", 'i' ,
303  "database lines ");
304 
305  /* Atoms & Molecules */
306  for( i=0; i < linesAdded2; i++ )
307  {
308  /* \todo 2 say which database in the comment */
309  PutLine(atmolEmis[i].tran, "lines from third-party databases", atmolEmis[i].tran->Hi->chLabel);
310  }
311 
312  /* next come the inner shell fluorescent lines */
313  i = StuffComment( "inner shell" );
314  linadd( 0., (realnum)i , "####", 'i' ,
315  "inner shell lines");
316 
317  /* the group of inner shell fluorescent lines */
318  for( i=0; i < t_yield::Inst().nlines(); ++i )
319  {
320  double xInten =
321  /* density of parent ion, cm-3 */
323  /* photo rate per atom per second, s-1 */
325  /* fluor yield - dimensionless */
326  t_yield::Inst().yield(i) *
327  /* photon energy - ryd, converted into ergs */
328  t_yield::Inst().energy(i) * EN1RYD;
329 
330  /* create label if initializing line stack */
331  if( LineSave.ipass == 0 )
332  {
333  /* only generate the line label if it is going to be used */
334  strcpy( chLabel , elementnames.chElementSym[t_yield::Inst().nelem(i)] );
335  strcat( chLabel , elementnames.chIonStage[t_yield::Inst().ion_emit(i)] );
336 # if 0
337  /* only print yields for atoms */
338  if( t_yield::Inst().ion(i) == 0 && t_yield::Inst().nelem(i) == ipIRON )
339  fprintf(ioQQQ,"DEBUGyeild\t%s\t%.3f\t%.3e\n",
340  /* line designation, energy in eV, yield */
341  chLabel , t_yield::Inst().energy(i)*EVRYD, t_yield::Inst().yield(i) );
342 # endif
343  }
344 
345  /* the group of inner shell fluorescent lines */
346  lindst(
347  /* intensity of line */
348  xInten,
349  /* wavelength of line in Angstroms */
350  (realnum)RYDLAM / t_yield::Inst().energy(i),
351  /* label */
352  chLabel ,
353  /* continuum array offset for line as set in ipoint */
354  t_yield::Inst().ipoint(i),
355  /* type of line - count as a recombination line */
356  'r',
357  /* include line in continuum? */
358  true ,
359  "inner shell line");
360  }
361 
362  /* >>chng 06 jan 03, confirm that number of lines never changes once we have
363  * created the labels */
364  {
365  static long nLineSave=-1 , ndLineSave=-1;
366  if( LineSave.ipass == 0 )
367  {
368  nLineSave = LineSave.nsum;
369  ndLineSave = LineSave.nsum;
370  }
371  else if( LineSave.ipass > 0 )
372  {
373  /* this can't happen */
374  if( nLineSave<= 0 || ndLineSave < 0 )
375  TotalInsanity();
376 
377  /* now make sure that we have the same number of lines as we had previously
378  * created labels. This would not pass if we did not add in exactly the same
379  * number of lines on each pass */
380  if( nLineSave != LineSave.nsum )
381  {
382  fprintf( ioQQQ, "DISASTER number of lines in LineSave.nsum changed between pass 0 and 1 - this is impossible\n" );
383  fprintf( ioQQQ, "DISASTER LineSave.nsum is %li and nLineSave is %li\n",
384  LineSave.nsum ,
385  nLineSave);
386  ShowMe();
387  cdEXIT(EXIT_FAILURE);
388  }
389  if( ndLineSave != LineSave.nsum )
390  {
391  fprintf( ioQQQ, "DISASTER number of lines in LineSave.nsum changed between pass 0 and 1 - this is impossible\n" );
392  fprintf( ioQQQ, "DISASTER LineSave.nsum is %li and ndLineSave is %li\n",
393  LineSave.nsum ,
394  ndLineSave);
395  ShowMe();
396  cdEXIT(EXIT_FAILURE);
397  }
398  }
399  }
400 
401  /* now do all molecules - do last since so many H2 lines */
402  lines_molecules();
403 
404  if( trace.lgTrace )
405  {
406  fprintf( ioQQQ, " lines returns\n" );
407  }
408  return;
409 }
410 
411 /*Drive_cdLine do the drive cdLine command */
412 STATIC void Drive_cdLine( void )
413 {
414  long int j;
415  bool lgMustPrintHeader = true;
416  double absval , rel;
417 
418  DEBUG_ENTRY( "Drive_cdLine()" );
419 
420  for( j=1; j < LineSave.nsum; j++ )
421  {
422  if( cdLine( LineSv[j].chALab , LineSv[j].wavelength , &absval , &rel ) <= 0 )
423  {
424  /* print header if first miss */
425  if( lgMustPrintHeader )
426  {
427  fprintf(ioQQQ,"n\tlab\twl\n");
428  lgMustPrintHeader = false;
429  }
430 
431  fprintf(ioQQQ,"%li\t%s\t%f\n", j, LineSv[j].chALab , LineSv[j].wavelength );
432  }
433  }
434  fprintf( ioQQQ, " Thanks for checking on the cdLine routine!\n" );
435  cdEXIT(EXIT_FAILURE);
436 }

Generated for cloudy by doxygen 1.8.1.2