cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
rt_line_driving.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*RT_line_driving derive radiative acceleration due to line absorption of incident continuum,
4  * return value is line radiative acceleration */
5 #include "cddefines.h"
6 #include "rt.h"
7 #include "iso.h"
8 #include "dense.h"
9 #include "taulines.h"
10 #include "h2.h"
11 #include "atmdat.h"
12 
13 /*RT_line_driving derive radiative acceleration due to line absorption of incident continuum,
14  * return value is line radiative acceleration */
15 double RT_line_driving(void)
16 {
17  long int ipHi,
18  nelem,
19  ipLo,
20  ipISO;
21 
22  double AllRest,
23  OneLine,
24  h2drive,
25  accel_iso[NISO];
26 
27  /* following used for debugging */
28  /* double
29  RestMax,
30  HeavMax,
31  hydromax;
32  long int
33  ipRestMax,
34  ihmax; */
35 
36  DEBUG_ENTRY( "RT_line_driving()" );
37 
38  /* this function finds the total rate the gas absorbs energy
39  * this result is divided by the calling routine to find the
40  * momentum absorbed by the gas, and eventually the radiative acceleration
41  *
42  * the energy absorbed by the line is
43  * Abundance * energy * A *(g_up/g_lo) * occnum * escape prob
44  * where occnum is the photon occupation number, and the g's are
45  * the ratios of statistical weights */
46 
47  /* total energy absorbed in this zone, per cubic cm
48  * do hydrogen first */
49 
50  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
51  {
52  accel_iso[ipISO] = 0;
53  for( nelem=ipISO; nelem < LIMELM; nelem++ )
54  {
55  if( (dense.IonHigh[nelem] >= nelem + 1-ipISO) )
56  {
57  for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
58  {
59  /* do not put in highest level since its not real */
60  for( ipLo=0; ipLo < ipHi - 1; ipLo++ )
61  {
62  /* do not include bogus lines */
63  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() > 0 )
64  {
65  OneLine = iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().pump()*
66  iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyErg()*
67  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().PopOpc();
68 
69  accel_iso[ipISO] += OneLine;
70  }
71  }
72  }
73 
74  if( iso_ctrl.lgDielRecom[ipISO] )
75  {
76  // SatelliteLines are indexed by lower level, summed over satellite levels
77  for( ipLo=0; ipLo < iso_sp[ipISO][nelem].numLevels_local; ipLo++ )
78  {
79  /* do not include bogus lines */
80  TransitionList::iterator tr = SatelliteLines[ipISO][nelem].begin()+ipSatelliteLines[ipISO][nelem][ipLo];
81  if((*tr).ipCont() > 0 )
82  {
83  OneLine = (*tr).Emis().pump()*
84  (*tr).EnergyErg()*
85  (*tr).Emis().PopOpc();
86 
87  accel_iso[ipISO] += OneLine;
88  }
89  }
90  }
91 
92  for( ipHi=iso_sp[ipISO][nelem].st[iso_sp[ipISO][nelem].numLevels_local-1].n()+1; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )
93  {
94  /* do not include bogus lines */
95  TransitionList::iterator tr = ExtraLymanLines[ipISO][nelem].begin()+ipExtraLymanLines[ipISO][nelem][ipHi];
96  if( (*tr).ipCont() > 0 )
97  {
98  OneLine = (*tr).Emis().pump()*
99  (*tr).EnergyErg()*
100  (*tr).Emis().PopOpc();
101 
102  accel_iso[ipISO] += OneLine;
103  }
104 
105  }
106  }
107  }
108  }
109 
110  /* all heavy element lines treated with g-bar
111  * these are the level 2 lines, f should be ok */
112  AllRest = 0.;
113  for( long i=0; i < nWindLine; i++ )
114  {
115  OneLine =
116  TauLine2[i].Emis().pump()*
117  TauLine2[i].EnergyErg()*
118  TauLine2[i].Emis().PopOpc();
119  AllRest += OneLine;
120  }
121  for( size_t i=0; i < UTALines.size(); i++ )
122  {
123  OneLine =
124  UTALines[i].Emis().pump()*
125  UTALines[i].EnergyErg()*
126  UTALines[i].Emis().PopOpc();
127  AllRest += OneLine;
128  }
129  for( size_t i=0; i < HFLines.size(); i++ )
130  {
131  OneLine =
132  HFLines[i].Emis().pump()*
133  HFLines[i].EnergyErg()*
134  HFLines[i].Emis().PopOpc();
135  AllRest += OneLine;
136  }
137 
138  for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
139  {
140  if( dBaseSpecies[ipSpecies].lgActive )
141  {
142  for (TransitionList::iterator tr=dBaseTrans[ipSpecies].begin();
143  tr != dBaseTrans[ipSpecies].end(); ++tr)
144  {
145  int ipHi = (*tr).ipHi();
146  if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local || (*tr).ipCont() <= 0)
147  continue;
148  OneLine = (*tr).EnergyErg()*(*tr).Emis().pump()*(*tr).Emis().PopOpc();
149  AllRest += OneLine;
150  }
151  }
152  }
153 
154  /* the H2 molecule */
155  h2drive = 0.;
156  for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
157  h2drive += (*diatom)->H2_Accel();
158 
159  // total radiative acceleration
160  double forlin_v = accel_iso[ipH_LIKE] + accel_iso[ipHE_LIKE] + h2drive + AllRest;
161  if( 0 )
162  {
163  fprintf(ioQQQ," wind te %e %e %e %e %e\n",
164  accel_iso[ipH_LIKE], accel_iso[ipHE_LIKE], h2drive, AllRest, forlin_v );
165  }
166 
167 
168  {
169  enum {DEBUG_LOC=false};
170  if( DEBUG_LOC)
171  {
172  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
173  {
174  for( nelem=ipISO; nelem < LIMELM; nelem++ )
175  {
176  if( (dense.IonHigh[nelem] >= nelem + 1-ipISO) )
177  {
178  for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
179  {
180  /* do not put in highest level since its not real */
181  for( ipLo=0; ipLo < ipHi - 1; ipLo++ )
182  {
183  /* do not include bogus lines */
184  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() > 0 )
185  {
186  OneLine = iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().pump()*
187  iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyErg()*
188  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().PopOpc();
189  if( OneLine / forlin_v > 0.03 )
190  {
191  fprintf(ioQQQ,"DEBUG OneLine %li %li %.2f %.2e\n",
192  ipISO,nelem,iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyAng() ,
193  OneLine/forlin_v);
194  }
195 
196  }
197  }
198  }
199  }
200  }
201  }
202 
203  for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
204  {
205  if( dBaseSpecies[ipSpecies].lgActive )
206  {
207  for (TransitionList::iterator tr=dBaseTrans[ipSpecies].begin();
208  tr != dBaseTrans[ipSpecies].end(); ++tr)
209  {
210  int ipHi = (*tr).ipHi();
211  if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local || (*tr).ipCont() <= 0)
212  continue;
213  OneLine = (*tr).EnergyErg()*(*tr).Emis().pump()*(*tr).Emis().PopOpc();
214  if( OneLine / forlin_v > 0.03 ||
215  (ipSpecies==atmdat.ipSpecIon[ipCARBON][3] && (*tr).ipHi()==1 )
216  )
217  {
218  fprintf(ioQQQ,"DEBUG OneLine %s %.2f %.2e\n",
219  dBaseSpecies[ipSpecies].chLabel, (*tr).EnergyAng() , OneLine/forlin_v);
220  }
221  }
222  }
223  }
224  }
225  }
226 
227  return( forlin_v );
228 }
t_atmdat atmdat
Definition: atmdat.cpp:6
size_t size(void) const
Definition: transition.h:331
realnum EnergyErg() const
Definition: transition.h:90
TransitionList UTALines("UTALines",&AnonStates)
const int ipHE_LIKE
Definition: iso.h:65
multi_arr< int, 3 > ipSatelliteLines
Definition: taulines.cpp:34
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
const int NISO
Definition: cddefines.h:311
long int IonHigh[LIMELM+1]
Definition: dense.h:130
TransitionList HFLines("HFLines",&AnonStates)
FILE * ioQQQ
Definition: cddefines.cpp:7
TransitionList TauLine2("TauLine2",&AnonStates)
long int nSpecies
Definition: taulines.cpp:22
t_dense dense
Definition: global.cpp:15
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
long int nLyman[NISO]
Definition: iso.h:352
EmissionList::reference Emis() const
Definition: transition.h:447
multi_arr< int, 3 > ipExtraLymanLines
Definition: taulines.cpp:24
double RT_line_driving(void)
bool lgDielRecom[NISO]
Definition: iso.h:385
vector< diatomics * > diatoms
Definition: h2.cpp:8
vector< vector< TransitionList > > SatelliteLines
Definition: taulines.cpp:35
long nWindLine
Definition: cdinit.cpp:19
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:473
const int ipH_LIKE
Definition: iso.h:64
vector< vector< TransitionList > > ExtraLymanLines
Definition: taulines.cpp:25
const int LIMELM
Definition: cddefines.h:308
vector< vector< long > > ipSpecIon
Definition: atmdat.h:455
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
double & PopOpc() const
Definition: emission.h:670
vector< species > dBaseSpecies
Definition: taulines.cpp:15
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
const int ipCARBON
Definition: cddefines.h:354
vector< TransitionList > dBaseTrans
Definition: taulines.cpp:18
long int numLevels_local
Definition: iso.h:529
EmissionList & Emis()
Definition: transition.h:363
vector< diatomics * >::iterator diatom_iter
Definition: h2.h:13
double & pump() const
Definition: emission.h:530