00001 /* This file is part of Cloudy and is copyright (C)1978-2010 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*CO_step fills in matrix for heavy elements molecular routines */ 00004 #include "cddefines.h" 00005 #include "mole.h" 00006 #include "mole_co_priv.h" 00007 /* Nick Abel between July and October of 2003 assisted Dr. Ferland in improving the heavy element 00008 * molecular network in Cloudy. Before this routine would predict negative abundances if 00009 * the fraction of carbon in the form of molecules came close to 100%. A reorganizing of 00010 * the reaction network detected several bugs. Treatment of "coupled reactions", 00011 * in which both densities in the reaction rate were being predicted by Cloudy, were also 00012 * added. Due to these improvements, Cloudy can now perform calculations 00013 * where 100% of the carbon is in the form of CO without predicting negative abundances 00014 * 00015 * Additional changes were made in November of 2003 so that our reaction 00016 * network would include all reactions from the TH85 paper. This involved 00017 * adding silicon to the chemical network. Also the reaction rates were 00018 * labeled to make identification with the reaction easier and the matrix 00019 * elements of atomic C, O, and Si are now done in a loop, which makes 00020 * the addition of future chemical species (like N or S) easy. 00021 * */ 00022 /* Robin Williams in August 2006 onwards reorganized the coding to cut down repetitions. 00023 * This isolated several further bugs, and allows a sigificant number of lines of 00024 * code to be eliminated. The balance of S2/S2+ amd ClO/ClO+ seems highly sensitive 00025 * (with small log scale results varying significantly if the order of arithmetic 00026 * operations is changed) -- I suspect this may imply a bug somewhere. 00027 * */ 00028 /*lint -e778 constant expression evaluatess to 0 in operation '-' */ 00029 /*=================================================================*/ 00030 00031 00032 void CO_step(void) 00033 { 00034 long int i, j, n, nt, ratei, ratej; 00035 struct COmole_rate_s *rate; 00036 double rate_tot, rate_deriv[MAXREACTANTS], rated, rk, rate_bval; 00037 00038 DEBUG_ENTRY("CO_step()"); 00039 /* zero out array used for formation rates */ 00040 for( i=0; i < mole.num_comole_calc; i++ ) 00041 { 00042 mole.b[i] = 0.; 00043 } 00044 for( j=0; j < mole.num_comole_tot; j++ ) 00045 { 00046 for( i=0; i < mole.num_comole_calc; i++ ) 00047 { 00048 mole.c[j][i] = 0.; 00049 } 00050 } 00051 00052 00053 /* Call all the routines that set up the matrix 00054 * CO_solve will call this routine, therefore all other matrix elements are 00055 * included here so that, when CO_solve is called, everything is accounted for */ 00056 00057 /* All now cross-validated with new treatment, switching causes only v. minor 00058 * differences in results */ 00059 /* Revised molecular network implementation */ 00060 /* Fills matrix elements for passive species -- these can be used to 00061 derive sources & sinks resulting from this part of the network */ 00062 nt = coreactions.n; 00063 for(n=0;n<nt;n++) 00064 { 00065 rate = coreactions.list[n]; 00066 rk = rate->rk; 00067 for(i=0;i<rate->nrates;i++) 00068 { 00069 rate_deriv[i] = rk; 00070 for(j=0;j<rate->nrates;j++) 00071 { 00072 if(i!=j) 00073 { 00074 rate_deriv[i] *= rate->rate_species[j]->hevmol; 00075 } 00076 } 00077 } 00078 00079 if(rate->nreactants != 1) 00080 { 00081 rate_tot = rate_deriv[0]*rate->rate_species[0]->hevmol; 00082 rate_bval = (rate->nreactants-1)*rate_tot; 00083 for(i=0;i<rate->nreactants;i++) 00084 { 00085 ratei = rate->reactants[i]->index; 00086 mole.b[ratei] -= rate_bval; 00087 } 00088 00089 for(i=0;i<rate->nproducts;i++) 00090 { 00091 ratei = rate->products[i]->index; 00092 mole.b[ratei] += rate_bval; 00093 } 00094 } 00095 00096 for(j=0;j<rate->nrates;j++) 00097 { 00098 ratej = rate->rate_species[j]->index; 00099 rated = rate_deriv[j]; 00100 for(i=0;i<rate->nreactants;i++) 00101 { 00102 mole.c[ratej][rate->reactants[i]->index] -= rated; 00103 } 00104 for(i=0;i<rate->nproducts;i++) 00105 { 00106 mole.c[ratej][rate->products[i]->index] += rated; 00107 } 00108 } 00109 } 00110 00111 }