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