cloudy trunk
Loading...
Searching...
No Matches
atom_fe2ovr.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/*atoms_fe2ovr compute FeII overlap with Lya */
4/*fe2par evaluate FeII partition function */
5#include "cddefines.h"
6#include "doppvel.h"
7#include "dense.h"
8#include "rfield.h"
9#include "iso.h"
10#include "phycon.h"
11#include "hydrogenic.h"
12#include "ipoint.h"
13#include "opacity.h"
14#include "radius.h"
15#include "atomfeii.h"
16#include "thirdparty.h"
17#include "taulines.h"
18
19const double WLAL = 1215.6845;
20
21/* There are 373 transitions:
22 * Wavelength (A)
23 * absorption oscillator strength
24 * Energy of lower level (Ryd)
25 * Statistical weight of lower level (g) */
28{
29 DEBUG_ENTRY( "t_fe2ovr_la()" );
30
31 const long VERSION_MAGIC = 20070717L;
32 static const char chFile[] = "fe2ovr_la.dat";
33
34 FILE *io = open_data( chFile, "r" );
35
36 bool lgErr = false;
37 long i = -1L;
38
39 lgErr = lgErr || ( fscanf( io, "%ld", &i ) != 1 );
40 if( lgErr || i != VERSION_MAGIC )
41 {
42 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chFile, i );
43 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
45 }
46
47 double help=0;
48 for( i=0; i < NFEII; i++ )
49 {
50 lgErr = lgErr || ( fscanf( io, "%le", &help ) != 1 );
51 fe2lam[i] = (realnum)help;
52 }
53 for( i=0; i < NFEII; i++ )
54 {
55 lgErr = lgErr || ( fscanf( io, "%le", &help ) != 1 );
56 fe2osc[i] = (realnum)help;
57 }
58 for( i=0; i < NFEII; i++ )
59 {
60 lgErr = lgErr || ( fscanf( io, "%le", &help ) != 1 );
61 fe2enr[i] = (realnum)help;
62 }
63 for( i=0; i < NFEII; i++ )
64 {
65 lgErr = lgErr || ( fscanf( io, "%le", &help ) != 1 );
66 fe2gs[i] = (realnum)help;
67 }
68 for( i=0; i < NFE2PR; i++ )
69 lgErr = lgErr || ( fscanf( io, "%le", &fe2pt[i] ) != 1 );
70 for( i=0; i < NFE2PR; i++ )
71 lgErr = lgErr || ( fscanf( io, "%le", &fe2pf[i] ) != 1 );
72
73 fclose( io );
74
75 ASSERT( !lgErr );
76 return;
77}
78
80{
81 DEBUG_ENTRY( "zero_opacity()" );
82
83 for( int i=0; i < NFEII; i++ )
84 {
85 Fe2PopLte[i] = 0.f;
86 feopc[i] = 0.f;
87 Fe2TauLte[i] = opac.taumin;
88 }
89 return;
90}
91
93{
94 DEBUG_ENTRY( "init_pointers()" );
95
96 for( int i=0; i < NFEII; i++ )
97 ipfe2[i] = ipoint(fe2enr[i]);
98 return;
99}
100
103{
104 DEBUG_ENTRY( "tau_inc()" );
105
106 for( int i=0; i < NFEII; i++ )
107 /* optical depths for Feii dest of lya when large feii not used */
108 Fe2TauLte[i] += feopc[i]*(realnum)radius.drad_x_fillfac;
109 return;
110}
111
114{
115 DEBUG_ENTRY( "atoms_fe2ovr()" );
116
117 long int i;
118
119 static long int nZoneEval;
120
121 double Fe2Partn,
122 displa,
123 hopc,
124 rate,
125 weight;
126
127 static double BigFeWidth,
128 BigHWidth;
129
130 /* wavelength of Lya in Angstroms */
131
132 /* compute efficiency of FeII emission overlapping with Ly-alpha
133 * implemented with Fred Hamann
134 *
135 * make Ly-a width monotonically increasing to avoid oscillation
136 * in deep regions of x-ray ionized clouds.
137 *
138 * do nothing if large FeII atom is used */
139 if( FeII.lgFeIILargeOn )
140 {
141 return;
142 }
143
144 if( nzone <= 1 )
145 {
146 BigHWidth = hydro.HLineWidth;
147 BigFeWidth = GetDopplerWidth(dense.AtomicWeight[ipIRON]);
148 nZoneEval = nzone;
149 }
150
151 /* do not do pumping if no population,line is thin, or turned off */
152 if( (dense.xIonDense[ipIRON][1] <= 0. || !FeII.lgLyaPumpOn) ||
153 hydro.HLineWidth <= 0. )
154 {
155 Fe2Partn = 0.;
156 hydro.dstfe2lya = 0.;
157
158 for( i=0; i < NFEII; i++ )
159 {
160 Fe2PopLte[i] = 0.;
161 }
162 return;
163 }
164
165 /* only evaluate this one time per zone to avoid oscillations
166 * deep in x-ray ionized clouds */
167 if( nZoneEval == nzone && nzone > 1 )
168 {
169 return;
170 }
171
172 BigHWidth = MAX2(BigHWidth,(double)hydro.HLineWidth);
173 BigFeWidth = MAX2(BigFeWidth,(double)GetDopplerWidth(dense.AtomicWeight[ipIRON]) );
174 nZoneEval = nzone;
175
176 /* check that data is linked in */
177 ASSERT( fe2lam[0] > 0. );
178
179 rate = 0.;
180 Fe2Partn = fe2par(phycon.te);
181 for( i=0; i < NFEII; i++ )
182 {
183 /* this is displacement from line center in units of Lya width */
184 displa = fabs(fe2lam[i]-WLAL)/WLAL*3e10/BigHWidth;
185 if( displa < 1.333 )
186 {
187 /* have variable weighting factor depending on distance away
188 * this comes form the Verner's fits to Adam's results */
189 weight = ( displa <= 0.66666 ) ? 1. : MAX2(0.,1.-(displa-0.666666)/0.66666);
190
191 Fe2PopLte[i] = (realnum)(fe2gs[i]/Fe2Partn*rfield.ContBoltz[ipfe2[i]-1]*
192 dense.xIonDense[ipIRON][1]);
193
194 feopc[i] = (realnum)(Fe2PopLte[i]*fe2osc[i]*0.0150*(fe2lam[i]*1e-8)/BigFeWidth);
195
196 /* Ly-alpha line-center opacity */
197 if( iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop() > 0. )
198 {
199 hopc = 7.60e-8*iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop()/
200 GetDopplerWidth(dense.AtomicWeight[ipHYDROGEN]);
201 }
202 else
203 {
204 hopc = 7.60e-8*dense.xIonDense[ipHYDROGEN][0]/GetDopplerWidth(dense.AtomicWeight[ipHYDROGEN]);
205 }
206 rate += (feopc[i]/SDIV((feopc[i] + hopc)))*(BigFeWidth/GetDopplerWidth(dense.AtomicWeight[ipHYDROGEN]))*
207 (1. - 1./(1. + 1.6*Fe2TauLte[i]))*weight;
208 }
209 }
210
211 /* dstfe2lya is total Lya deexcitation probability due to line overlap */
212 hydro.dstfe2lya = (realnum)rate;
213 return;
214}
215
217double t_fe2ovr_la::fe2par(double te)
218{
219 DEBUG_ENTRY( "fe2par()" );
220
221 double fe2par_v;
222
223 /* function to evaluate partition function for FeII
224 *
225 * Temperature (K) */
226
227 if( te <= fe2pt[0] )
228 fe2par_v = fe2pf[0];
229 else if( te >= fe2pt[NFE2PR-1] )
230 fe2par_v = fe2pf[NFE2PR-1];
231 else
232 {
233 long i = hunt_bisect(fe2pt,NFE2PR,te);
234 double slope = (fe2pf[i+1] - fe2pf[i])/(fe2pt[i+1] - fe2pt[i]);
235 double du = slope*(te - fe2pt[i]);
236 fe2par_v = fe2pf[i] + du;
237 }
238 return( fe2par_v );
239}
const double WLAL
const int NFE2PR
Definition atomfeii.h:270
t_FeII FeII
Definition atomfeii.cpp:5
const int NFEII
Definition atomfeii.h:268
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
const int ipIRON
Definition cddefines.h:330
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long int ipfe2[NFEII]
Definition atomfeii.h:283
realnum fe2gs[NFEII]
Definition atomfeii.h:281
realnum fe2osc[NFEII]
Definition atomfeii.h:279
void atoms_fe2ovr(void)
void zero_opacity()
double fe2pf[NFE2PR]
Definition atomfeii.h:291
realnum fe2enr[NFEII]
Definition atomfeii.h:280
realnum feopc[NFEII]
Definition atomfeii.h:286
realnum fe2lam[NFEII]
Definition atomfeii.h:278
realnum Fe2TauLte[NFEII]
Definition atomfeii.h:287
double fe2pt[NFE2PR]
Definition atomfeii.h:290
void init_pointers()
realnum Fe2PopLte[NFEII]
Definition atomfeii.h:288
double fe2par(double te)
long ipoint(double energy_ryd)
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
t_dense dense
Definition dense.cpp:24
realnum GetDopplerWidth(realnum massAMU)
t_hydro hydro
Definition hydrogenic.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH1s
Definition iso.h:27
const int ipH_LIKE
Definition iso.h:62
t_opac opac
Definition opacity.cpp:5
t_phycon phycon
Definition phycon.cpp:6
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
long hunt_bisect(const T x[], long n, T xval)
Definition thirdparty.h:270