cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
lines_service.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 /*PutLine enter local line intensity into the intensity stack for eventual printout */
4 /*PutExtra enter and 'extra' intensity source for some line */
5 /*linadd enter lines into the line storage array, called once per zone */
6 /*DumpLine print various information about an emission line vector,
7  * used in debugging, print to std out, ioQQQ */
8 /*TexcLine derive excitation temperature of line from contents of line array */
9 /*GetGF convert Einstein A into oscillator strength */
10 /*emit_frac returns fraction of populations the produce emission */
11 /*abscf convert gf into absorption coefficient */
12 /*chIonLbl use information in line array to generate a null terminated ion label in "Fe 2" */
13 /*chLineLbl use information in line transfer arrays to generate a line label */
14 /*RefIndex calculates the index of refraction of air using the line energy in wavenumbers,
15  * used to convert vacuum wavelengths to air wavelengths. */
16 /*eina convert a gf into an Einstein A */
17 /*WavlenErrorGet - find difference between two wavelengths */
18 /*PutCS enter a collision strength into an individual line vector */
19 /*lindst add local line intensity to line luminosity stack */
20 /*PntForLine generate pointer for forbidden line */
21 /*TransitionZero zeros out TransitionZero */
22 /*EmLineJunk set all elements of transition struc to dangerous values */
23 /*EmLineZero set all elements of transition struc to zero */
24 /*LineConvRate2CS convert down coll rate back into electron cs in case other parts of code need this for reference */
25 /*lgTauGood returns true is we have not overrun optical depth scale */
26 /*OccupationNumberLine - derive the photon occupation number at line center for any line */
27 /*outline - adds line photons to reflin and outlin */
28 /*MakeCS compute collision strength by g-bar approximations */
29 /*gbar1 compute g-bar collision strength using Mewe approximations */
30 /*gbar0 compute g-bar gaunt factor for neutrals */
31 /*totlin sum total intensity of cooling, recombination, or intensity lines */
32 /*FndLineHt search through line heat arrays to find the strongest heat source */
33 /*ConvRate2CS convert down coll rate back into electron cs in case other parts of code need this for reference */
34 #include "cddefines.h"
35 #include "physconst.h"
36 #include "phycon.h"
37 #include "lines.h"
38 #include "radius.h"
39 #include "geometry.h"
40 #include "elementnames.h"
41 #include "rt.h"
42 #include "dense.h"
43 #include "rfield.h"
44 #include "opacity.h"
45 #include "ipoint.h"
46 #include "iso.h"
47 #include "taulines.h"
48 #include "hydrogenic.h"
49 #include "lines_service.h"
50 
51 /*outline - adds line photons to reflin and outlin */
52 void outline( transition *t )
53 {
54  long int ip = t->ipCont-1;
55  double xInWrd = t->Emis->phots*t->Emis->FracInwd;
56 
57  DEBUG_ENTRY( "outline()" );
58 
59  ASSERT( t->Emis->phots >= 0. );
60  ASSERT( t->Emis->FracInwd >= 0. );
61  ASSERT( radius.BeamInIn >= 0. );
62  ASSERT( radius.BeamInOut >= 0. );
63 
64  /* the reflected part */
65  rfield.reflin[ip] += (realnum)(xInWrd*radius.BeamInIn);
66 
67  /* inward beam that goes out since sphere set */
68  rfield.outlin[ip] += (realnum)(xInWrd*radius.BeamInOut*opac.tmn[ip]*
69  t->Emis->ColOvTot);
70 
71  /* outward part */
73  rfield.outlin[ip] += (realnum)(t->Emis->phots*
74  (1. - t->Emis->FracInwd)*radius.BeamOutOut* t->Emis->ColOvTot);
75  return;
76 }
77 
78 /*emit_frac returns fraction of populations the produce emission */
79 double emit_frac( transition *t )
80 {
81  DEBUG_ENTRY( "emit_frac()" );
82 
83  ASSERT( t->ipCont > 0 );
84 
85  /* collisional deexcitation and destruction by background opacities
86  * are loss of photons without net emission */
87  double deexcit_loss = t->Coll.cs * dense.cdsqte + t->Emis->Aul*t->Emis->Pdest;
88  /* this is what is observed */
89  double rad_deexcit = t->Emis->Aul*(t->Emis->Pelec_esc + t->Emis->Pesc);
90  return rad_deexcit/(deexcit_loss + rad_deexcit);
91 }
92 
93 /*DumpLine print various information about an emission line vector,
94  * used in debugging, print to std out, ioQQQ */
96 {
97  char chLbl[11];
98 
99  DEBUG_ENTRY( "DumpLine()" );
100 
101  ASSERT( t->ipCont > 0 );
102 
103  /* routine to print contents of line arrays */
104  strcpy( chLbl, chLineLbl(t) );
105 
106  fprintf( ioQQQ,
107  " %10.10s Te%.2e eden%.1e CS%.2e Aul%.1e Tex%.2e cool%.1e het%.1e conopc%.1e albdo%.2e\n",
108  chLbl,
109  phycon.te,
110  dense.eden,
111  t->Coll.cs,
112  t->Emis->Aul,
113  TexcLine(t),
114  t->Coll.cool,
115  t->Coll.heat ,
116  opac.opacity_abs[t->ipCont-1],
117  opac.albedo[t->ipCont-1]);
118 
119  fprintf( ioQQQ,
120  " Tin%.1e Tout%.1e Esc%.1e eEsc%.1e DesP%.1e Pump%.1e OTS%.1e PopL,U %.1e %.1e PopOpc%.1e\n",
121  t->Emis->TauIn,
122  t->Emis->TauTot,
123  t->Emis->Pesc,
124  t->Emis->Pelec_esc,
125  t->Emis->Pdest,
126  t->Emis->pump,
127  t->Emis->ots,
128  t->Lo->Pop,
129  t->Hi->Pop ,
130  t->Emis->PopOpc );
131  return;
132 }
133 
134 
135 /*OccupationNumberLine - derive the photon occupation number at line center for any line */
137 {
138  double OccupationNumberLine_v;
139 
140  DEBUG_ENTRY( "OccupationNumberLine()" );
141 
142  ASSERT( t->ipCont > 0 );
143 
144  /* routine to evaluate line photon occupation number */
145  if( t->Lo->Pop > SMALLFLOAT )
146  {
147  /* the lower population with correction for stimulated emission */
148  double PopLo_corr = t->Lo->Pop / t->Lo->g - t->Hi->Pop / t->Hi->g;
149  OccupationNumberLine_v = ( t->Hi->Pop / t->Hi->g )/SDIV(PopLo_corr ) *
150  (1. - t->Emis->Pesc);
151  }
152  else
153  {
154  OccupationNumberLine_v = 0.;
155  }
156  return( OccupationNumberLine_v );
157 }
158 
159 /*TexcLine derive excitation temperature of line from contents of line array */
160 double TexcLine(transition * t)
161 {
162  double TexcLine_v;
163 
164  DEBUG_ENTRY( "TexcLine()" );
165 
166  /* routine to evaluate line excitation temp using contents of line array
167  * */
168  if( t->Hi->Pop * t->Lo->Pop > 0. )
169  {
170  TexcLine_v = ( t->Hi->Pop / t->Hi->g )/( t->Lo->Pop / t->Lo->g );
171  TexcLine_v = log(TexcLine_v);
172  /* protect against infinite temp limit */
173  if( fabs(TexcLine_v) > SMALLFLOAT )
174  {
175  TexcLine_v = - t->EnergyK / TexcLine_v;
176  }
177  }
178  else
179  {
180  TexcLine_v = 0.;
181  }
182  return( TexcLine_v );
183 }
184 
185 /*eina convert a gf into an Einstein A */
186 double eina(double gf,
187  double enercm,
188  double gup)
189 {
190  double eina_v;
191 
192  DEBUG_ENTRY( "eina()" );
193 
194  /* derive the transition prob, given the following
195  * call to function is gf, energy in cm^-1, g_up
196  * gf is product of g and oscillator strength
197  * eina = ( gf / 1.499e-8 ) / (wl/1e4)**2 / gup */
198  eina_v = (gf/gup)*TRANS_PROB_CONST*POW2(enercm);
199  return( eina_v );
200 }
201 
202 /*GetGF convert Einstein A into oscillator strength */
203 double GetGF(double trans_prob,
204  double enercm,
205  double gup)
206 {
207  double GetGF_v;
208 
209  DEBUG_ENTRY( "GetGF()" );
210 
211  ASSERT( enercm > 0. );
212  ASSERT( trans_prob > 0. );
213  ASSERT( gup > 0.);
214 
215  /* derive the transition prob, given the following
216  * call to function is gf, energy in cm^-1, g_up
217  * gf is product of g and oscillator strength
218  * trans_prob = ( GetGF/gup) / 1.499e-8 / ( 1e4/enercm )**2 */
219  GetGF_v = trans_prob*gup/TRANS_PROB_CONST/POW2(enercm);
220  return( GetGF_v );
221 }
222 
223 /*abscf convert gf into absorption coefficient */
224 double abscf(double gf,
225  double enercm,
226  double gl)
227 {
228  double abscf_v;
229 
230  DEBUG_ENTRY( "abscf()" );
231 
232  ASSERT(gl > 0. && enercm > 0. && gf > 0. );
233 
234  /* derive line absorption coefficient, given the following:
235  * gf, enercm, g_low
236  * gf is product of g and oscillator strength */
237  abscf_v = 1.4974e-6*(gf/gl)*(1e4/enercm);
238  return( abscf_v );
239 }
240 
241 /*chIonLbl use information in line array to generate a null terminated ion label in "Fe 2" */
242 void chIonLbl(char *chIonLbl_v , transition * t )
243 {
244  /*static char chIonLbl_v[5];*/
245 
246  DEBUG_ENTRY( "chIonLbl()" );
247 
248  /* function to use information within the line array
249  * to generate an ion label, giving element and
250  * ionization stage
251  * */
252  if( t->Hi->nelem < 0 )
253  {
254  /* this line is to be ignored */
255  strcpy( chIonLbl_v, "Dumy" );
256  }
257  else if( t->Hi->nelem-1 >= LIMELM )
258  {
259  /* this is one of the molecules, either 12CO or 13CO */
260 
261  /* >>chng 02 may 15, from to chElementNameShort to go mole right */
262  strcpy( chIonLbl_v , elementnames.chElementNameShort[t->Hi->nelem-1] );
263 
264  /* chIonStage is four char null terminated string, starting with "_1__"
265  strcat( chIonLbl_v, "CO");*/
266  }
267 
268  else
269  {
270  ASSERT( t->Hi->nelem >= 1 );
271  /* ElmntSym.chElementSym is null terminated, 2 ch + null, string giving very
272  * short form of element name */
273  strcpy( chIonLbl_v , elementnames.chElementSym[t->Hi->nelem -1] );
274 
275  /* chIonStage is four char null terminated string, starting with "_1__" */
276  strcat( chIonLbl_v, elementnames.chIonStage[t->Hi->IonStg-1]);
277  }
278  /* chIonLbl is four char null terminated string */
279  return/*( chIonLbl_v )*/;
280 }
281 
282 /*chLineLbl use information in line transfer arrays to generate a line label */
283 /* this label is null terminated */
284 /* ContCreatePointers has test this with full range of wavelengths */
286 {
287  static char chLineLbl_v[11];
288 
289  DEBUG_ENTRY( "chLineLbl()" );
290 
291 
292  /* function to use information within the line array
293  * to generate a line label, giving element and
294  * ionization stage
295  * rhs are set in large block data */
296 
297  /* NB this function is profoundly slow due to sprintf statement
298  * also - it cannot be evaluated within a write statement itself*/
299 
300  if( t->WLAng > (realnum)INT_MAX )
301  {
302  sprintf( chLineLbl_v, "%2.2s%2.2s%5i%c",
305  (int)(t->WLAng/1e8), 'c' );
306  }
307  else if( t->WLAng > 9999999. )
308  {
309  /* wavelength is very large, convert to centimeters */
310  sprintf( chLineLbl_v, "%2.2s%2.2s%5.2f%c",
313  t->WLAng/1e8, 'c' );
314  }
315  else if( t->WLAng > 999999. )
316  {
317  /* wavelength is very large, convert to microns */
318  sprintf( chLineLbl_v, "%2.2s%2.2s%5i%c",
321  (int)(t->WLAng/1e4), 'm' );
322  }
323  else if( t->WLAng > 99999. )
324  {
325  /* wavelength is very large, convert to microns */
326  sprintf( chLineLbl_v, "%2.2s%2.2s%5.1f%c",
329  t->WLAng/1e4, 'm' );
330  }
331  else if( t->WLAng > 9999. )
332  {
333  sprintf( chLineLbl_v, "%2.2s%2.2s%5.2f%c",
334  elementnames.chElementSym[ t->Hi->nelem -1],
336  t->WLAng/1e4, 'm' );
337  }
338  else if( t->WLAng >= 100. )
339  {
340  sprintf( chLineLbl_v, "%2.2s%2.2s%5i%c",
341  elementnames.chElementSym[ t->Hi->nelem -1],
343  (int)t->WLAng, 'A' );
344  }
345  /* the following two formats should be changed for the
346  * wavelength to get more precision */
347  else if( t->WLAng >= 10. )
348  {
349  sprintf( chLineLbl_v, "%2.2s%2.2s%5.1f%c",
350  elementnames.chElementSym[ t->Hi->nelem -1],
352  t->WLAng, 'A' );
353  }
354  else
355  {
356  sprintf( chLineLbl_v, "%2.2s%2.2s%5.2f%c",
357  elementnames.chElementSym[ t->Hi->nelem -1],
359  t->WLAng, 'A' );
360  }
361  /* make sure that string ends with null character - this should
362  * be redundant */
363  ASSERT( chLineLbl_v[10]==0 );
364  return( chLineLbl_v );
365 }
366 
367 /*RefIndex calculates the index of refraction of air using the line energy in wavenumbers,
368  * used to convert vacuum wavelengths to air wavelengths. */
369 double RefIndex(double EnergyWN )
370 {
371  double RefIndex_v,
372  WaveMic,
373  xl,
374  xn;
375 
376  DEBUG_ENTRY( "RefIndex()" );
377 
378  ASSERT( EnergyWN > 0. );
379 
380  /* the wavelength in microns */
381  WaveMic = 1.e4/EnergyWN;
382 
383  /* only do index of refraction if longward of 2000A */
384  if( WaveMic > 0.2 )
385  {
386  /* longward of 2000A
387  * xl is 1/WaveMic^2 */
388  xl = 1.0/WaveMic/WaveMic;
389  /* use a formula from
390  *>>refer air index refraction Allen, C.W. 1973, Astrophysical quantities,
391  *>>refercon 3rd Edition (AQ), p.124 */
392  xn = 255.4/(41. - xl);
393  xn += 29498.1/(146. - xl);
394  xn += 64.328;
395  RefIndex_v = xn/1.e6 + 1.;
396  }
397  else
398  {
399  RefIndex_v = 1.;
400  }
401  ASSERT( RefIndex_v > 0. );
402  return( RefIndex_v );
403 }
404 
405 /*PutCS enter a collision strength into an individual line vector */
406 void PutCS(double cs,
407  transition * t)
408 {
409 
410  DEBUG_ENTRY( "PutCS()" );
411 
412  /* collision strength must not be negative, had been test for being positive,
413  * but called with zero - did not check why 98 jul 5 */
414  ASSERT( cs >= 0. );
415 
416  t->Coll.cs = (realnum)cs;
417  return;
418 }
419 
420 /*WavlenErrorGet - given the real wavelength in A for a line
421  * routine will find the error expected between the real
422  * wavelength and the wavelength printed in the output, with 4 sig figs,
423  * function returns difference between exact and 4 sig fig wl, so
424  * we have found correct line is fabs(d wl) < return */
426 {
427  double a;
428  realnum errorwave;
429 
430  DEBUG_ENTRY( "WavlenErrorGet()" );
431 
432  ASSERT( LineSave.sig_figs <= 6 );
433 
434  if( wavelength > 0. )
435  {
436  /* normal case, positive (non zero) wavelength */
437  a = log10( wavelength+FLT_EPSILON);
438  a = floor(a);
439  }
440  else
441  {
442  /* might be called with wl of zero, this is that case */
443  /* errorwave = 1e-4f; */
444  a = 0.;
445  }
446 
447  errorwave = 5.f * (realnum)pow( 10., a - (double)LineSave.sig_figs );
448  return errorwave;
449 }
450 
451 /*linadd enter lines into the line storage array, called once per zone for each line*/
452 void linadd(
453  double xInten, /* xInten - local emissivity per unit vol, no fill fac */
454  realnum wavelength, /* realnum wavelength */
455  const char *chLab,/* string label for ion */
456  char chInfo, /* character type of entry for line - given below */
457  /* 'c' cooling, 'h' heating, 'i' info only, 'r' recom line */
458  const char *chComment )
459 {
460  DEBUG_ENTRY( "linadd()" );
461 
462  /* main routine to actually enter lines into the line storage array
463  * called at top level within routine lines
464  * called series of times in routine PutLine for lines transferred
465  */
466 
467  /* three values, -1 is just counting, 0 if init, 1 for calculation */
468  if( LineSave.ipass > 0 )
469  {
470  /* not first pass, sum lines only
471  * total emission from vol */
472  LineSv[LineSave.nsum].sumlin[0] += xInten*radius.dVeff;
473  /* local emissivity in line */
474  LineSv[LineSave.nsum].emslin[0] = xInten;
475  /* no need to increment or set [1] version since this is called with no continuum
476  * index, don't know what to do */
477  if( wavelength > 0 )
478  {
479  /* only put informational lines, like "Q(H) 4861", in this stack
480  * many continua have a wavelength of zero and are proper intensities,
481  * it would be wrong to predict their transferred intensity */
484  }
485  }
486 
487  else if( LineSave.ipass == 0 )
488  {
489  /* first call to stuff lines in array, confirm that label is one of
490  * the four correct ones */
491  /* >>chng 04 apr 24, last two had been != so this test never really happened; PvH */
492  ASSERT( (chInfo == 'c') || (chInfo == 'h') || (chInfo == 'i') || (chInfo == 'r' ) );
493 
494  /* then save it into array */
495  LineSv[LineSave.nsum].chSumTyp = (char)chInfo;
496 
497  /* number of lines ok, set parameters for first pass */
498  LineSv[LineSave.nsum].sumlin[0] = 0.;
499  LineSv[LineSave.nsum].emslin[0] = 0.;
501  LineSv[LineSave.nsum].sumlin[1] = 0.;
502  LineSv[LineSave.nsum].emslin[1] = 0.;
503  LineSv[LineSave.nsum].chComment = chComment;
504  /* check that null is correct, string overruns have
505  * been a problem in the past */
506  ASSERT( strlen( chLab )<5 );
507  strcpy( LineSv[LineSave.nsum].chALab, chLab );
508  }
509 
510  /* increment the line counter */
511  ++LineSave.nsum;
512 
513  /* routine can be called with negative LineSave.ipass, in this case
514  * we are just counting number of lines for current setup */
515  return;
516 }
517 
518 
519 /* this is energy in Ryd as set by call to PntForLine */
520 static double EnergyRyd;
521 
522 /*emergent_line find absorption due to continuous opacity */
524  /* emissivity [erg cm-3 s-1] in inward direction */
525  double emissivity_in ,
526  /* emissivity [erg cm-3 s-1] in outward direction */
527  double emissivity_out ,
528  /* array index for continuum frequency */
529  long int ipCont )
530 {
531 
532  double emergent_in , emergent_out;
533  long int i = ipCont-1;
534 
535  DEBUG_ENTRY( "emergent_line()" );
536 
537  ASSERT( i >= 0 && i < rfield.nupper-1 );
538 
539  /* do nothing if first iteration since we do not know the outward-looking
540  * optical depths. In version C07.02.00 we assumed an infinite optical
541  * depth in the outward direction, which would be appropriate for a
542  * HII region on the surface of a molecular cloud. This converged onto
543  * the correct solution in later iterations, but on the first iteration
544  * this underestimated total emission if the infinite cloud were not
545  * present. With C07.02.01 we make no assuptions about what is in the
546  * outward direction and simply use the local emission.
547  * Behavior is unchanged on later iterations */
548  if( iteration == 1 )
549  {
550  /* first iteration - do not know outer optical depths so assume very large
551  * optical depths */
552  emergent_in = emissivity_in*opac.E2TauAbsFace[i];
553  emergent_out = emissivity_out;
554  }
555  else
556  {
557  if( geometry.lgSphere )
558  {
559  /* second or later iteration in closed or spherical geometry */
560  /* inwardly directed emission must get to central hole then across entire
561  * far side of shell */
562  emergent_in = emissivity_in * opac.E2TauAbsFace[i] *opac.E2TauAbsTotal[i];
563 
564  /* E2 is outwardly directed emission to get to outer edge of cloud */
565  emergent_out = emissivity_out * opac.E2TauAbsOut[i];
566  }
567  else
568  {
569  /* open geometry in second or later iteration, outer optical depths are known
570  * this is light emitted into the outer direction and backscattered
571  * into the inner */
572  double reflected = emissivity_out * opac.albedo[i] * (1.-opac.E2TauAbsOut[i]);
573  /* E2 is to get to central hole */
574  emergent_in = (emissivity_in + reflected) * opac.E2TauAbsFace[i];
575  /* E2 is to get to outer edge */
576  emergent_out = (emissivity_out - reflected) * opac.E2TauAbsOut[i];
577  }
578  }
579  /* return the net emissivity times a vol element */
580  return( (emergent_in + emergent_out) );
581 }
582 
583 /*lindst add line with destruction and outward */
584 void lindst(
585  double xInten,
587  const char *chLab,
588  long int ipnt,
589  char chInfo,
590  bool lgOutToo,
591  const char *chComment )
592 {
593  double saveemis;
594  DEBUG_ENTRY( "lindst()" );
595 
596  /* >>chng 06 feb 08, add test on xInten positive, no need to evaluate
597  * for majority of zero */
598  if( LineSave.ipass > 0 && xInten > 0. )
599  {
600  /* LineSave.ipass > 0, integration across simulation, sum lines only
601  * emissivity, emission per unit vol, for this zone */
602  LineSv[LineSave.nsum].emslin[0] = xInten;
603  /* integrated intensity or luminosity, the emissivity times the volume */
604  LineSv[LineSave.nsum].sumlin[0] += xInten*radius.dVeff;
605 
606  /* add line to outward beam
607  * there are lots of lines that are sums of other lines, or
608  * just for info of some sort. These have flag lgOutToo false.
609  * Note that the EnergyRyd variable only has a rational
610  * value if PntForLine was called just before this routine - in all
611  * cases where this did not happen the flag is false. */
612  if( lgOutToo )
613  {
614  rfield.outlin[ipnt-1] += (realnum)(xInten/(rfield.anu[ipnt-1]*EN1RYD)*
615  radius.dVolOutwrd*opac.ExpZone[ipnt-1]);
616  rfield.reflin[ipnt-1] += (realnum)(xInten/(rfield.anu[ipnt-1]*EN1RYD)*
618  }
619  /* main production pass, sum lines only */
620  if( ipnt <= rfield.nflux )
621  {
622  /* emergent_line accounts for destruction by absorption outside
623  * the line-forming region */
624  saveemis = emergent_line(
625  xInten*rt.fracin , xInten*(1.-rt.fracin) , ipnt );
626  LineSv[LineSave.nsum].emslin[1] = saveemis;
627  LineSv[LineSave.nsum].sumlin[1] += saveemis*radius.dVeff;
628  }
629  }
630  else if( LineSave.ipass == 0 )
631  {
632  ASSERT( ipnt > 0 );
633 
634  /* first call to stuff lines in array, confirm that label is one of
635  * the four correct ones */
636  /* >>chng 04 apr 24, last two had been != so this test never really happened; PvH */
637  ASSERT( (chInfo == 'c') || (chInfo == 'h') || (chInfo == 'i') || (chInfo == 'r' ) );
638  LineSv[LineSave.nsum].chSumTyp = (char)chInfo;
639  /* number of lines ok, set parameters for first pass */
640  LineSv[LineSave.nsum].sumlin[0] = 0.;
641  LineSv[LineSave.nsum].emslin[0] = 0.;
643  LineSv[LineSave.nsum].emslin[1] = 0.;
644  LineSv[LineSave.nsum].sumlin[1] = 0.;
645  LineSv[LineSave.nsum].chComment = chComment;
646 
647  /* check that null is correct, string overruns have
648  * been a problem in the past */
649  ASSERT( strlen( chLab )<5 );
650  strcpy( LineSv[LineSave.nsum].chALab, chLab );
651  }
652 
653  /* increment the line counter */
654  ++LineSave.nsum;
655 
656  EnergyRyd = 0.;
657  return;
658 }
659 
660 /*PntForLine generate pointer for forbidden line */
662  /* wavelength of transition in Angstroms */
663  double wavelength,
664  /* label for this line */
665  const char *chLabel,
666  /* this is array index on the f, not c scale,
667  * for the continuum cell holding the line */
668  long int *ipnt)
669 {
670  /*
671  * maximum number of forbidden lines - this is a good bet since
672  * new lines do not go into this group, and lines are slowly
673  * moving to level 1
674  */
675 # define MAXFORLIN 1000
676  static long int ipForLin[MAXFORLIN]={0};
677 
678  /* number of forbidden lines entered into continuum array */
679  static long int nForLin;
680 
681  DEBUG_ENTRY( "PntForLine()" );
682 
683  /* must be 0 or greater */
684  ASSERT( wavelength >= 0. );
685 
686  if( wavelength == 0. )
687  {
688  /* zero is special flag to initialize */
689  nForLin = 0;
690  /* ipLineEnergy will only put in line label if nothing already there */
691  EnergyRyd = 0.;
692  }
693  else
694  {
695 
696  if( LineSave.ipass > 0 )
697  {
698  /* not first pass, sum lines only */
699  *ipnt = ipForLin[nForLin];
700  }
701  else if( LineSave.ipass == 0 )
702  {
703  /* check if number of lines in arrays exceeded */
704  if( nForLin >= MAXFORLIN )
705  {
706  fprintf( ioQQQ, "PROBLEM %5ld lines is too many for PntForLine.\n",
707  nForLin );
708  fprintf( ioQQQ, " Increase the value of maxForLine everywhere in the code.\n" );
709  cdEXIT(EXIT_FAILURE);
710  }
711 
712  /* ipLineEnergy will only put in line label if nothing already there */
714  ipForLin[nForLin] = ipLineEnergy(EnergyRyd,chLabel , 0);
715  *ipnt = ipForLin[nForLin];
716  }
717  else
718  {
719  /* this is case where we are only counting lines */
720  *ipnt = 0;
721  }
722  ++nForLin;
723  }
724  return;
725 }
726 
728 
729 /*PutLine enter local line intensity into the intensity stack for eventual printout */
730 void PutLine(transition * t, const char *chComment, const char *chLabelTemp)
731 {
732  char chLabel[5];
733  realnum wl;
734  double xIntensity,
735  other,
736  xIntensity_in;
737 
738  DEBUG_ENTRY( "PutLine()" );
739 
740  /* routine to use line array data to generate input
741  * for emission line array */
742  ASSERT( t->ipCont > 0. );
743 
744  strcpy( chLabel, chLabelTemp );
745  chLabel[4] = '\0';
746 
747  /* if ipass=0 then we must generate label info since first pass
748  * gt.0 then only need line intensity data */
749  if( LineSave.ipass == 0 )
750  {
751  /* save wavelength */
752  wl = t->WLAng;
753 
754  xIntensity = 0.;
755  }
756  else
757  {
758  /* both the counting and integrating modes comes here */
759  /*>>chng 06 mar 10, these are not actually used so set to impossible values */
760  /*strcpy( chLabel, " " );*/
761  chLabel[0] = '\n';
762  wl = 0.;
763  /* total line intensity or luminosity
764  * these may not be defined in initial calls so define here */
765  xIntensity = t->Emis->xIntensity + ExtraInten;
766  }
767  /* initial counting case, where ipass == -1, just ignored above, call linadd below */
768 
769  /* ExtraInten is option that allows extra intensity (i.e., recomb)
770  * to be added to this line with Call PutExtra( exta )
771  * in main lines this extra
772  * contribution must be identified explicitly */
773  ExtraInten = 0.;
774  /*linadd(xIntensity,wl,chLabel,'i');*/
775  /*lindst add line with destruction and outward */
776  rt.fracin = t->Emis->FracInwd;
777  lindst(xIntensity,
778  wl,
779  chLabel,
780  t->ipCont,
781  /* this is information only - has been counted in cooling already */
782  'i',
783  /* do not add to outward beam, also done separately */
784  false,
785  chComment);
786  rt.fracin = 0.5;
787 
788  /* inward part of line - do not move this away from previous lines
789  * since xIntensity is used here */
790  xIntensity_in = xIntensity*t->Emis->FracInwd;
791  linadd(xIntensity_in,wl,"Inwd",'i',chComment);
792 
793  /* cooling part of line */
794  other = t->Coll.cool;
795  linadd(other,wl,"Coll",'i',chComment);
796 
797  /* fluorescent excited part of line */
798  other = t->Emis->PopOpc * t->Emis->pump * (1.-t->Emis->ColOvTot) * t->EnergyErg;
799  linadd(other,wl,"Pump",'i',chComment);
800 
801  /* heating part of line */
802  other = t->Coll.heat;
803  linadd(other,wl,"Heat",'i',chComment);
804  return;
805 }
806 
807 /*PutLine enter local line intensity into the intensity stack for eventual printout */
809  const char *chComment)
810 {
811  char chLabel[5];
812  realnum wl;
813  double xIntensity,
814  other,
815  xIntensity_in;
816 
817  DEBUG_ENTRY( "PutLine()" );
818 
819  /* routine to use line array data to generate input
820  * for emission line array */
821  ASSERT( t->ipCont > 0. );
822 
823  /* if ipass=0 then we must generate label info since first pass
824  * gt.0 then only need line intensity data */
825  if( LineSave.ipass == 0 )
826  {
827  /* these variables not used by linadd if ipass ne 0 */
828  /* chIonLbl is function that generates a null terminated 4 char string, of form "C 2" */
829  chIonLbl(chLabel, t);
830 
831  /* save wavelength */
832  wl = t->WLAng;
833 
834  xIntensity = 0.;
835  }
836  else
837  {
838  /* both the counting and integrating modes comes here */
839  /*>>chng 06 mar 10, these are not actually used so set to impossible values */
840  /*strcpy( chLabel, " " );*/
841  chLabel[0] = '\n';
842  wl = 0.;
843  /* total line intensity or luminosity
844  * these may not be defined in initial calls so define here */
845  xIntensity = t->Emis->xIntensity + ExtraInten;
846  }
847  /* initial counting case, where ipass == -1, just ignored above, call linadd below */
848 
849  /* ExtraInten is option that allows extra intensity (i.e., recomb)
850  * to be added to this line with Call PutExtra( exta )
851  * in main lines this extra
852  * contribution must be identified explicitly */
853  ExtraInten = 0.;
854  /*linadd(xIntensity,wl,chLabel,'i');*/
855  /*lindst add line with destruction and outward */
856  rt.fracin = t->Emis->FracInwd;
857  lindst(xIntensity,
858  wl,
859  chLabel,
860  t->ipCont,
861  /* this is information only - has been counted in cooling already */
862  'i',
863  /* do not add to outward beam, also done separately */
864  false,
865  chComment);
866  rt.fracin = 0.5;
867 
868  /* inward part of line - do not move this away from previous lines
869  * since xIntensity is used here */
870  xIntensity_in = xIntensity*t->Emis->FracInwd;
871  linadd(xIntensity_in,wl,"Inwd",'i',chComment);
872 
873  /* cooling part of line */
874  other = t->Coll.cool;
875  linadd(other,wl,"Coll",'i',chComment);
876 
877  /* fluorescent excited part of line */
878  other = t->Emis->PopOpc * t->Emis->pump * (1.-t->Emis->ColOvTot) * t->EnergyErg;
879  linadd(other,wl,"Pump",'i',chComment);
880 
881  /* heating part of line */
882  other = t->Coll.heat;
883  linadd(other,wl,"Heat",'i',chComment);
884  return;
885 }
886 
887 /* ==================================================================== */
888 /*PutExtra enter and 'extra' intensity source for some line */
889 void PutExtra(double Extra)
890 {
891 
892  DEBUG_ENTRY( "PutExtra()" );
893 
894  ExtraInten = (realnum)Extra;
895  return;
896 }
897 
899 {
900 
901  DEBUG_ENTRY( "TransitionJunk()" );
902 
903  /* wavelength, usually in A, used for printout */
904  t->WLAng = -FLT_MAX;
905 
906  /* transition energy in degrees kelvin*/
907  t->EnergyK = -FLT_MAX;
908  /* transition energy in ergs */
909  t->EnergyErg = -FLT_MAX;
910  /* transition energy in wavenumbers */
911  t->EnergyWN = -FLT_MAX;
912 
913  /* array offset for radiative transition within continuum array
914  * is negative if transition is non-radiative. */
915  t->ipCont = -10000;
916 
917  CollisionJunk( &t->Coll );
918 
919  /* set these equal to NULL first. That will cause the code to crash if
920  * the variables are ever used before being deliberately set. */
921  t->Emis = &DummyEmis;
922 
923  t->Lo = NULL;
924  t->Hi = NULL;
925  return;
926 }
927 
928 
929 /*EmLineJunk set all elements of transition struc to dangerous values */
931 {
932 
933  DEBUG_ENTRY( "EmLineJunk()" );
934 
935  /* optical depth in continuum to ill face */
936  t->TauCon = -FLT_MAX;
937 
938  /* inward and total line optical depths */
939  t->TauIn = -FLT_MAX;
940  t->TauTot = -FLT_MAX;
941 
942  /* type of redistribution function, */
943  t->iRedisFun = INT_MIN;
944 
945  /* array offset for line within fine opacity */
946  t->ipFine = -10000;
947 
948  /* inward fraction */
949  t->FracInwd = -FLT_MAX;
950 
951  /* continuum pumping rate */
952  t->pump = -FLT_MAX;
953 
954  /* line intensity */
955  t->xIntensity = -FLT_MAX;
956 
957  /* number of photons emitted per sec in the line */
958  t->phots = -FLT_MAX;
959 
960  /* gf value */
961  t->gf = -FLT_MAX;
962 
963  /* escape and destruction probs */
964  t->Pesc = -FLT_MAX;
965  t->Pdest = -FLT_MAX;
966  t->Pelec_esc = -FLT_MAX;
967 
968  /* damping constant, and number related to it */
969  t->dampXvel = -FLT_MAX;
970  t->damp = -FLT_MAX;
971 
972  /* ratio of collisional to radiative excitation*/
973  t->ColOvTot = -FLT_MAX;
974 
975  /* line opacity */
976  t->opacity = -FLT_MAX;
977 
978  t->PopOpc = -FLT_MAX;
979 
980  /* transition prob, Einstein A upper to lower */
981  t->Aul = -FLT_MAX;
982 
983  /* ots rate */
984  t->ots = -FLT_MAX;
985  return;
986 }
987 
988 /*CollisionJunk set all elements of transition struc to dangerous values */
990 {
991 
992  DEBUG_ENTRY( "CollisionJunk()" );
993 
995  t->ColUL = -FLT_MAX;
996 
997  /* Coll->cooling and Coll->heating due to collisional excitation */
998  t->cool = -FLT_MAX;
999  t->heat = -FLT_MAX;
1000 
1001  /* collision strengths for transition */
1002  t->cs = -FLT_MAX;
1003  for( long i=0; i<ipNCOLLIDER; i++ )
1004  t->csi[i] = -FLT_MAX;
1005  return;
1006 }
1007 
1008 /*StateJunk set all elements of transition struc to dangerous values */
1010 {
1011 
1012  DEBUG_ENTRY( "StateJunk()" );
1013 
1014  /* t->chLabel = ''; */
1015 
1017  t->g = -FLT_MAX;
1018 
1020  t->Pop = -FLT_MAX;
1021 
1023  t->IonStg = -10000;
1024 
1026  t->nelem = -10000;
1027  return;
1028 }
1029 
1030 /*TransitionZero zeros out TransitionZero */
1032 {
1033 
1034  DEBUG_ENTRY( "TransitionZero()" );
1035 
1036  CollisionZero( &t->Coll );
1037 
1038  StateZero( t->Lo );
1039  StateZero( t->Hi );
1040  EmLineZero( t->Emis );
1041 
1042  return;
1043 }
1044 
1045 /*EmLineZero zeros out the emission line structure */
1047 {
1048 
1049  DEBUG_ENTRY( "EmLineZero()" );
1050 
1051  /* total optical depth in all overlapping lines to illuminated face,
1052  * used for pumping */
1053  t->TauCon = opac.taumin;
1054 
1055  /* inward and total line optical depths */
1056  /* >>chng 03 feb 14, from 0 to opac.taumin */
1057  t->TauIn = opac.taumin;
1058 
1059  /* total optical depths */
1060  t->TauTot = 1e20f;
1061 
1062  /* inward fraction */
1063  /* >>chng 03 feb 14, from 0 to 1 */
1064  t->FracInwd = 1.;
1065 
1066  /* continuum pumping rate */
1067  t->pump = 0.;
1068 
1069  /* line intensity */
1070  t->xIntensity = 0.;
1071 
1072  /* number of photons emitted per sec in the line */
1073  t->phots = 0.;
1074 
1075  /* escape and destruction probs */
1076  /* >>chng 03 feb 14, change from 0 to 1 */
1077  t->Pesc = 1.;
1078  t->Pdest = 0.;
1079  t->Pelec_esc = 0.;
1080 
1081  /* ratio of collisional to radiative excitation*/
1082  t->ColOvTot = 0.;
1083 
1084  /* pop that enters net opacity */
1085  t->PopOpc = 0.;
1086 
1087  /* ots rate */
1088  t->ots = 0.;
1089  return;
1090 }
1091 
1092 /*CollisionZero zeros out the structure */
1094 {
1095 
1096  DEBUG_ENTRY( "CollisionZero()" );
1097 
1098  /* Coll->cooling and Coll->heating due to collisional excitation */
1099  t->cool = 0.;
1100  t->heat = 0.;
1101  return;
1102 }
1103 
1104 /*StateZero zeros out the structure */
1106 {
1107 
1108  DEBUG_ENTRY( "StateZero()" );
1109 
1111  t->Pop = 0.;
1112  return;
1113 }
1114 
1115 /*LineConvRate2CS convert down coll rate back into electron cs in case other parts of code need this for reference */
1117 {
1118 
1119  DEBUG_ENTRY( "LineConvRate2CS()" );
1120 
1121  /* return is collision strength, convert from collision rate from
1122  * upper to lower, this assumes pure electron collisions, but that will
1123  * also be assumed by anything that uses cs, for self-consistency */
1124  t->Coll.cs = rate * t->Hi->g / (realnum)dense.cdsqte;
1125 
1126  /* change assert to non-negative - there can be cases (Iin H2) where cs has
1127  * underflowed to 0 on some platforms */
1128  ASSERT( t->Coll.cs >= 0. );
1129  return;
1130 }
1131 
1132 /*ConvRate2CS convert down coll rate back into electron cs in case other parts of code need this for reference */
1133 double ConvRate2CS( realnum gHi , realnum rate )
1134 {
1135 
1136  double cs;
1137 
1138  DEBUG_ENTRY( "ConvRate2CS()" );
1139 
1140  /* return is collision strength, convert from collision rate from
1141  * upper to lower, this assumes pure electron collisions, but that will
1142  * also be assumed by anything that uses cs, for self-consistency */
1143  cs = rate * gHi / dense.cdsqte;
1144 
1145  /* change assert to non-negative - there can be cases (Iin H2) where cs has
1146  * underflowed to 0 on some platforms */
1147  ASSERT( cs >= 0. );
1148  return cs;
1149 }
1150 
1151 /* returns true is we have not overrun optical depth scale */
1153 {
1154  bool lgGoodTau;
1155 
1156  DEBUG_ENTRY( "lgTauGood()" );
1157 
1158  if( (t->Emis->TauTot*0.9 - t->Emis->TauIn < 0. && t->Emis->TauIn > 0.) &&
1159  ! fp_equal( t->Emis->TauTot, opac.taumin ) )
1160  {
1161  /* do not do anything if we have overrun the scale, */
1162  lgGoodTau = false;
1163  }
1164  else
1165  {
1166  lgGoodTau = true;
1167  }
1168  return lgGoodTau;
1169 }
1170 
1171 /*gbar0 compute g-bar gaunt factor for neutrals */
1172 STATIC void gbar0(double ex,
1173  realnum *g)
1174 {
1175  double a,
1176  b,
1177  c,
1178  d,
1179  y;
1180 
1181  DEBUG_ENTRY( "gbar0()" );
1182 
1183  /* written by Dima Verner
1184  *
1185  * Calculation of the effective Gaunt-factor by use of
1186  * >>refer gbar cs Van Regemorter, H., 1962, ApJ 136, 906
1187  * fits for neutrals
1188  * Input parameters:
1189  * ex - energy ryd - now K
1190  * t - temperature in K
1191  * Output parameter:
1192  * g - effective Gaunt factor
1193  * */
1194 
1195  /* y = ex*157813.7/te */
1196  y = ex/phycon.te;
1197  if( y < 0.01 )
1198  {
1199  *g = (realnum)(0.29*(log(1.0+1.0/y) - 0.4/POW2(y + 1.0))/exp(y));
1200  }
1201  else
1202  {
1203  if( y > 10.0 )
1204  {
1205  *g = (realnum)(0.066/sqrt(y));
1206  }
1207  else
1208  {
1209  a = 1.5819068e-02;
1210  b = 1.3018207e00;
1211  c = 2.6896230e-03;
1212  d = 2.5486007e00;
1213  d = log(y/c)/d;
1214  *g = (realnum)(a + b*exp(-0.5*POW2(d)));
1215  }
1216  }
1217  return;
1218 }
1219 
1220 /*gbar1 compute g-bar collision strength using Mewe approximations */
1221 STATIC void gbar1(double ex,
1222  realnum *g)
1223 {
1224  double y;
1225 
1226  DEBUG_ENTRY( "gbar1()" );
1227 
1228  /* *written by Dima Verner
1229  *
1230  * Calculation of the effective Gaunt-factor by use of
1231  * >>refer gbar cs Mewe,R., 1972, A&A 20, 215
1232  * fits for permitted transitions in ions MgII, CaII, FeII (delta n = 0)
1233  * Input parameters:
1234  * ex - excitation energy in Ryd - now K
1235  * te - temperature in K
1236  * Output parameter:
1237  * g - effective Gaunt factor
1238  */
1239 
1240  /* y = ex*157813.7/te */
1241  y = ex/phycon.te;
1242  *g = (realnum)(0.6 + 0.28*(log(1.0+1.0/y) - 0.4/POW2(y + 1.0)));
1243  return;
1244 }
1245 
1246 /*MakeCS compute collision strength by g-bar approximations */
1248 {
1249  long int ion;
1250  double Abun,
1251  cs;
1252  realnum
1253  gbar;
1254 
1255  DEBUG_ENTRY( "MakeCS()" );
1256 
1257  /* routine to get cs from various approximations */
1258 
1259  /* check if abundance greater than 0 */
1260  ion = t->Hi->IonStg;
1261 
1262  Abun = dense.xIonDense[ t->Hi->nelem -1 ][ ion-1 ];
1263  if( Abun <= 0. )
1264  {
1265  gbar = 1.;
1266  }
1267  else
1268  {
1269  /* check if neutral or ion */
1270  if( ion == 1 )
1271  {
1272  /* neutral - compute gbar for eventual cs */
1273  gbar0(t->EnergyK,&gbar);
1274  }
1275  else
1276  {
1277  /* ion - compute gbar for eventual cs */
1278  gbar1(t->EnergyK,&gbar);
1279  }
1280  }
1281 
1282  /* above was g-bar, convert to cs */
1283  cs = gbar*(14.5104/WAVNRYD)*t->Emis->gf/t->EnergyWN;
1284 
1285  /* stuff the cs in place */
1286  t->Coll.cs = (realnum)cs;
1287  return;
1288 }
1289 
1290 /*totlin sum total intensity of cooling, recombination, or intensity lines */
1291 double totlin(
1292  /* chInfor is 1 char,
1293  'i' information,
1294  'r' recombination or
1295  'c' collision */
1296  int chInfo)
1297 {
1298  long int i;
1299  double totlin_v;
1300 
1301  DEBUG_ENTRY( "totlin()" );
1302 
1303  /* routine goes through set of entered line
1304  * intensities and picks out those which have
1305  * types agreeing with chInfo. Valid types are
1306  * 'c', 'r', and 'i'
1307  *begin sanity check */
1308  if( (chInfo != 'i' && chInfo != 'r') && chInfo != 'c' )
1309  {
1310  fprintf( ioQQQ, " TOTLIN does not understand chInfo=%c\n",
1311  chInfo );
1312  cdEXIT(EXIT_FAILURE);
1313  }
1314  /*end sanity check */
1315 
1316  /* now find sum of lines of type chInfo */
1317  totlin_v = 0.;
1318  for( i=0; i < LineSave.nsum; i++ )
1319  {
1320  if( LineSv[i].chSumTyp == chInfo )
1321  {
1322  totlin_v += LineSv[i].sumlin[0];
1323  }
1324  }
1325  return( totlin_v );
1326 }
1327 
1328 
1329 /*FndLineHt search through line heat arrays to find the strongest heat source */
1330 void FndLineHt(long int *level,
1331  /* this is the index of the strongest line in the array on the c scale */
1332  long int *ipStrong,
1333  double *Strong)
1334 {
1335  long int i;
1336 
1337  DEBUG_ENTRY( "FndLineHt()" );
1338 
1339  *Strong = 0.;
1340  *level = 0;
1341 
1342  /* do the level 1 lines, 0 is dummy line, <=nLevel1 is correct for c scale */
1343  for( i=1; i <= nLevel1; i++ )
1344  {
1345  /* check if a line was the major heat agent */
1346  if( TauLines[i].Coll.heat > *Strong )
1347  {
1348  *ipStrong = i;
1349  *level = 1;
1350  *Strong = TauLines[i].Coll.heat;
1351  }
1352  }
1353 
1354  /* now do the level 2 lines */
1355  for( i=0; i < nWindLine; i++ )
1356  {
1357  if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
1358  {
1359  /* check if a line was the major heat agent */
1360  if( TauLine2[i].Coll.heat > *Strong )
1361  {
1362  *ipStrong = i;
1363  *level = 2;
1364  *Strong = TauLine2[i].Coll.heat;
1365  }
1366  }
1367  }
1368 
1369  /* now do co carbon monoxide lines */
1370  for( i=0; i < nCORotate; i++ )
1371  {
1372  /* check if a line was the major heat agent */
1373  if( C12O16Rotate[i].Coll.heat > *Strong )
1374  {
1375  *ipStrong = i;
1376  *level = 3;
1377  *Strong = C12O16Rotate[i].Coll.heat;
1378  }
1379  }
1380  for( i=0; i < nCORotate; i++ )
1381  {
1382  /* check if a line was the major heat agent */
1383  if( C13O16Rotate[i].Coll.heat > *Strong )
1384  {
1385  *ipStrong = i;
1386  *level = 4;
1387  *Strong = C13O16Rotate[i].Coll.heat;
1388  }
1389  }
1390 
1391  /* now do the hyperfine structure lines */
1392  for( i=0; i < nHFLines; i++ )
1393  {
1394  /* check if a line was the major heat agent */
1395  if( HFLines[i].Coll.heat > *Strong )
1396  {
1397  *ipStrong = i;
1398  *level = 5;
1399  *Strong = HFLines[i].Coll.heat;
1400  }
1401  }
1402 
1403  /*Chianti & Leiden Atoms & Molecules */
1404  for( i=0; i <linesAdded2; i++)
1405  {
1406  /* check if a line was the major heat agent */
1407  if( atmolEmis[i].tran->Coll.heat > *Strong )
1408  {
1409  *ipStrong = i;
1410  *level = 6;
1411  *Strong = atmolEmis[i].tran->Coll.heat;
1412  }
1413  }
1414  return;
1415 }
1416 
1418 {
1419  DEBUG_ENTRY( "AddState2Stack()" );
1420 
1421  ASSERT( !lgStatesAdded );
1422 
1423  currentState = new quantumState;
1424 
1426 
1427  if( statesAdded == 0 )
1428  {
1430  GenericStates->next = NULL;
1432  }
1433  else
1434  {
1438  }
1439 
1440  statesAdded++;
1441 
1442  return currentState;
1443 }
1444 
1445 emission *AddLine2Stack( bool lgRadiativeTrans )
1446 {
1447  DEBUG_ENTRY( "AddLine2Stack()" );
1448 
1449  if( !lgRadiativeTrans )
1450  {
1451  return &DummyEmis;
1452  }
1453  else
1454  {
1455  ASSERT( lgLinesAdded == false );
1456 
1457  currentLine = new emission;
1458 
1460 
1461  if( linesAdded == 0 )
1462  {
1464  GenericLines->next = NULL;
1466  }
1467  else
1468  {
1469  /*
1470  \todo 2 Does doing EmLineZero here defeat the purpose of EmLineJunk?
1471  * maybe we should pass full set of Emis components, fill everything in
1472  * here, and THEN use EmLineZero? */
1474 
1476  lastLine = lastLine->next;
1477  }
1478 
1479  linesAdded++;
1480  return currentLine;
1481  }
1482 }

Generated for cloudy by doxygen 1.8.3.1