cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 doxygen 1.8.4