cloudy trunk
Loading...
Searching...
No Matches
prt_lines_general.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/*lines_general put general information and energetics into line intensity stack */
4#include "cddefines.h"
5#include "taulines.h"
6#include "coolheavy.h"
7#include "hydrogenic.h"
8#include "dense.h"
9#include "thermal.h"
10#include "continuum.h"
11#include "geometry.h"
12#include "dynamics.h"
13#include "rt.h"
14#include "iso.h"
15#include "rfield.h"
16#include "trace.h"
17#include "ionbal.h"
18#include "lines_service.h"
19#include "radius.h"
20#include "lines.h"
21
22void lines_general(void)
23{
24 long int i,
25 ipHi,
26 ipLo,
27 nelem,
28 ipnt;
29
30 double
31 hbetac,
32 HeatMetal ,
33 ee511,
34 hlalph;
35
36 DEBUG_ENTRY( "lines_general()" );
37
38 if( trace.lgTrace )
39 {
40 fprintf( ioQQQ, " lines_general called\n" );
41 }
42
43 i = StuffComment( "general properties" );
44 linadd( 0., (realnum)i , "####", 'i',
45 " start of general properties");
46
47 /* total H-beta from multi-level atom */
48 nelem = ipHYDROGEN;
49 ipLo = ipH2p;
50
51 // this can be changed with the atom levels command but must be at
52 // least 3
53 ASSERT( iso_sp[ipH_LIKE][nelem].n_HighestResolved_max >=3 );
54
55 if( iso_sp[ipH_LIKE][nelem].n_HighestResolved_max >= 4 )
56 {
57 ipHi = ipH4p;
58 hbetac =
59 (iso_sp[ipH_LIKE][nelem].trans(ipH4p,ipH2s).Emis().Aul() *
60 (iso_sp[ipH_LIKE][nelem].trans(ipH4p,ipH2s).Emis().Pesc() +
61 iso_sp[ipH_LIKE][nelem].trans(ipH4p,ipH2s).Emis().Pelec_esc() ) *
62 iso_sp[ipH_LIKE][nelem].st[ipH4p].Pop() +
63 iso_sp[ipH_LIKE][nelem].trans(ipH4s,ipH2p).Emis().Aul() *
64 (iso_sp[ipH_LIKE][nelem].trans(ipH4s,ipH2p).Emis().Pesc() +
65 iso_sp[ipH_LIKE][nelem].trans(ipH4s,ipH2p).Emis().Pelec_esc() ) *
66 iso_sp[ipH_LIKE][nelem].st[ipH4s].Pop() +
67 iso_sp[ipH_LIKE][nelem].trans(ipH4d,ipH2p).Emis().Aul() *
68 (iso_sp[ipH_LIKE][nelem].trans(ipH4d,ipH2p).Emis().Pesc() +
69 iso_sp[ipH_LIKE][nelem].trans(ipH4d,ipH2p).Emis().Pelec_esc() ) *
70 iso_sp[ipH_LIKE][nelem].st[ipH4d].Pop()) *
71 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).EnergyErg();
72 }
73 else
74 {
75 // atom levels command does not allow < 3
76 ASSERT( iso_sp[ipH_LIKE][nelem].n_HighestResolved_max == 3 );
77 ipHi = 6;
78 hbetac =
79 (iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH2s).Emis().Aul() *
80 (iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH2s).Emis().Pesc() +
81 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH2s).Emis().Pelec_esc() ) +
82 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH2p).Emis().Aul() *
83 (iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH2p).Emis().Pesc() +
84 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH2p).Emis().Pelec_esc() ) ) *
85 iso_sp[ipH_LIKE][nelem].st[ipHi].Pop() *
86 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).EnergyErg();
87 }
88
89 /* these lines added to outlin in metdif - following must be false
90 * this passes array index for line energy in continuum mesh - in rest
91 * of code this is set by a previous call to PntForLine, this index
92 * is on the f not c scale */
93 rt.fracin = iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH2s).Emis().FracInwd();
94 lindst(hbetac,iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH2s).WLAng(),"TOTL",
95 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH2s).ipCont(),'t',false,
96 " H I Balmer beta predicted by model atom " );
97 rt.fracin = 0.5;
98
99 if( iso_sp[ipH_LIKE][nelem].n_HighestResolved_max < 4 )
100 {
101 // we need to have something for Hb "H 1" and "Inwd"
102 lindst(hbetac,iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH2s).WLAng(),"H 1",
103 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH2s).ipCont(),'t',false,
104 " H I Balmer beta predicted by model atom " );
105 // we need to have something for Hb "H 1" and "Inwd"
106 lindst(hbetac/2.,iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH2s).WLAng(),"Inwd",
107 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH2s).ipCont(),'t',false,
108 " H I Balmer beta predicted by model atom " );
109 }
110
111 /* total Ly-a from multi-level atom */
112 ipHi = ipH2p;
113 ipLo = ipH1s;
114 hlalph =
115 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul()*
116 iso_sp[ipH_LIKE][nelem].st[ipHi].Pop()*
117 (iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pesc() +
118 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pelec_esc() )*
119 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).EnergyErg();
120
121 rt.fracin = iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().FracInwd();
122 lindst(hlalph,iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).WLAng(),"TOTL",
123 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).ipCont(),'t',false ,
124 " H I Lya predicted from model atom ");
125 rt.fracin = 0.5;
126
127 /* this entry only works correctly if the APERTURE command is not in effect */
128 if( geometry.iEmissPower == 2 )
129 {
130 linadd(continuum.totlsv/radius.dVeffAper,0,"Inci",'i',
131 "total luminosity in incident continuum");
132 /* ipass is flag to indicate whether to only set up line array
133 * (ipass=0) or actually evaluate lines intensities (ipass=1) */
134 if( LineSave.ipass > 0 )
135 {
136 continuum.totlsv = 0.;
137 }
138 }
139
140 linadd(thermal.htot,0,"TotH",'i',
141 " total heating, all forms, information since individuals added later ");
142
143 linadd(thermal.ctot,0,"TotC",'i',
144 " total cooling, all forms, information since individuals added later ");
145
146 linadd(thermal.heating[0][0],0,"BFH1",'h',
147 " hydrogen photoionization heating, ground state only ");
148
149 linadd(thermal.heating[0][1],0,"BFHx",'h',
150 " net hydrogen photoionization heating less rec cooling, all excited states normally zero, positive if excited states are net heating ");
151
152 linadd(thermal.heating[0][22],0,"Line",'h',
153 " heating due to induced lines absorption of continuum ");
154 if( thermal.htot > 0. )
155 {
156 if( thermal.heating[0][22]/thermal.htot > thermal.HeatLineMax )
157 {
158 thermal.HeatLineMax = (realnum)(thermal.heating[0][22]/thermal.htot);
159 }
160 }
161
162 linadd(thermal.heating[1][0]+thermal.heating[1][1]+thermal.heating[1][2],0,"BFHe",'h',
163 " total helium photoionization heating, all stages ");
164
165 HeatMetal = 0.;
166 /* some sums that will be printed in the stack */
167 for( nelem=2; nelem<LIMELM; ++nelem)
168 {
169 /* we now have final solution for this element */
170 for( i=dense.IonLow[nelem]; i < dense.IonHigh[nelem]; i++ )
171 {
172 ASSERT( i < LIMELM );
173 /* total metal photo heating for LINES */
174 HeatMetal += thermal.heating[nelem][i];
175 }
176 }
177
178 linadd(HeatMetal,0,"TotM",'h',
179 " total heavy element photoionization heating, all stages ");
180
181 linadd(thermal.heating[0][21],0,"pair",'h',
182 " heating due to pair production ");
183
184 /* ipass is flag to indicate whether to only set up line array
185 * (ipass=0) or actually evaluate lines intensities (ipass=1) */
186 if( LineSave.ipass > 0 )
187 {
188 /* this will be max local heating due to bound compton */
189 ionbal.CompHeating_Max = MAX2( ionbal.CompHeating_Max , ionbal.CompRecoilHeatLocal/thermal.htot);
190 }
191 else
192 {
193 ionbal.CompHeating_Max = 0.;
194 }
195
196 linadd(ionbal.CompRecoilHeatLocal,0,"Cbnd",'h',
197 " heating due to bound compton scattering ");
198
199 linadd(rfield.cmheat,0,"ComH",'h',
200 " Compton heating ");
201
202 linadd(CoolHeavy.tccool,0,"ComC",'c',
203 " total Compton cooling ");
204
205 /* record max local heating due to advection */
206 dynamics.HeatMax = MAX2( dynamics.HeatMax , dynamics.Heat() /thermal.htot );
207 /* record max local cooling due to advection */
208 dynamics.CoolMax = MAX2( dynamics.CoolMax , dynamics.Cool() /thermal.htot );
209
210 linadd(dynamics.Cool() , 0 , "advC" , 'i',
211 " cooling due to advection " );
212
213 linadd(dynamics.Heat() , 0 , "advH" , 'i' ,
214 " heating due to advection ");
215
216 linadd( thermal.char_tran_heat ,0,"CT H",'h',
217 " heating due to charge transfer ");
218
219 linadd( thermal.char_tran_cool ,0,"CT C",'c',
220 " cooling due to charge transfer ");
221
222 linadd(thermal.heating[1][6],0,"CR H",'h',
223 " cosmic ray heating ");
224
225 linadd(thermal.heating[0][20],0,"extH",'h',
226 " extra heat added to this zone, from HEXTRA command ");
227
228 linadd(CoolHeavy.cextxx,0,"extC",'c',
229 " extra cooling added to this zone, from CEXTRA command ");
230
231 // 511 keV annihilation line, counts as recombination line since
232 // neither heating nor cooling, but does remove energy
233 ee511 = (dense.gas_phase[ipHYDROGEN] + 4.*dense.gas_phase[ipHELIUM])*ionbal.PairProducPhotoRate[0]*2.*8.20e-7;
234 PntForLine(2.427e-2,"e-e+",&ipnt);
235 lindst(ee511,(realnum)2.427e-2,"e-e+",ipnt,'r',true,
236 " 511keV annihilation line " );
237
238 linadd(CoolHeavy.expans,0,"Expn",'c',
239 " expansion cooling, only non-zero for wind ");
240
241 linadd(iso_sp[ipH_LIKE][ipHYDROGEN].RadRecCool,0,"H FB",'i',
242 " H radiative recombination cooling ");
243
244 linadd(MAX2(0.,iso_sp[ipH_LIKE][ipHYDROGEN].FreeBnd_net_Cool_Rate),0,"HFBc",'c',
245 " net free-bound cooling ");
246
247 linadd(MAX2(0.,-iso_sp[ipH_LIKE][ipHYDROGEN].FreeBnd_net_Cool_Rate),0,"HFBh",'h',
248 " net free-bound heating ");
249
250 linadd(iso_sp[ipH_LIKE][ipHYDROGEN].RecomInducCool_Rate,0,"Hind",'c',
251 " cooling due to induced rec of hydrogen ");
252
253 linadd(CoolHeavy.cyntrn,0,"Cycn",'c',
254 " cyclotron cooling ");
255
256 // cooling due to database species
257 for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
258 {
259 // label may be too long for linadd
260 char chLabel[5];
261 strncpy( chLabel, dBaseStates[ipSpecies][0].chLabel(), 4 );
262 chLabel[4] = '\0';
263 // this is information, 'i', since individual lines
264 // have been added as cooling or heating
265 linadd(dBaseSpecies[ipSpecies].CoolTotal,0, chLabel,'i',
266 " net cooling due to database species");
267 }
268
269 return;
270}
271
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
const int LIMELM
Definition cddefines.h:258
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
t_continuum continuum
Definition continuum.cpp:5
t_CoolHeavy CoolHeavy
Definition coolheavy.cpp:5
t_dense dense
Definition dense.cpp:24
t_dynamics dynamics
Definition dynamics.cpp:44
t_geometry geometry
Definition geometry.cpp:5
t_ionbal ionbal
Definition ionbal.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH1s
Definition iso.h:27
const int ipH4p
Definition iso.h:34
const int ipH4d
Definition iso.h:35
const int ipH2p
Definition iso.h:29
const int ipH4s
Definition iso.h:33
const int ipH2s
Definition iso.h:28
const int ipH_LIKE
Definition iso.h:62
t_LineSave LineSave
Definition lines.cpp:5
long int StuffComment(const char *chComment)
void linadd(double xInten, realnum wavelength, const char *chLab, char chInfo, const char *chComment)
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)
void lines_general(void)
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
vector< qList > dBaseStates
Definition taulines.cpp:15
species * dBaseSpecies
Definition taulines.cpp:14
t_thermal thermal
Definition thermal.cpp:5
t_trace trace
Definition trace.cpp:5