cloudy trunk
Loading...
Searching...
No Matches
mole.h
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3
4#ifndef MOLE_H_
5#define MOLE_H_
6
7/* mole.h */
8
9#include "count_ptr.h"
10#include "elementnames.h"
11#include "transition.h"
12
13#define SMALLABUND 1e-24
14
16
17class chem_atom;
18
20 explicit chem_element(); // Do not implement
21 chem_element &operator=(const chem_element&); // Do not implement
22public:
23 explicit chem_element(int Z, const char*label) : Z(Z), label(label)
24 {}
25 ~chem_element() throw()
26 {}
27 const int Z;
28 const string label;
29 map<int, count_ptr<chem_atom> > isotopes;
30 //(first -> Atomic A; second -> chem_atom )
31 //(first -> -1 for bogus isotope, i.e. where no
32 // isotopes have been explicitly defined)
33};
34
35typedef map<int, count_ptr<chem_atom> >::iterator isotopes_i;
36
37class chem_atom {
38public:
39 // Link back to basic element for convenience -- not a count_ptr
40 // as this would lead to a reference cycle (and hence the
41 // destructors never being called). Many-to-one relation suggests
42 // that the weak link should be this way around.
44 int A; /* mass number */
45 vector<int> ipMl; /* Atom and ion species in molecule arrays */
46 realnum mass_amu; /* mass of isotope in AMU */
47 double frac; /* fraction of element in this isotope */
48
49 bool lgMeanAbundance( void ) const
50 {
51 return ( A < 0 );
52 }
53
54 /* Chemical symbols for elements */
55 string label( void ) const
56 {
57 if( lgMeanAbundance() )
58 return el->label;
59 else if( el->Z==1 && A==2 )
60 {
61 // Deuterium is a special case
62 return "D\0";
63 }
64 else
65 {
66 char str[4];
67 sprintf(str,"^%d",A);
68 return ( str + el->label );
69 }
70 }
71
72 int compare ( const chem_atom &b ) const
73 {
74 // sort by proton number first
75 if ( el->Z < b.el->Z )
76 return -1;
77 else if ( el->Z > b.el->Z )
78 return 1;
79
80 if (mass_amu < b.mass_amu)
81 return -1;
82 else if (mass_amu > b.mass_amu)
83 return 1;
84 else if (A < b.A )
85 return -1;
86 else
87 return 0;
88 }
89};
90inline bool operator< (const chem_atom &a, const chem_atom &b)
91{
92 return a.compare(b) < 0;
93}
94inline bool operator> (const chem_atom &a, const chem_atom &b)
95{
96 return a.compare(b) > 0;
97}
98inline bool operator<= (const chem_atom &a, const chem_atom &b)
99{
100 return a.compare(b) <= 0;
101}
102inline bool operator>= (const chem_atom &a, const chem_atom &b)
103{
104 return a.compare(b) >= 0;
105}
106inline bool operator== (const chem_atom &a, const chem_atom &b)
107{
108 return a.compare(b) == 0;
109}
110inline bool operator!= (const chem_atom &a, const chem_atom &b)
111{
112 return !(a == b);
113}
114
115typedef vector< count_ptr<chem_atom> > ChemAtomList;
119extern chem_atom *null_atom;
120
122{
123public:
125 const count_ptr<chem_atom>& b) const
126 {
127 return *a < *b;
128 }
129};
130
131/* Structure containing molecule data, initially only CO */
132class molecule {
133public:
134 typedef map<const count_ptr<chem_atom>, int,
136
139 bool isEnabled; /* Is it enabled? */
140
141 /* Species physical data */
142 string label;
144 int charge;
145 bool lgExcit;
147 int n_nuclei(void) const
148 {
149 int num = 0;
150 for (nAtomsMap::const_iterator el = nAtom.begin();
151 el != nAtom.end(); ++el)
152 {
153 num += el->second;
154 }
155 return num;
156 }
157 bool isMonatomic(void) const
158 {
159 if (nAtom.size() == 1 && nAtom.begin()->second == 1)
160 return true;
161 return false;
162 }
163
166
167 /* Parameters as computational object */
170
171 chem_atom *heavyAtom(void) //const
172 {
173 for( nAtomsMap::reverse_iterator it=nAtom.rbegin(); it!=nAtom.rend(); ++it )
174 {
175 if (0 != it->second )
176 {
177 return it->first.get_ptr();
178 }
179 }
180 return null_atom;
181 }
182
183 int compare(const molecule &mol2) const
184 {
185 nAtomsMap::const_reverse_iterator it1, it2;
186
187 for( it1 = nAtom.rbegin(), it2 = mol2.nAtom.rbegin();
188 it1 != nAtom.rend() && it2 != mol2.nAtom.rend(); ++it1, ++it2 )
189 {
190 if( *(it1->first) > *(it2->first) )
191 return 1;
192 else if( *(it1->first) < *(it2->first) )
193 return -1;
194 else if( it1->second > it2->second)
195 return 1;
196 else if( it1->second < it2->second)
197 return -1;
198 }
199
200 if( it1 != nAtom.rend() && it2 == mol2.nAtom.rend() )
201 return 1;
202 else if( it1 == nAtom.rend() && it2 != mol2.nAtom.rend() )
203 return -1;
204 else
205 ASSERT( it1 == nAtom.rend() && it2 == mol2.nAtom.rend() );
206
207 // sort by label if falls through to here
208 return ( label.compare(mol2.label) );
209
210 }
211};
212
213/* iterators on nAtom */
214typedef molecule::nAtomsMap::iterator nAtoms_i;
215typedef molecule::nAtomsMap::reverse_iterator nAtoms_ri;
216typedef molecule::nAtomsMap::const_reverse_iterator nAtoms_cri;
217
219extern void mole_drive(void);
220
222extern void mole_create_react(void);
223
224class mole_reaction;
225
226mole_reaction *mole_findrate_s(const char buf[]);
227
228extern void mole_print_species_reactions( molecule *speciesToPrint );
229
230extern molecule *null_mole;
231
232extern molecule *findspecies(const char buf[]);
233
258
265
267
268public:
269 void init(void);
270
271 void make_species(void);
272
274 void zero(void);
275
278
281
284
287
290
295
300
305
309
310 // flag saying whether to model isotopes (and isotopologues) of a given element
311 vector<bool> lgTreatIsotopes;
312
315
316 typedef vector<count_ptr<molecule> > MoleculeList;
318
319 static void sort(MoleculeList::iterator start,
320 MoleculeList::iterator end);
321 static void sort(molecule **start, molecule **end);
322};
323
325
326class molezone {
327public:
329 {
330 init();
331 }
332 void init (void)
333 {
334 location = NULL;
335 levels = NULL;
336 lines = NULL;
337 zero();
338 }
339 void zero (void)
340 {
341 src = 0.;
342 snk = 0.;
343 den = 0.;
344 column = 0.;
345 nAtomLim = -1;
346 xFracLim = 0.;
347 column_old = 0.;
348 }
349 double *location;
350
352 double src, snk;
353
356
357 /* Current zone data */
358 double den;
362
363 /* Historical solution data */
365};
366
367extern molezone *null_molezone;
368
370{
371public:
372 void set_location( long nelem, long ion, double *dense );
373 void set_isotope_abundances( void );
374 double sink_rate_tot(const char chSpecies[]) const;
375 double sink_rate_tot(const molecule* const sp) const;
376 double sink_rate(const molecule* const sp, const mole_reaction& rate) const;
377 double sink_rate(const molecule* const sp, const char buf[]) const;
378 double source_rate_tot(const char chSpecies[]) const;
379 double source_rate_tot(const molecule* const sp) const;
382 double dissoc_rate(const char chSpecies[]) const;
383 double chem_heat(void) const;
384 double findrk(const char buf[]) const;
385 double findrate(const char buf[]) const;
386
388
390 double elec;
391
394 double **source , **sink;
395
396 realnum ***xMoleChTrRate;/***[LIMELM][LIMELM+1][LIMELM+1];*/
397
398 valarray<class molezone> species;
399
400 vector<double> reaction_rks;
401 vector<double> old_reaction_rks;
403};
404
405extern t_mole_local mole;
406
407extern molezone *findspecieslocal(const char buf[]);
408
409extern void mole_punch(FILE *punit, const char speciesname[], const char args[], bool lgHeader, bool lgData, double depth);
410
411extern void total_molecule_elems(realnum total[LIMELM]);
412extern void total_molecule_deut(realnum &total);
413
414extern realnum total_molecules(void);
415
417
418extern void mole_make_list(void);
419extern void mole_make_groups(void);
420
422
423bool lgDifferByExcitation( const molecule &mol1, const molecule &mol2 );
424
425extern void mole_update_species_cache(void);
426
427void mole_update_sources(void);
428
429void mole_rk_bigchange(void);
430
433 vector< int >& numAtoms,
434 string atom_old,
435 string atom_new,
436 string embellishments,
437 vector<string>& newLabels );
438
439bool parse_species_label( const char label[], ChemAtomList &atomsLeftToRight, vector<int> &numAtoms, string &embellishments );
440bool parse_species_label( const char mylab[], ChemAtomList &atomsLeftToRight, vector<int> &numAtoms, string &embellishments,
441 bool &lgExcit, int &charge, bool &lgGas_Phase );
442
443#endif /* MOLE_H_ */
444
t_atoms atoms
Definition atoms.cpp:5
#define ASSERT(exp)
Definition cddefines.h:578
const int LIMELM
Definition cddefines.h:258
float realnum
Definition cddefines.h:103
int compare(const chem_atom &b) const
Definition mole.h:72
string label(void) const
Definition mole.h:55
int A
Definition mole.h:44
bool lgMeanAbundance(void) const
Definition mole.h:49
vector< int > ipMl
Definition mole.h:45
double frac
Definition mole.h:47
chem_element * el
Definition mole.h:43
realnum mass_amu
Definition mole.h:46
map< int, count_ptr< chem_atom > > isotopes
Definition mole.h:29
chem_element(int Z, const char *label)
Definition mole.h:23
~chem_element()
Definition mole.h:25
const int Z
Definition mole.h:27
chem_element & operator=(const chem_element &)
const string label
Definition mole.h:28
bool operator()(const count_ptr< chem_atom > &a, const count_ptr< chem_atom > &b) const
Definition mole.h:124
int parentIndex
Definition mole.h:138
chem_atom * heavyAtom(void)
Definition mole.h:171
realnum form_enthalpy
Definition mole.h:164
bool lgGas_Phase
Definition mole.h:146
map< const count_ptr< chem_atom >, int, element_pointer_value_less > nAtomsMap
Definition mole.h:135
int charge
Definition mole.h:144
bool isEnabled
Definition mole.h:139
string parentLabel
Definition mole.h:137
bool lgExcit
Definition mole.h:145
nAtomsMap nAtom
Definition mole.h:143
realnum mole_mass
Definition mole.h:165
string label
Definition mole.h:142
int compare(const molecule &mol2) const
Definition mole.h:183
enum mole_state state
Definition mole.h:168
int groupnum
Definition mole.h:169
int n_nuclei(void) const
Definition mole.h:147
bool isMonatomic(void) const
Definition mole.h:157
int index
Definition mole.h:169
qList * levels
Definition mole.h:354
realnum column
Definition mole.h:359
double snk
Definition mole.h:352
double * location
Definition mole.h:349
void init(void)
Definition mole.h:332
double den
Definition mole.h:358
double src
Definition mole.h:352
int nAtomLim
Definition mole.h:360
void zero(void)
Definition mole.h:339
realnum xFracLim
Definition mole.h:361
realnum column_old
Definition mole.h:364
TransitionList * lines
Definition mole.h:355
molezone()
Definition mole.h:328
bool lgLeidenHack
Definition mole.h:286
bool lgNoHeavyMole
Definition mole.h:280
void zero(void)
Definition mole.cpp:36
bool lgH2Ozer
Definition mole.h:283
bool lgFederman
Definition mole.h:288
bool lgNeutrals
Definition mole.h:304
void make_species(void)
int num_compacted
Definition mole.h:314
vector< count_ptr< molecule > > MoleculeList
Definition mole.h:316
bool lgStancil
Definition mole.h:289
MoleculeList list
Definition mole.h:317
int num_calc
Definition mole.h:314
bool lgNoMole
Definition mole.h:277
bool lgNonEquilChem
Definition mole.h:294
void init(void)
Definition mole.cpp:11
static void sort(MoleculeList::iterator start, MoleculeList::iterator end)
vector< bool > lgTreatIsotopes
Definition mole.h:311
int num_total
Definition mole.h:314
bool lgGrain_mole_deplete
Definition mole.h:308
bool lgProtElim
Definition mole.h:299
double dissoc_rate(const char chSpecies[]) const
double ** source
Definition mole.h:394
vector< double > reaction_rks
Definition mole.h:400
double sink_rate(const molecule *const sp, const mole_reaction &rate) const
double source_rate_tot(const char chSpecies[]) const
double ** sink
Definition mole.h:394
double grain_density
Definition mole.h:387
double findrate(const char buf[]) const
void set_location(long nelem, long ion, double *dense)
double findrk(const char buf[]) const
vector< double > old_reaction_rks
Definition mole.h:401
double grain_saturation
Definition mole.h:387
double grain_area
Definition mole.h:387
double chem_heat(void) const
long old_zone
Definition mole.h:402
void set_isotope_abundances(void)
double elec
Definition mole.h:390
valarray< class molezone > species
Definition mole.h:398
double sink_rate_tot(const char chSpecies[]) const
realnum *** xMoleChTrRate
Definition mole.h:396
t_dense dense
Definition dense.cpp:24
t_mole_global mole_global
Definition mole.cpp:6
t_mole_local mole
Definition mole.cpp:7
molecule::nAtomsMap::reverse_iterator nAtoms_ri
Definition mole.h:215
void mole_create_react(void)
realnum total_molecules_gasphase(void)
bool operator<(const chem_atom &a, const chem_atom &b)
Definition mole.h:90
chem_element * null_element
ChemAtomList atom_list
molecule::nAtomsMap::const_reverse_iterator nAtoms_cri
Definition mole.h:216
void mole_make_groups(void)
chem_atom * null_atom
realnum total_molecules(void)
void mole_print_species_reactions(molecule *speciesToPrint)
void create_isotopologues_one(ChemAtomList &atoms, vector< int > &numAtoms, string atom_old, string atom_new, string embellishments, vector< string > &newLabels)
molezone * findspecieslocal(const char buf[])
void mole_punch(FILE *punit, const char speciesname[], const char args[], bool lgHeader, bool lgData, double depth)
void mole_cmp_num_in_out_reactions(void)
void mole_rk_bigchange(void)
bool operator>=(const chem_atom &a, const chem_atom &b)
Definition mole.h:102
void mole_drive(void)
mole_reaction * mole_findrate_s(const char buf[])
void mole_update_sources(void)
mole_state
Definition mole.h:15
@ MOLE_NULL
Definition mole.h:15
@ MOLE_ACTIVE
Definition mole.h:15
@ MOLE_PASSIVE
Definition mole.h:15
void mole_update_species_cache(void)
void mole_make_list(void)
bool operator==(const chem_atom &a, const chem_atom &b)
Definition mole.h:106
ChemAtomList unresolved_atom_list
molecule::nAtomsMap::iterator nAtoms_i
Definition mole.h:214
bool operator!=(const chem_atom &a, const chem_atom &b)
Definition mole.h:110
bool lgDifferByExcitation(const molecule &mol1, const molecule &mol2)
bool operator>(const chem_atom &a, const chem_atom &b)
Definition mole.h:94
bool operator<=(const chem_atom &a, const chem_atom &b)
Definition mole.h:98
bool parse_species_label(const char label[], ChemAtomList &atomsLeftToRight, vector< int > &numAtoms, string &embellishments)
molecule * findspecies(const char buf[])
molezone * null_molezone
map< int, count_ptr< chem_atom > >::iterator isotopes_i
Definition mole.h:35
void total_molecule_elems(realnum total[LIMELM])
void total_molecule_deut(realnum &total)
molecule * null_mole
vector< count_ptr< chem_atom > > ChemAtomList
Definition mole.h:115
STATIC void start(long, realnum[], realnum[], long, long[], realnum *, int *)
char ** chSpecies
Definition taulines.cpp:13