cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
opacity_add1subshell.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2017 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*OpacityAdd1Subshell add opacity due to single shell to main opacity array*/
4 /*OpacityAdd1SubshellInduc add opacity of individual species, including stimulated emission */
5 #include "cddefines.h"
6 #include "rfield.h"
7 #include "hydrogenic.h"
8 #include "opacity.h"
9 
11  /*ipOpac is opacity index within opac opacity offset for this species */
12  long int ipOpac,
13  /* lower freq limit to opacity range on energy mesh */
14  long int ipLowLim,
15  /* upper limit to opacity range on energy mesh */
16  long int ipUpLim,
17  /* abundance, we bail if zero */
18  realnum abundance,
19  /* either static 's' or volitile 'v' */
20  char chStat )
21 {
22  long int i ,
23  ipOffset,
24  limit;
25 
26  DEBUG_ENTRY( "OpacityAdd1Subshell()" );
27 
28  /* code spends roughly 20% of its time in this loop*/
29 
30  ASSERT( chStat == 's' || chStat == 'v' );
31 
32  ipOffset = ipOpac - ipLowLim;
33  ASSERT( ipLowLim > 0 );
34  /* >>chng 02 aug 13, negative offset is ok, it is only this plus ipLowLim that matters */
35  /*ASSERT( ipOffset >= 0 );*/
36 
37  limit = MIN2(ipUpLim,rfield.nflux);
38 
39  /* do nothing if abundance is zero, or if static opacities do not
40  * need to be redone */
41  if( abundance <= 0. || (chStat=='s' && !opac.lgRedoStatic) )
42  {
43  return;
44  }
45 
46  /* volative (outer shell, constantly reevaluated) or static opacity? */
47  if( chStat=='v' )
48  {
49  for( i=ipLowLim-1; i < limit; i++ )
50  {
51  opac.opacity_abs[i] += opac.OpacStack[i+ipOffset]*abundance;
52  }
53  }
54  else
55  {
56  for( i=ipLowLim-1; i < limit; i++ )
57  {
58  opac.OpacStatic[i] += opac.OpacStack[i+ipOffset]*abundance;
59  }
60  }
61  return;
62 }
63 
64 /*OpacityAdd1SubshellInduc add opacity of individual species, including stimulated emission */
66  /* pointer to opacity within opacity stack */
67  long int ipOpac,
68  /* pointer to low energy in continuum array for this opacity band */
69  long int ipLowEnergy,
70  /* pointer to high energy in continuum array for this opacity */
71  long int ipHiEnergy,
72  /* this abundance of this species, may be zero */
73  double abundance,
74  /* the departure coef, may be infinite or zero */
75  double DepartCoef ,
76  /* either 'v' for volitile or 's' for static opacities */
77  char chStat )
78 {
79  long int i,
80  iup,
81  k;
82 
83  DEBUG_ENTRY( "OpacityAdd1SubshellInduc()" );
84 
85  /* add opacity of individual species, including stimulated emission
86  * abundance is the density of the lower level (cm^-3)
87  * DepartCoef is its departure coefficient, can be zero */
88 
89  /* this is opacity offset, must be positive */
90  ASSERT( ipOpac > 0 );
91 
92  /* check that chStat is either 'v' or 's' */
93  ASSERT( chStat == 'v' || chStat == 's' );
94 
95  /* do nothing if abundance is zero, or if static opacities do not
96  * need to be redone */
97  if( abundance <= 0. || (chStat=='s' && !opac.lgRedoStatic) )
98  {
99  return;
100  }
101 
102  k = ipOpac - ipLowEnergy;
103 
104  /* DepartCoef is dep coef, rfield.lgInducProcess is turned off with 'no indcued' command */
105  if( (DepartCoef > 1e-35 && rfield.lgInducProcess) && hydro.lgHInducImp )
106  {
107  iup = MIN2(ipHiEnergy,rfield.nflux);
108  /* >>>chng 99 apr 29, following was present, caused pdr to make opac at impossible energy*/
109  /*iup = MAX2(ipLowEnergy,iup);*/
110  double DepartCoefInv = 1./DepartCoef;
111  if( chStat == 'v' )
112  {
113  /* volitile opacities, always reevaluate */
114  for( i=ipLowEnergy-1; i < iup; i++ )
115  {
116  opac.opacity_abs[i] += opac.OpacStack[i+k]*abundance*
117  MAX2(0. , 1.- rfield.ContBoltz[i]*DepartCoefInv);
118  }
119  }
120  else
121  {
122  /* static opacities, save in special array */
123  for( i=ipLowEnergy-1; i < iup; i++ )
124  {
125  opac.OpacStatic[i] += opac.OpacStack[i+k]*abundance*
126  MAX2(0. , 1.- rfield.ContBoltz[i]*DepartCoefInv);
127  }
128  }
129  }
130 
131  else
132  {
133  /* DepartCoef is the departure coef, which can go to zero,
134  * neglect stimulated emission in this case */
135  iup = MIN2(ipHiEnergy,rfield.nflux);
136  /* >>>chng 99 apr 29, following was present, caused pdr to make opac at impossible energy*/
137  /*iup = MAX2(ipLowEnergy,iup);*/
138  if( chStat == 'v' )
139  {
140  for( i=ipLowEnergy-1; i < iup; i++ )
141  {
142  opac.opacity_abs[i] += opac.OpacStack[i+k]*abundance;
143  }
144  }
145  else
146  {
147  for( i=ipLowEnergy-1; i < iup; i++ )
148  {
149  opac.OpacStatic[i] += opac.OpacStack[i+k]*abundance;
150  }
151  }
152  }
153 
154  return;
155 }
double * OpacStack
Definition: opacity.h:163
double * opacity_abs
Definition: opacity.h:103
t_opac opac
Definition: opacity.cpp:5
double * OpacStatic
Definition: opacity.h:122
bool lgHInducImp
Definition: hydrogenic.h:131
void OpacityAdd1Subshell(long int ipOpac, long int ipLowLim, long int ipUpLim, realnum abundance, char chStat)
#define MIN2(a, b)
Definition: cddefines.h:807
t_rfield rfield
Definition: rfield.cpp:9
float realnum
Definition: cddefines.h:124
t_hydro hydro
Definition: hydrogenic.cpp:5
double * ContBoltz
Definition: rfield.h:126
#define ASSERT(exp)
Definition: cddefines.h:617
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
bool lgRedoStatic
Definition: opacity.h:159
#define MAX2(a, b)
Definition: cddefines.h:828
bool lgInducProcess
Definition: rfield.h:235
void OpacityAdd1SubshellInduc(long int ipOpac, long int low, long int ihi, double a, double b, char chStat)
long int nflux
Definition: rfield.h:48