cloudy
trunk
Main Page
Related Pages
Data Structures
Files
File List
Globals
All
Data Structures
Files
Functions
Variables
Typedefs
Enumerations
Enumerator
Friends
Macros
Pages
source
eden_sum.cpp
Go to the documentation of this file.
1
/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2
* others. For conditions of distribution and use see copyright notice in license.txt */
3
/*eden_sum sum free electron density over all species, sets variable erredn.EdenTrue
4
*called by ConvBase - ConvEdenIoniz actually updates the electron density
5
* returns 0 if all is ok, 1 if need to abort calc */
6
#include "
cddefines.h
"
7
#include "
hmi.h
"
8
#include "
trace.h
"
9
#include "
grainvar.h
"
10
#include "
rfield.h
"
11
#include "
mole.h
"
12
#include "
dense.h
"
13
#include "
conv.h
"
14
15
/*eden_sum sum free electron density over all species, sets variable erredn.EdenTrue
16
*called by ConvBase - ConvEdenIoniz actually updates the electron density
17
* returns 0 if all is ok, 1 if need to abort calc */
18
int
eden_sum
(
void
)
19
{
20
long
int
i,
21
ion,
22
nelem;
23
24
realnum
sum_all_ions ,
25
sum_metals ,
26
hmole_eden,
27
eden_ions[
LIMELM
];
28
29
DEBUG_ENTRY
(
"eden_sum()"
);
30
31
/* EdenExtra is normally zero, set with EDEN command, to add extra e- */
32
dense
.
EdenTrue
=
dense
.
EdenExtra
;
33
34
/* sum over all ions */
35
sum_all_ions = 0.;
36
sum_metals = 0.;
37
for
( nelem=
ipHYDROGEN
; nelem <
LIMELM
; nelem++ )
38
{
39
if
( nelem==
ipLITHIUM
)
40
sum_metals = 0.;
41
eden_ions[nelem] =
dense
.
xIonDense
[nelem][1];
42
/* >>chng 02 jul 15, upper limit now includes H-like, so all species
43
* in this one loop (except hydrogen and helium ) */
44
/* >>chng 02 apr 26, upper limit had been nelem, not nelem+1, so H-like
45
* metals were missed - ok for H and He since not part of this loop before */
46
for
( ion=2; ion <= nelem+1; ion++ )
47
{
48
/* >>chng 96 oct 27, save all contributors to electron density */
49
eden_ions[nelem] += ion*
dense
.
xIonDense
[nelem][ion];
50
}
51
52
sum_all_ions += eden_ions[nelem];
53
sum_metals += eden_ions[nelem];
54
}
55
dense
.
EdenTrue
+= sum_all_ions;
56
ASSERT
(
dense
.
EdenTrue
>= 0. );
57
58
/* add on molecules */
59
co
.
comole_eden
= 0.;
60
for
( i=0; i <
mole
.
num_comole_calc
; i++ )
61
{
62
if
(
COmole
[i]->n_nuclei != 1)
63
{
64
co
.
comole_eden
+=
COmole
[i]->
hevmol
*
COmole
[i]->
nElec
;
65
}
66
}
67
68
/* electrons contributed by molecules */
69
dense
.
EdenTrue
+=
co
.
comole_eden
;
70
ASSERT
(
dense
.
EdenTrue
>= 0. );
71
72
/* >>chng 03 nov 28, loop over all H molecules, had just been H- */
73
hmole_eden = 0.;
74
for
(i=0;i<
N_H_MOLEC
;i++)
75
{
76
/* hmi.nElectron is zero for H+ since counted as an ion, -1 for H-, etc */
77
hmole_eden +=
hmi
.
Hmolec
[i]*
hmi
.
nElectron
[i];
78
}
79
/* >>chng 01 jan 08, following logic added to stop H- from forcing
80
* negative electron density during very first search, when H- is badly off */
81
/* >>chng 03 nov 28, add test for lgSearch, had been absent */
82
/* >>chng 04 feb 20, recoil_ion.in shows this important,
83
* update test so prevent neg eden */
84
if
( (-hmole_eden) >
dense
.
EdenTrue
/4. &&
conv
.
lgSearch
)
85
{
86
/*dense.EdenTrue += MIN2(dense.EdenTrue*0.05,hmole_eden);*/
87
dense
.
EdenTrue
*= 0.9;
88
}
89
else
if
( (-hmole_eden) >
dense
.
EdenTrue
)
90
{
91
/* >>chng 04 mar 14, add this branch to prevent
92
* negative electron density. occurred in pdr
93
* with large h2, after hmole failure */
94
fprintf(
ioQQQ
,
" PROBLEM sum eden from hmole too neg, set to limt. EdenTrue:%.3e hmole_eden:%.3e \n"
,
95
dense
.
EdenTrue
,
96
hmole_eden);
97
dense
.
EdenTrue
=
dense
.
EdenTrue
/2.;
98
}
99
else
100
{
101
/* account for electrons on all H-bearing molecules */
102
dense
.
EdenTrue
+= hmole_eden;
103
}
104
ASSERT
(
dense
.
EdenTrue
>= 0. );
105
106
/* this variable is set with the set eden command,
107
* is supposed to override physical electron density */
108
if
(
dense
.
EdenSet
> 0. )
109
{
110
dense
.
EdenTrue
=
dense
.
EdenSet
;
111
dense
.
eden_from_metals
= 1.;
112
}
113
else
114
{
115
/* fraction of electrons from si, s, mg, al */
116
/*dense.eden_from_metals = sum_all_ions / SDIV(dense.EdenTrue);*/
117
dense
.
eden_from_metals
= sum_metals /
SDIV
(
dense
.
EdenTrue
);
118
/*fprintf(ioQQQ," debuggg ionfrac %.2f %.3f\n",fnzone,dense.eden_from_metals );*/
119
}
120
/* dense.eden itself is actually changed in ConvEdenIoniz */
121
122
/* >>chng 00 dec 19, include electrons on grains in total sum */
123
/* >>chng 02 jul 15, only add grain electrons when we have stable soln -
124
* everett.in had very negative grain elec total, so drove total eden negative,
125
* but due to bad inital electron density - do not add grain electrons until we
126
* are close to a solution */
127
/* >>chng 04 sep 26, allow nEdenTrue to become -ve */
128
if
(
dense
.
EdenTrue
+
gv
.
TotalEden
*
gv
.
lgGrainElectrons
>= 0. )
129
{
130
/* gv.lgGrainElectrons - should grain electron source/sink be included in overall electron sum?
131
* default is true, set false with no grain electrons command */
132
dense
.
EdenTrue
+=
gv
.
TotalEden
*
gv
.
lgGrainElectrons
;
133
}
134
else
135
{
136
/* >>chng 05 jan 24, only print if not in search phase */
137
if
( !
conv
.
lgSearch
)
138
fprintf(
ioQQQ
,
139
" PROBLEM eden grain neg limt %.2f eden %.4e new eden bef grn %.4e"
140
"grain eden %.4e set edentrue to %.4e Search?%c\n"
,
141
fnzone
,
142
dense
.
eden
,
143
dense
.
EdenTrue
,
144
gv
.
TotalEden
,
145
dense
.
eden
*0.9,
146
TorF
(
conv
.
lgSearch
));
147
148
/*dense.EdenTrue = dense.eden*0.9;*/
149
dense
.
EdenTrue
+=
gv
.
TotalEden
*
gv
.
lgGrainElectrons
;
150
}
151
152
if
(
trace
.
lgTrace
||
trace
.
lgESOURCE
)
153
{
154
fprintf(
ioQQQ
,
155
" eden_sum zn:%.2f current:%.4e new true:%.4e ions:%.4e comole:%.4e hmole:%.4e grain:%.4e extra:%.4e sum:%.4e LaOTS:%.4e\n"
,
156
fnzone
,
157
dense
.
eden
,
158
dense
.
EdenTrue
,
159
sum_all_ions ,
160
co
.
comole_eden
,
161
hmole_eden ,
162
gv
.
TotalEden
*
gv
.
lgGrainElectrons
,
163
dense
.
EdenExtra
,
164
sum_all_ions +
co
.
comole_eden
+ hmole_eden +
gv
.
TotalEden
*
gv
.
lgGrainElectrons
+
dense
.
EdenExtra
,
165
rfield
.
otslin
[2182] );
166
167
if
(
trace
.
lgNeBug
)
168
{
169
for
(nelem=
ipHYDROGEN
; nelem <
LIMELM
; nelem++)
170
{
171
if
( nelem==0 )
172
{
173
fprintf(
ioQQQ
,
" eden_sum H -Ne:"
);
174
}
175
else
if
( nelem==10 )
176
{
177
fprintf(
ioQQQ
,
"\n eden_sum Na-Ca:"
);
178
}
179
else
if
( nelem==20 )
180
{
181
fprintf(
ioQQQ
,
"\n eden_sum Sc-Zn:"
);
182
}
183
fprintf(
ioQQQ
,
" %.4e"
, eden_ions[nelem] );
184
if
( nelem==29 )
185
{
186
fprintf(
ioQQQ
,
"\n"
);
187
}
188
/*if( nelem==9 || nelem==19 || nelem==29 )
189
fprintf( ioQQQ, "\n " );*/
190
}
191
}
192
}
193
194
/* abort if negative electron density */
195
/*>>chng 04 sep 25, allow negative true elec density so solver can work */
196
if
( 0 &&
dense
.
EdenTrue
< 0. )
197
{
198
fprintf(
ioQQQ
,
"eden_sum finds non-positive electron density.\n"
);
199
fprintf(
ioQQQ
,
200
" eden_sum: EdenTrue to%12.4e fm%12.4e Ne(H):%10.2e Ne(He):%10.2e Ne(C):\n"
,
201
dense
.
EdenTrue
,
202
dense
.
eden
,
203
dense
.
xIonDense
[
ipHYDROGEN
][1],
204
dense
.
xIonDense
[
ipHELIUM
][1] + 2.*
dense
.
xIonDense
[
ipHELIUM
][2] );
205
ShowMe
();
206
cdEXIT
(EXIT_FAILURE);
207
}
208
else
if
(
dense
.
EdenTrue
== 0. )
209
{
210
fprintf(
ioQQQ
,
"\nDISASTER PROBLEM eden_sum finds an electron density of zero, this is unphysical.\n"
);
211
fprintf(
ioQQQ
,
"There appears to be no source of ionization.\n"
);
212
fprintf(
ioQQQ
,
"Consider adding background cosmic rays with COSMIC RAY BACKGROUND and BACKGROUND commands.\n\n"
);
213
lgAbort
=
true
;
214
215
return
1;
216
}
217
218
{
219
/*@-redef@*/
220
enum
{DEBUG_LOC=
false
};
221
/*@+redef@*/
222
if
( DEBUG_LOC )
223
{
224
fprintf(
ioQQQ
,
"esumdebugg\t%li\t%.2e\t%.2e\n"
,
225
nzone
,
226
dense
.
eden
,
dense
.
EdenTrue
);
227
}
228
}
229
230
/* >>chng 05 jan 05, don't let elec den be zero - logs are taken */
231
dense
.
eden
=
MAX2
(
SMALLFLOAT
,
dense
.
eden
);
232
233
return
0;
234
}
Generated for cloudy by
1.8.1.2