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
mole_co_step.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
/*CO_step fills in matrix for heavy elements molecular routines */
4
#include "
cddefines.h
"
5
#include "
mole.h
"
6
#include "
mole_co_priv.h
"
7
/* Nick Abel between July and October of 2003 assisted Dr. Ferland in improving the heavy element
8
* molecular network in Cloudy. Before this routine would predict negative abundances if
9
* the fraction of carbon in the form of molecules came close to 100%. A reorganizing of
10
* the reaction network detected several bugs. Treatment of "coupled reactions",
11
* in which both densities in the reaction rate were being predicted by Cloudy, were also
12
* added. Due to these improvements, Cloudy can now perform calculations
13
* where 100% of the carbon is in the form of CO without predicting negative abundances
14
*
15
* Additional changes were made in November of 2003 so that our reaction
16
* network would include all reactions from the TH85 paper. This involved
17
* adding silicon to the chemical network. Also the reaction rates were
18
* labeled to make identification with the reaction easier and the matrix
19
* elements of atomic C, O, and Si are now done in a loop, which makes
20
* the addition of future chemical species (like N or S) easy.
21
* */
22
/* Robin Williams in August 2006 onwards reorganized the coding to cut down repetitions.
23
* This isolated several further bugs, and allows a sigificant number of lines of
24
* code to be eliminated. The balance of S2/S2+ amd ClO/ClO+ seems highly sensitive
25
* (with small log scale results varying significantly if the order of arithmetic
26
* operations is changed) -- I suspect this may imply a bug somewhere.
27
* */
28
/*lint -e778 constant expression evaluatess to 0 in operation '-' */
29
/*=================================================================*/
30
31
32
void
CO_step
(
void
)
33
{
34
long
int
i, j, n, nt, ratei, ratej;
35
struct
COmole_rate_s
*rate;
36
double
rate_tot, rate_deriv[
MAXREACTANTS
], rated,
rk
, rate_bval;
37
38
DEBUG_ENTRY
(
"CO_step()"
);
39
/* zero out array used for formation rates */
40
for
( i=0; i <
mole
.
num_comole_calc
; i++ )
41
{
42
mole
.
b
[i] = 0.;
43
}
44
for
( j=0; j <
mole
.
num_comole_tot
; j++ )
45
{
46
for
( i=0; i <
mole
.
num_comole_calc
; i++ )
47
{
48
mole
.
c
[j][i] = 0.;
49
}
50
}
51
52
53
/* Call all the routines that set up the matrix
54
* CO_solve will call this routine, therefore all other matrix elements are
55
* included here so that, when CO_solve is called, everything is accounted for */
56
57
/* All now cross-validated with new treatment, switching causes only v. minor
58
* differences in results */
59
/* Revised molecular network implementation */
60
/* Fills matrix elements for passive species -- these can be used to
61
derive sources & sinks resulting from this part of the network */
62
nt =
coreactions
.
n
;
63
for
(n=0;n<nt;n++)
64
{
65
rate =
coreactions
.
list
[n];
66
rk = rate->
rk
;
67
for
(i=0;i<rate->
nrates
;i++)
68
{
69
rate_deriv[i] =
rk
;
70
for
(j=0;j<rate->
nrates
;j++)
71
{
72
if
(i!=j)
73
{
74
rate_deriv[i] *= rate->
rate_species
[j]->
hevmol
;
75
}
76
}
77
}
78
79
if
(rate->
nreactants
!= 1)
80
{
81
rate_tot = rate_deriv[0]*rate->
rate_species
[0]->
hevmol
;
82
rate_bval = (rate->
nreactants
-1)*rate_tot;
83
for
(i=0;i<rate->
nreactants
;i++)
84
{
85
ratei = rate->
reactants
[i]->
index
;
86
mole
.
b
[ratei] -= rate_bval;
87
}
88
89
for
(i=0;i<rate->
nproducts
;i++)
90
{
91
ratei = rate->
products
[i]->
index
;
92
mole
.
b
[ratei] += rate_bval;
93
}
94
}
95
96
for
(j=0;j<rate->
nrates
;j++)
97
{
98
ratej = rate->
rate_species
[j]->
index
;
99
rated = rate_deriv[j];
100
for
(i=0;i<rate->
nreactants
;i++)
101
{
102
mole
.
c
[ratej][rate->
reactants
[i]->
index
] -= rated;
103
}
104
for
(i=0;i<rate->
nproducts
;i++)
105
{
106
mole
.
c
[ratej][rate->
products
[i]->
index
] += rated;
107
}
108
}
109
}
110
111
}
Generated for cloudy by
1.8.1.2