cloudy trunk
Loading...
Searching...
No Matches
lines_service.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/*GetGF convert Einstein A into oscillator strength */
4/*abscf convert gf into absorption coefficient */
5/*RefIndex calculates the index of refraction of air using the line energy in wavenumbers,
6 * used to convert vacuum wavelengths to air wavelengths. */
7/*eina convert a gf into an Einstein A */
8/*WavlenErrorGet - find difference between two wavelengths */
9/*linadd enter lines into the line storage array, called once per zone */
10/*lindst add local line intensity to line luminosity stack */
11/*PntForLine generate pointer for forbidden line */
12/*totlin sum total intensity of cooling, recombination, or intensity lines */
13/*FndLineHt search through line heat arrays to find the strongest heat source */
14/*ConvRate2CS convert down coll rate back into electron cs in case other parts of code need this for reference */
15#include "cddefines.h"
16#include "lines_service.h"
17#include "dense.h"
18#include "geometry.h"
19#include "hydrogenic.h"
20#include "ipoint.h"
21#include "iso.h"
22#include "lines.h"
23#include "trace.h"
24#include "opacity.h"
25#include "physconst.h"
26#include "radius.h"
27#include "rfield.h"
28#include "rt.h"
29#include "taulines.h"
30#include "thirdparty.h"
31
33{
34 DEBUG_ENTRY( "LineStackCreate()" );
35
36 // set up emission line intensity stack
37 /* there are three types of calls to lines()
38 * ipass = -1, first call, only count number of lines
39 * ipass = 0, second pass, save labels and wavelengths
40 * ipass = 1, integrate intensity*/
41 LineSave.ipass = -1;
42 lines();
43 ASSERT( LineSave.nsum > 0 );
44
45 /* in a grid or MPI run this may not be the first time here,
46 * return old memory and grab new appropriate for this size,
47 * since number of lines to be stored can change */
48 if( LineSv != NULL )
49 free( LineSv );
50 if( LineSvSortWL != NULL )
51 free( LineSvSortWL );
52
53 /* this is the large main line array */
54 LineSv = (LinSv*)MALLOC((unsigned)LineSave.nsum*sizeof( LinSv ) );
55 LineSvSortWL = (LinSv*)MALLOC((unsigned)LineSave.nsum*sizeof( LinSv ));
56
57 /* this will be used as sanity check in future iterations and grid points */
58 LineSave.nsumAllocated = LineSave.nsum;
59
60 /* this is only done on first iteration since will integrate over time */
61 for( long i=0; i < LineSave.nsum; i++ )
62 {
63 for( long j =0; j<4; ++j )
64 LineSv[i].SumLine[j] = 0;
65 }
66
67 /* there are three calls to lines()
68 * first call with ipass = -1 was above, only counted number
69 * of lines and malloced space. This is second call and will do
70 * one-time creation of line labels */
71 LineSave.ipass = 0;
72 lines();
73 /* has to be positive */
74 ASSERT( LineSave.nsum > 0);
75 /* in the future calls to lines will result in integrations */
76 LineSave.ipass = 1;
77
78 if( trace.lgTrace )
79 fprintf( ioQQQ, "%7ld lines printed in main line array\n",
80 LineSave.nsum );
81}
82
83/*eina convert a gf into an Einstein A */
84double eina(double gf,
85 double enercm,
86 double gup)
87{
88 double eina_v;
89
90 DEBUG_ENTRY( "eina()" );
91
92 /* derive the transition prob, given the following
93 * call to function is gf, energy in cm^-1, g_up
94 * gf is product of g and oscillator strength
95 * eina = ( gf / 1.499e-8 ) / (wl/1e4)**2 / gup */
96 eina_v = (gf/gup)*TRANS_PROB_CONST*POW2(enercm);
97 return( eina_v );
98}
99
100/*GetGF convert Einstein A into oscillator strength */
101double GetGF(double trans_prob,
102 double enercm,
103 double gup)
104{
105 double GetGF_v;
106
107 DEBUG_ENTRY( "GetGF()" );
108
109 ASSERT( enercm > 0. );
110 ASSERT( trans_prob > 0. );
111 ASSERT( gup > 0.);
112
113 /* derive the transition prob, given the following
114 * call to function is gf, energy in cm^-1, g_up
115 * gf is product of g and oscillator strength
116 * trans_prob = ( GetGF/gup) / 1.499e-8 / ( 1e4/enercm )**2 */
117 GetGF_v = trans_prob*gup/TRANS_PROB_CONST/POW2(enercm);
118 return( GetGF_v );
119}
120
121/*abscf convert gf into absorption coefficient */
122double abscf(double gf,
123 double enercm,
124 double gl)
125{
126 double abscf_v;
127
128 DEBUG_ENTRY( "abscf()" );
129
130 ASSERT(gl > 0. && enercm > 0. && gf > 0. );
131
132 /* derive line absorption coefficient, given the following:
133 * gf, enercm, g_low
134 * gf is product of g and oscillator strength */
135 abscf_v = 1.4974e-6*(gf/gl)*(1e4/enercm);
136 return( abscf_v );
137}
138
139/*RefIndex calculates the index of refraction of air using the line energy in wavenumbers,
140 * used to convert vacuum wavelengths to air wavelengths. */
141double RefIndex(double EnergyWN )
142{
143 double RefIndex_v,
144 WaveMic,
145 xl,
146 xn;
147
148 DEBUG_ENTRY( "RefIndex()" );
149
150 ASSERT( EnergyWN > 0. );
151
152 /* the wavelength in microns */
153 WaveMic = 1.e4/EnergyWN;
154
155 /* only do index of refraction if longward of 2000A */
156 if( WaveMic > 0.2 )
157 {
158 /* longward of 2000A
159 * xl is 1/WaveMic^2 */
160 xl = 1.0/WaveMic/WaveMic;
161 /* use a formula from
162 *>>refer air index refraction Allen, C.W. 1973, Astrophysical quantities,
163 *>>refercon 3rd Edition (AQ), p.124 */
164 xn = 255.4/(41. - xl);
165 xn += 29498.1/(146. - xl);
166 xn += 64.328;
167 RefIndex_v = xn/1.e6 + 1.;
168 }
169 else
170 {
171 RefIndex_v = 1.;
172 }
173 ASSERT( RefIndex_v >= 1. );
174 return( RefIndex_v );
175}
176
177/*WavlenErrorGet - given the real wavelength in A for a line
178 * routine will find the error expected between the real
179 * wavelength and the wavelength printed in the output, with 4 sig figs,
180 * function returns difference between exact and 4 sig fig wl, so
181 * we have found correct line is fabs(d wl) < return */
183{
184 double a;
185 realnum errorwave;
186
187 DEBUG_ENTRY( "WavlenErrorGet()" );
188
189 ASSERT( LineSave.sig_figs <= 6 );
190
191 if( wavelength > 0. )
192 {
193 /* normal case, positive (non zero) wavelength */
194 a = log10( wavelength+FLT_EPSILON);
195 a = floor(a);
196 }
197 else
198 {
199 /* might be called with wl of zero, this is that case */
200 /* errorwave = 1e-4f; */
201 a = 0.;
202 }
203
204 errorwave = 5.f * (realnum)pow( 10., a - (double)LineSave.sig_figs );
205 return errorwave;
206}
207
208/*linadd enter lines into the line storage array, called once per zone for each line*/
210 double xInten, /* xInten - local emissivity per unit vol, no fill fac */
211 realnum wavelength, /* realnum wavelength */
212 const char *chLab,/* string label for ion */
213 // ipnt offset of line in continuum mesh
214 long int ipnt,
215 char chInfo, /* character type of entry for line - given below */
216 /* 'c' cooling, 'h' heating, 'i' info only, 'r' recom line , 't' transferred line */
217 // *chComment string explaining line
218 const char *chComment,
219 // lgAdd says whether we've come in via linadd (true) or lindst (false)
220 bool lgAdd)
221{
222 DEBUG_ENTRY( "lincom()" );
223
224 /* main routine to actually enter lines into the line storage array
225 * called at top level within routine lines
226 * called series of times in routine PutLine for lines transferred
227 */
228
229 /* three values, -1 is just counting, 0 if init, 1 for calculation */
230 if( LineSave.ipass > 0 )
231 {
232 /* not first pass, sum lines only
233 * total emission from vol */
234 /* LineSave.ipass > 0, integration across simulation, sum lines only
235 * emissivity, emission per unit vol, for this zone */
236 LineSv[LineSave.nsum].SumLine[0] += xInten*radius.dVeffAper;
237 /* local emissivity in line */
238 /* integrated intensity or luminosity, the emissivity times the volume */
239 LineSv[LineSave.nsum].emslin[0] = xInten;
240
241 if (lgAdd)
242 {
243 if (wavelength > 0 && chInfo == 't' )
244 {
245 /* no need to increment or set [1] version since this is called with no continuum
246 * index, don't know what to do */
247 /* only put informational lines, like "Q(H) 4861", in this stack
248 * many continua have a wavelength of zero and are proper intensities,
249 * it would be wrong to predict their transferred intensity */
250 LineSv[LineSave.nsum].emslin[1] = LineSv[LineSave.nsum].emslin[0];
251 LineSv[LineSave.nsum].SumLine[1] = LineSv[LineSave.nsum].SumLine[0];
252 }
253 }
254 else
255 {
256 if ( ipnt <= rfield.nflux && chInfo == 't' )
257 {
258 /* emergent_line accounts for destruction by absorption outside
259 * the line-forming region */
260 const double saveemis = emergent_line(
261 xInten*rt.fracin , xInten*(1.-rt.fracin) , ipnt );
262 LineSv[LineSave.nsum].emslin[1] = saveemis;
263 LineSv[LineSave.nsum].SumLine[1] += saveemis*radius.dVeffAper;
264 }
265 }
266 }
267
268 else if( LineSave.ipass == 0 )
269 {
270 /* first call to stuff lines in array, confirm that label is one of
271 * the four correct ones */
272 ASSERT( (chInfo == 'c') || (chInfo == 'h') || (chInfo == 'i') || (chInfo == 'r' ) || (chInfo == 't') );
273 /* then save it into array */
274 LineSv[LineSave.nsum].chSumTyp = (char)chInfo;
275 LineSv[LineSave.nsum].emslin[0] = 0.;
276 LineSv[LineSave.nsum].emslin[1] = 0.;
277 LineSv[LineSave.nsum].chComment = chComment;
278 /* check that null is correct, string overruns have
279 * been a problem in the past */
280 ASSERT( strlen( chLab )<5 );
281 strcpy( LineSv[LineSave.nsum].chALab, chLab );
282
283 if (lgAdd)
284 {
285 LineSv[LineSave.nsum].wavelength = wavelength;
286 }
287 else
288 {
289 // number of lines ok, set parameters for first pass
290 // negative wavelengh means it is just label, possibly not correct
291 LineSv[LineSave.nsum].wavelength = fabs(wavelength);
292 LineSv[LineSave.nsum].SumLine[0] = 0.;
293 LineSv[LineSave.nsum].SumLine[1] = 0.;
294
295 // check that line wavelength and continuum index agree to some extent
296 // this check cannot be very precise because some lines have
297 // "wavelengths" that are set by common usage rather than the correct
298 // wavelength derived from energy and index of refraction of air
299 ASSERT( ipnt > 0 );
300# ifndef NDEBUG
301 double error = MAX2(0.1*rfield.AnuOrg[ipnt-1] , rfield.widflx[ipnt-1] );
302 ASSERT( wavelength<=0 ||
303 fabs( rfield.AnuOrg[ipnt-1] - RYDLAM / wavelength) < error );
304# endif
305 }
306 }
307
308 /* increment the line counter */
309 ++LineSave.nsum;
310
311 /* routine can be called with negative LineSave.ipass, in this case
312 * we are just counting number of lines for current setup */
313}
314
315/*linadd enter lines into the line storage array, called once per zone for each line*/
317 double xInten, /* xInten - local emissivity per unit vol, no fill fac */
318 realnum wavelength, /* realnum wavelength */
319 const char *chLab,/* string label for ion */
320 char chInfo, /* character type of entry for line - given below */
321 /* 'c' cooling, 'h' heating, 'i' info only, 'r' recom line */
322 const char *chComment )
323{
324 DEBUG_ENTRY( "linadd()" );
325
326 // Values added to get common interface with lindst
327 const long int ipnt = LONG_MAX;
328
329 lincom( xInten, wavelength, chLab, ipnt, chInfo, chComment, true );
330}
331
332
333/*emergent_line find emission from surface of cloud after correcting for
334 * extinction due to continuous opacity for inward & outward directed emission */
336 /* emiemission in inward direction */
337 double emissivity_in ,
338 /* emission in outward direction */
339 double emissivity_out ,
340 /* array index for continuum frequency on fortran scale */
341 long int ipCont )
342{
343
344 double emergent_in , emergent_out;
345 long int i = ipCont-1;
346
347 DEBUG_ENTRY( "emergent_line()" );
348
349 ASSERT( i >= 0 && i < rfield.nupper-1 );
350
351 /* do nothing if first iteration since we do not know the outward-looking
352 * optical depths. In version C07.02.00 we assumed an infinite optical
353 * depth in the outward direction, which would be appropriate for a
354 * HII region on the surface of a molecular cloud. This converged onto
355 * the correct solution in later iterations, but on the first iteration
356 * this underestimated total emission if the infinite cloud were not
357 * present. With C07.02.01 we make no assuptions about what is in the
358 * outward direction and simply use the local emission.
359 * Behavior is unchanged on later iterations */
360 if( iteration == 1 )
361 {
362 /* first iteration - do not know outer optical depths so assume very large
363 * optical depths */
364 emergent_in = emissivity_in*opac.E2TauAbsFace[i];
365 emergent_out = emissivity_out;
366 }
367 else
368 {
369 if( geometry.lgSphere )
370 {
371 /* second or later iteration in closed or spherical geometry */
372 /* inwardly directed emission must get to central hole then across entire
373 * far side of shell */
374 emergent_in = emissivity_in * opac.E2TauAbsFace[i] *opac.E2TauAbsTotal[i];
375
376 /* E2 is outwardly directed emission to get to outer edge of cloud */
377 emergent_out = emissivity_out * opac.E2TauAbsOut[i];
378 }
379 else
380 {
381 /* open geometry in second or later iteration, outer optical depths are known
382 * this is light emitted into the outer direction and backscattered
383 * into the inner */
384 double reflected = emissivity_out * opac.albedo[i] * (1.-opac.E2TauAbsOut[i]);
385 /* E2 is to get to central hole */
386 emergent_in = (emissivity_in + reflected) * opac.E2TauAbsFace[i];
387 /* E2 is to get to outer edge */
388 emergent_out = (emissivity_out - reflected) * opac.E2TauAbsOut[i];
389 }
390 }
391 /* return the net emission that makes it to the surface */
392 return( emergent_in + emergent_out );
393}
394
395/* outline_base - calls outline_base_bin after deciding whether to add impulse or resolved line */
396void outline_base(double dampXvel, double damp, bool lgTransStackLine, long int ip, double phots, realnum inwd,
397 double nonScatteredFraction)
398{
399 DEBUG_ENTRY( "outline_base()" );
400
401 static const bool DO_PROFILE = false;
402
403 if( !DO_PROFILE || !rfield.lgDoLineTrans )
404 outline_base_bin(lgTransStackLine, ip, phots, inwd, nonScatteredFraction);
405 else
406 {
407 ASSERT( damp > 0. );
408 double LineWidth = dampXvel/damp;
409 LineWidth = MIN2( 0.1 * SPEEDLIGHT, LineWidth );
410 double sigma = (LineWidth/SPEEDLIGHT);
411 long ip3SigmaRed = ipoint( MAX2( rfield.emm, rfield.anu[ip] - 3.*sigma*rfield.anu[ip] ) );
412 long ip3SigmaBlue = ipoint( MIN2( rfield.egamry, rfield.anu[ip] + 3.*sigma*rfield.anu[ip] ) );
413 ASSERT( ip3SigmaBlue >= ip3SigmaRed );
414 long numBins = ip3SigmaBlue - ip3SigmaRed + 1;
415
416 if( numBins < 3 )
417 outline_base_bin(lgTransStackLine, ip, phots, inwd, nonScatteredFraction);
418 else
419 {
420 valarray<realnum> x(numBins);
421 valarray<realnum> profile(numBins);
422
423 for( long ipBin=ip3SigmaRed; ipBin<=ip3SigmaBlue; ipBin++ )
424 x[ipBin] = (rfield.anu[ip] - rfield.anu[ipBin])/rfield.anu[ip]/sigma;
425 // this profile must have unit normalization to conserve energy -> use U(a,v)
426 VoigtU(damp,get_ptr(x),get_ptr(profile),numBins);
427
428 for( long ipBin=ip3SigmaRed; ipBin<=ip3SigmaBlue; ipBin++ )
429 outline_base_bin(lgTransStackLine, ipBin, phots*profile[ipBin-ip3SigmaRed], inwd, nonScatteredFraction);
430 }
431 }
432}
433
434/*outline_base_bin - adds line photons to bins of reflin and outlin */
435void outline_base_bin(bool lgTransStackLine, long int ip, double phots, realnum inwd,
436 double nonScatteredFraction)
437{
438 DEBUG_ENTRY( "outline_base_bin()" );
439
440 if (lgTransStackLine)
441 {
442 rfield.DiffuseLineEmission[ip] +=
443 (realnum)phots;
444
445 /* the reflected part */
446 rfield.reflin[0][ip] +=
447 (realnum)(inwd*phots*radius.BeamInIn);
448
449 /* inward beam that goes out since sphere set */
450 rfield.outlin[0][ip] +=
451 (realnum)(inwd*phots*radius.BeamInOut*opac.tmn[ip]*nonScatteredFraction);
452
453 /* outward part */
454 rfield.outlin[0][ip] +=
455 (realnum)((1.-inwd)*phots*radius.BeamOutOut*opac.tmn[ip]*nonScatteredFraction);
456 }
457 else
458 {
459 rfield.reflin[0][ip] +=
460 (realnum)(phots*radius.dVolReflec);
461
462 rfield.outlin[0][ip] +=
463 (realnum)(phots*radius.dVolOutwrd*opac.ExpZone[ip]);
464 }
465}
466
467/*lindst add line with destruction and outward */
469 // xInten - local emissivity per unit vol
470 double xInten,
471 // wavelength of line in Angstroms
473 // *chLab string label for ion
474 const char *chLab,
475 // ipnt offset of line in continuum mesh
476 long int ipnt,
477 // chInfo character type of entry for line - 'c' cooling, 'h' heating, 'i' info only, 'r' recom line
478 char chInfo,
479 // lgOutToo should line be included in outward beam?
480 bool lgOutToo,
481 // *chComment string explaining line
482 const char *chComment )
483{
484 DEBUG_ENTRY( "lindst()" );
485
486 // do not add information lines to outward beam
487 ASSERT( !lgOutToo || chInfo!='i' );
488
489 lincom(xInten, wavelength, chLab, ipnt, chInfo, chComment, false );
490
491 if( LineSave.ipass > 0 )
492 {
493 /* >>chng 06 feb 08, add test on xInten positive, no need to evaluate
494 * for majority of zero */
495 if (lgOutToo && xInten > 0.)
496 {
497 /* add line to outward beam
498 * there are lots of lines that are sums of other lines, or
499 * just for info of some sort. These have flag lgOutToo false.
500 * Note that the EnergyRyd variable only has a rational
501 * value if PntForLine was called just before this routine - in all
502 * cases where this did not happen the flag is false. */
503 const bool lgTransStackLine = false;
504 const long int ip = ipnt - 1;
505 const double phots = xInten/(rfield.anu[ipnt-1]*EN1RYD);
506 const realnum inwd = (realnum)(1.0-(1.+geometry.covrt)/2.);
507 const double nonScatteredFraction = 1.;
508
509 outline_base_bin(lgTransStackLine, ip, phots, inwd, nonScatteredFraction);
510 }
511 }
512}
513
514/*lindst add line with destruction and outward */
516 double dampXvel,
517 double damp,
518 // xInten - local emissivity per unit vol
519 double xInten,
520 // wavelength of line in Angstroms
522 // *chLab string label for ion
523 const char *chLab,
524 // ipnt offset of line in continuum mesh
525 long int ipnt,
526 // chInfo character type of entry for line - 'c' cooling, 'h' heating, 'i' info only, 'r' recom line
527 char chInfo,
528 // lgOutToo should line be included in outward beam?
529 bool lgOutToo,
530 // *chComment string explaining line
531 const char *chComment )
532{
533 DEBUG_ENTRY( "lindst()" );
534
535 // do not add information lines to outward beam
536 ASSERT( !lgOutToo || chInfo!='i' );
537
538 lincom(xInten, wavelength, chLab, ipnt, chInfo, chComment, false );
539
540 if( LineSave.ipass > 0 )
541 {
542 /* >>chng 06 feb 08, add test on xInten positive, no need to evaluate
543 * for majority of zero */
544 if (lgOutToo && xInten > 0.)
545 {
546 /* add line to outward beam
547 * there are lots of lines that are sums of other lines, or
548 * just for info of some sort. These have flag lgOutToo false.
549 * Note that the EnergyRyd variable only has a rational
550 * value if PntForLine was called just before this routine - in all
551 * cases where this did not happen the flag is false. */
552 const bool lgTransStackLine = false;
553 const long int ip = ipnt - 1;
554 const double phots = xInten/(rfield.anu[ipnt-1]*EN1RYD);
555 const realnum inwd = (realnum)(1.0-(1.+geometry.covrt)/2.);
556 const double nonScatteredFraction = 1.;
557
558 outline_base(dampXvel, damp, lgTransStackLine, ip, phots, inwd, nonScatteredFraction);
559 }
560 }
561}
562
563/*lindst add line with destruction and outward */
565 const TransitionProxy& t,
566 // *chLab string label for ion
567 const char *chLab,
568 // chInfo character type of entry for line - 'c' cooling, 'h' heating, 'i' info only, 'r' recom line
569 char chInfo,
570 // lgOutToo should line be included in outward beam?
571 bool lgOutToo,
572 // *chComment string explaining line
573 const char *chComment )
574{
575 DEBUG_ENTRY( "lindst()" );
576
577 lindst( t.Emis().dampXvel(), t.Emis().damp(), t.Emis().xIntensity(), t.WLAng(), chLab, t.ipCont(), chInfo,
578 lgOutToo, chComment );
579
580}
581
582/*PntForLine generate pointer for forbidden line */
584 /* wavelength of transition in Angstroms */
585 double wavelength,
586 /* label for this line */
587 const char *chLabel,
588 /* this is array index on the f, not c scale,
589 * for the continuum cell holding the line */
590 long int *ipnt)
591{
592 /*
593 * maximum number of forbidden lines - this is a good bet since
594 * new lines do not go into this group, and lines are slowly
595 * moving to level 1
596 */
597 const int MAXFORLIN = 1000;
598 static long int ipForLin[MAXFORLIN]={0};
599
600 /* number of forbidden lines entered into continuum array */
601 static long int nForLin;
602
603 DEBUG_ENTRY( "PntForLine()" );
604
605 /* must be 0 or greater */
606 ASSERT( wavelength >= 0. );
607
608 if( wavelength == 0. )
609 {
610 /* zero is special flag to initialize */
611 nForLin = 0;
612 }
613 else
614 {
615
616 if( LineSave.ipass > 0 )
617 {
618 /* not first pass, sum lines only */
619 *ipnt = ipForLin[nForLin];
620 }
621 else if( LineSave.ipass == 0 )
622 {
623 /* check if number of lines in arrays exceeded */
624 if( nForLin >= MAXFORLIN )
625 {
626 fprintf( ioQQQ, "PROBLEM %5ld lines is too many for PntForLine.\n",
627 nForLin );
628 fprintf( ioQQQ, " Increase the value of maxForLine everywhere in the code.\n" );
630 }
631
632 /* ipLineEnergy will only put in line label if nothing already there */
633 const double EnergyRyd = RYDLAM/wavelength;
634 ipForLin[nForLin] = ipLineEnergy(EnergyRyd,chLabel , 0);
635 *ipnt = ipForLin[nForLin];
636 }
637 else
638 {
639 /* this is case where we are only counting lines */
640 *ipnt = 0;
641 }
642 ++nForLin;
643 }
644 return;
645}
646
647/*ConvRate2CS convert down coll rate back into electron cs in case other parts of code need this for reference */
648double ConvRate2CS( realnum gHi , realnum rate )
649{
650
651 double cs;
652
653 DEBUG_ENTRY( "ConvRate2CS()" );
654
655 /* return is collision strength, convert from collision rate from
656 * upper to lower, this assumes pure electron collisions, but that will
657 * also be assumed by anything that uses cs, for self-consistency */
658 cs = rate * gHi / dense.cdsqte;
659
660 /* change assert to non-negative - there can be cases (Iin H2) where cs has
661 * underflowed to 0 on some platforms */
662 ASSERT( cs >= 0. );
663 return cs;
664}
665
666/*ConvCrossSect2CollStr convert collisional deexcitation cross section for into collision strength */
667double ConvCrossSect2CollStr( double CrsSectCM2, double gLo, double E_ProjectileRyd, double reduced_mass_grams )
668{
669 double CollisionStrength;
670
671 DEBUG_ENTRY( "ConvCrossSect2CollStr()" );
672
673 ASSERT( CrsSectCM2 >= 0. );
674 ASSERT( gLo >= 0. );
675 ASSERT( E_ProjectileRyd >= 0. );
676 ASSERT( reduced_mass_grams >= 0. );
677
678 CollisionStrength = CrsSectCM2 * gLo * E_ProjectileRyd / (PI*BOHR_RADIUS_CM*BOHR_RADIUS_CM);
679
680 // this part is being tested.
681#if 0
682 CollisionStrength *= reduced_mass_grams / ELECTRON_MASS;
683#endif
684
685 ASSERT( CollisionStrength >= 0. );
686 return CollisionStrength;
687}
688
689/*totlin sum total intensity of cooling, recombination, or intensity lines */
690double totlin(
691 /* chInfor is 1 char,
692 'i' information,
693 'r' recombination or
694 'c' collision */
695 int chInfo)
696{
697 long int i;
698 double totlin_v;
699
700 DEBUG_ENTRY( "totlin()" );
701
702 /* routine goes through set of entered line
703 * intensities and picks out those which have
704 * types agreeing with chInfo. Valid types are
705 * 'c', 'r', and 'i'
706 *begin sanity check */
707 if( (chInfo != 'i' && chInfo != 'r') && chInfo != 'c' )
708 {
709 fprintf( ioQQQ, " TOTLIN does not understand chInfo=%c\n",
710 chInfo );
712 }
713 /*end sanity check */
714
715 /* now find sum of lines of type chInfo */
716 totlin_v = 0.;
717 for( i=0; i < LineSave.nsum; i++ )
718 {
719 if( LineSv[i].chSumTyp == chInfo )
720 {
721 totlin_v += LineSv[i].SumLine[0];
722 }
723 }
724 return( totlin_v );
725}
726
727
728/*FndLineHt search through line heat arrays to find the strongest heat source */
729const TransitionProxy FndLineHt(long int *level)
730{
731 long int i;
733 DEBUG_ENTRY( "FndLineHt()" );
734
735 double Strong = -1.;
736 *level = 0;
737
738 /* do the level 1 lines, 0 is dummy line, <=nLevel1 is correct for c scale */
739 for( i=1; i <= nLevel1; i++ )
740 {
741 /* check if a line was the major heat agent */
742 if( TauLines[i].Coll().heat() > Strong )
743 {
744 *level = 1;
745 t = TauLines[i];
746 Strong = TauLines[i].Coll().heat();
747 }
748 }
749
750 /* now do the level 2 lines */
751 for( i=0; i < nWindLine; i++ )
752 {
753 if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
754 {
755 /* check if a line was the major heat agent */
756 if( TauLine2[i].Coll().heat() > Strong )
757 {
758 *level = 2;
759 t = TauLine2[i];
760 Strong = TauLine2[i].Coll().heat();
761 }
762 }
763 }
764
765 /* now do the hyperfine structure lines */
766 for( i=0; i < nHFLines; i++ )
767 {
768 /* check if a line was the major heat agent */
769 if( HFLines[i].Coll().heat() > Strong )
770 {
771 *level = 3;
772 t = HFLines[i];
773 Strong = HFLines[i].Coll().heat();
774 }
775 }
776
777 /* lines from external databases */
778 for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
779 {
780 for( EmissionList::iterator em=dBaseTrans[ipSpecies].Emis().begin();
781 em != dBaseTrans[ipSpecies].Emis().end(); ++em)
782 {
783 /* check if a line was the major heat agent */
784 if( (*em).Tran().Coll().heat() > Strong )
785 {
786 *level = 4;
787 t = (*em).Tran();
788 Strong = t.Coll().heat();
789 }
790 }
791 }
792
793 fixit(); // all other line stacks need to be included here.
794 // can we just sweep over line stack? Is that ready yet?
795
796 ASSERT( t.associated() );
797 return t;
798}
799
FILE * ioQQQ
Definition cddefines.cpp:7
long int iteration
Definition cddefines.cpp:16
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define MIN2
Definition cddefines.h:761
#define MALLOC(exp)
Definition cddefines.h:501
#define POW2
Definition cddefines.h:929
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int NISO
Definition cddefines.h:261
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
T * get_ptr(T *v)
Definition cddefines.h:1079
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void fixit(void)
Definition service.cpp:991
LinSv * LineSvSortWL
Definition cdinit.cpp:71
long nWindLine
Definition cdinit.cpp:19
LinSv * LineSv
Definition cdinit.cpp:70
double & heat() const
Definition collision.h:194
EmissionProxy::iterator iterator
Definition emission.h:317
realnum & damp() const
Definition emission.h:563
double & xIntensity() const
Definition emission.h:483
realnum & dampXvel() const
Definition emission.h:553
CollisionProxy Coll() const
Definition transition.h:424
bool associated() const
Definition transition.h:50
realnum & WLAng() const
Definition transition.h:429
long & ipCont() const
Definition transition.h:450
EmissionList::reference Emis() const
Definition transition.h:408
long ipoint(double energy_ryd)
long ipLineEnergy(double energy, const char *chLabel, long ipIonEnergy)
t_dense dense
Definition dense.cpp:24
t_geometry geometry
Definition geometry.cpp:5
t_LineSave LineSave
Definition lines.cpp:5
void lines(void)
Definition prt_lines.cpp:34
struct t_tag_LineSv LinSv
const TransitionProxy FndLineHt(long int *level)
void linadd(double xInten, realnum wavelength, const char *chLab, char chInfo, const char *chComment)
void outline_base(double dampXvel, double damp, bool lgTransStackLine, long int ip, double phots, realnum inwd, double nonScatteredFraction)
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)
double RefIndex(double EnergyWN)
double eina(double gf, double enercm, double gup)
double ConvRate2CS(realnum gHi, realnum rate)
double abscf(double gf, double enercm, double gl)
void outline_base_bin(bool lgTransStackLine, long int ip, double phots, realnum inwd, double nonScatteredFraction)
double ConvCrossSect2CollStr(double CrsSectCM2, double gLo, double E_ProjectileRyd, double reduced_mass_grams)
double GetGF(double trans_prob, double enercm, double gup)
double totlin(int chInfo)
void LineStackCreate()
realnum WavlenErrorGet(realnum wavelength)
double emergent_line(double emissivity_in, double emissivity_out, long int ipCont)
STATIC void lincom(double xInten, realnum wavelength, const char *chLab, long int ipnt, char chInfo, const char *chComment, bool lgAdd)
static realnum * wavelength
t_opac opac
Definition opacity.cpp:5
UNUSED const double TRANS_PROB_CONST
Definition physconst.h:237
UNUSED const double PI
Definition physconst.h:29
UNUSED const double BOHR_RADIUS_CM
Definition physconst.h:222
UNUSED const double SPEEDLIGHT
Definition physconst.h:100
UNUSED const double ELECTRON_MASS
Definition physconst.h:91
UNUSED const double EN1RYD
Definition physconst.h:179
UNUSED const double RYDLAM
Definition physconst.h:176
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
t_rt rt
Definition rt.cpp:5
long int nSpecies
Definition taulines.cpp:21
TransitionList TauLine2("TauLine2", &AnonStates)
vector< TransitionList > dBaseTrans
Definition taulines.cpp:17
TransitionList HFLines("HFLines", &AnonStates)
long int nLevel1
Definition taulines.cpp:28
long int nHFLines
Definition taulines.cpp:31
TransitionList TauLines("TauLines", &AnonStates)
void VoigtU(realnum a, const realnum v[], realnum y[], int n)
Definition thirdparty.h:364
t_trace trace
Definition trace.cpp:5