cloudy trunk
Loading...
Searching...
No Matches
mean.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/*MeanInc increment mean ionization fractions and temperatures over computed structure */
4/*MeanZero zero mean of ionization fractions / temperatures arrays */
5/*MeanIon derive mean ionization fractions or mean temperatures for some element */
6#include "cddefines.h"
7#include "physconst.h"
8#include "radius.h"
9#include "dense.h"
10#include "hyperfine.h"
11#include "magnetic.h"
12#include "hmi.h"
13#include "phycon.h"
14#include "geometry.h"
15#include "mean.h"
16
18
20{
21 DEBUG_ENTRY( "t_mean()" );
22
23 // dim 3 for average over radius, area, and volume
24 mean.xIonMean.reserve(3);
25 for( int j=0; j < 3; ++j )
26 {
27 // compute average for every element
28 mean.xIonMean.reserve(j,LIMELM);
29 for( int nelem=0; nelem < LIMELM; ++nelem )
30 {
31 int limit = max(3,nelem+2);
32 // compute average for every ionization stage (incl. H2 for hydrogen)
33 mean.xIonMean.reserve(j,nelem,limit);
34 for( int ion=0; ion < limit; ++ion )
35 // dim 2 for storing Sum(quant*norm) and Sum(norm)
36 mean.xIonMean.reserve(j,nelem,ion,2);
37 }
38 }
39 mean.xIonMean.alloc();
40 mean.xIonEdenMean.alloc( mean.xIonMean.clone() );
41 mean.TempIonMean.alloc( mean.xIonMean.clone() );
42 mean.TempIonEdenMean.alloc( mean.xIonMean.clone() );
43 mean.TempB_HarMean.alloc(3,2);
44 mean.TempHarMean.alloc(3,2);
45 mean.TempH_21cmSpinMean.alloc(3,2);
46 mean.TempMean.alloc(3,2);
47 mean.TempEdenMean.alloc(3,2);
48}
49
50/*MeanZero zero mean of ionization fractions / temperatures arrays */
52{
53 DEBUG_ENTRY( "t_mean::MeanZero()" );
54
55 /* MeanZero is called at start of calculation by zero, and by
56 * startenditer to initialize the variables */
57
58 mean.xIonMean.zero();
59 mean.xIonEdenMean.zero();
60 mean.TempIonMean.zero();
61 mean.TempIonEdenMean.zero();
62 mean.TempB_HarMean.zero();
63 mean.TempHarMean.zero();
64 mean.TempH_21cmSpinMean.zero();
65 mean.TempMean.zero();
66 mean.TempEdenMean.zero();
67
68 return;
69}
70
71/*MeanInc increment mean ionization fractions and temperatures over computed structure */
73{
74 /* this routine is called by radius_increment during the calculation to update the
75 * sums that will become the rad, area, and vol weighted mean ionizations */
76
77 DEBUG_ENTRY( "t_mean::MeanInc()" );
78
79 /* Jacobian used in the integrals below */
80 double Jac[3] = { radius.drad_x_fillfac, radius.darea_x_fillfac, radius.dVeffVol };
81
82 for( int d=0; d < 3; ++d )
83 {
84 double dc;
85 for( int nelem=0; nelem < LIMELM; nelem++ )
86 {
87 int limit = max(3,nelem+2);
88 /* this normalizes xIonMean and xIonEdenMean,
89 * use gas_phase which includes molecular parts */
90 double norm = dense.gas_phase[nelem]*Jac[d];
91 for( int ion=0; ion < limit; ion++ )
92 {
93 if( nelem == ipHYDROGEN && ion == 2 )
94 /* hydrogen is special case since must include H2,
95 * and must mult by 2 to conserve total H - will need to div
96 * by two if column density ever used */
97 dc = 2.*hmi.H2_total*Jac[d];
98 else
99 dc = dense.xIonDense[nelem][ion]*Jac[d];
100 mean.xIonMean[d][nelem][ion][0] += dc;
101 mean.xIonMean[d][nelem][ion][1] += norm;
102 mean.TempIonMean[d][nelem][ion][0] += dc*phycon.te;
103 mean.TempIonMean[d][nelem][ion][1] += dc;
104
105 dc *= dense.eden;
106 mean.xIonEdenMean[d][nelem][ion][0] += dc;
107 mean.xIonEdenMean[d][nelem][ion][1] += norm*dense.eden;
108 mean.TempIonEdenMean[d][nelem][ion][0] += dc*phycon.te;
109 mean.TempIonEdenMean[d][nelem][ion][1] += dc;
110 }
111 }
112
113 /* =============================================================
114 *
115 * these are some special quantities for the mean
116 */
117
118 /* used to get magnetic field weighted wrt harmonic mean spin temperature,
119 * * H0 - as measured by 21cm temperature */
120 dc = ( hyperfine.Tspin21cm > SMALLFLOAT ) ? dense.xIonDense[ipHYDROGEN][0]*Jac[d]/phycon.te : 0.;
121 /* mean magnetic field weighted wrt 21 cm opacity, n(H0)/T */
122 mean.TempB_HarMean[d][0] += sqrt(fabs(magnetic.pressure)*PI8) * dc;
123 /* this assumes that field is tangled */
124 mean.TempB_HarMean[d][1] += dc;
125
126 /* used to get harmonic mean temperature with respect to H,
127 * for comparison with 21cm temperature */
128 dc = dense.xIonDense[ipHYDROGEN][0]*Jac[d];
129 /* harmonic mean gas kinetic temp, for comparison with 21 cm obs */
130 /*HEADS UP - next are the inverse of equation 3 of
131 * >>refer Tspin 21 cm Abel, N.P., Brogan, C.L., Ferland, G.J., O'Dell, C.R.,
132 * >>refercon Shaw, G. & Troland, T.H. 2004, ApJ, 609, 247 */
133 mean.TempHarMean[d][0] += dc;
134 mean.TempHarMean[d][1] += dc/phycon.te;
135
136 /* harmonic mean of computed 21 cm spin temperature - this is what 21 cm actually measures */
137 mean.TempH_21cmSpinMean[d][0] += dc;
138 mean.TempH_21cmSpinMean[d][1] += dc / SDIV( hyperfine.Tspin21cm );
139
140 dc = Jac[d];
141 mean.TempMean[d][0] += dc*phycon.te;
142 mean.TempMean[d][1] += dc;
143
144 dc = Jac[d]*dense.eden;
145 mean.TempEdenMean[d][0] += dc*phycon.te;
146 mean.TempEdenMean[d][1] += dc;
147 }
148 return;
149}
150
151/*MeanIon derive mean ionization fractions or mean temperatures for some element */
153 /* either 'i' or 't' for ionization or temperature */
154 char chType,
155 /* atomic number on C scale */
156 long int nelem,
157 /* type of average: 0=radius, 1=area, 2=volume */
158 long int dim,
159 /* this will say how many ion stages in arlog have non-zero values */
160 long int *n,
161 /* array of values, log both cases,
162 * hydrogen is special case since [2] will be H2 */
163 realnum arlog[],
164 /* true, include electron density, false do not*/
165 bool lgDensity ) const
166{
167 DEBUG_ENTRY( "t_mean::MeanIon()" );
168
169 /* limit is on C scale, such that ion < limit */
170 int limit = max( 3, nelem+2 );
171
172 /* fills in array arlog with log of ionization fractions
173 * n is set to number of non-zero abundances
174 * n set to 0 if element turned off
175 *
176 * first call MeanZero to zero out arrays
177 * next call MeanInc in zone calc to enter ionziation fractions or temperature
178 * finally this routine computes actual means
179 * */
180 if( !dense.lgElmtOn[nelem] )
181 {
182 /* no ionization for this element */
183 for( int ion=0; ion < limit; ion++ )
184 arlog[ion] = -30.f;
185 *n = 0;
186 return;
187 }
188
189 /* set high ion stages, with zero abundances, to -30 */
190 *n = limit;
191 while( *n > 0 && mean.xIonMean[0][nelem][*n-1][0] <= 0. )
192 {
193 arlog[*n-1] = -30.f;
194 --*n;
195 }
196
197 double meanv, normv;
198 for( int ion=0; ion < *n; ion++ )
199 {
200 /* mean ionization */
201 if( chType == 'i' )
202 {
203 if( lgDensity )
204 {
205 meanv = mean.xIonEdenMean[dim][nelem][ion][0];
206 normv = mean.xIonEdenMean[dim][nelem][ion][1];
207 }
208 else
209 {
210 meanv = mean.xIonMean[dim][nelem][ion][0];
211 normv = mean.xIonMean[dim][nelem][ion][1];
212 }
213 arlog[ion] = ( meanv > 0. ) ? (realnum)log10(max(1e-30,meanv/normv)) : -30.f;
214 }
215 /* mean temperature */
216 else if( chType == 't' )
217 {
218 if( lgDensity )
219 {
220 meanv = mean.TempIonEdenMean[dim][nelem][ion][0];
221 normv = mean.TempIonEdenMean[dim][nelem][ion][1];
222 }
223 else
224 {
225 meanv = mean.TempIonMean[dim][nelem][ion][0];
226 normv = mean.TempIonMean[dim][nelem][ion][1];
227 }
228 arlog[ion] = ( normv > SMALLFLOAT ) ? (realnum)log10(max(1e-30,meanv/normv)) : -30.f;
229 }
230 else
231 {
232 fprintf(ioQQQ," MeanIon called with insane job: %c \n",chType);
234 }
235 }
236 return;
237}
FILE * ioQQQ
Definition cddefines.cpp:7
const int LIMELM
Definition cddefines.h:258
float realnum
Definition cddefines.h:103
long max(int a, long b)
Definition cddefines.h:775
NORETURN void TotalInsanity(void)
Definition service.cpp:886
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
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
t_hmi hmi
Definition hmi.cpp:5
t_hyperfine hyperfine
Definition hyperfine.cpp:5
t_magnetic magnetic
Definition magnetic.cpp:17
t_mean mean
Definition mean.cpp:17
t_mean mean
Definition mean.cpp:17
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double PI8
Definition physconst.h:38
t_radius radius
Definition radius.cpp:5
Definition mean.h:10
void MeanIon(char chType, long nelem, long dim, long *n, realnum arlog[], bool lgDensity) const
Definition mean.cpp:152
t_mean()
Definition mean.cpp:19
void MeanZero()
Definition mean.cpp:51
void MeanInc()
Definition mean.cpp:72