00001
00002
00003
00004 #ifndef MOLE_H_
00005 #define MOLE_H_
00006
00007
00008
00009 #include "count_ptr.h"
00010 #include "elementnames.h"
00011 #include "transition.h"
00012
00013 #define SMALLABUND 1e-24
00014
00015 enum mole_state {MOLE_NULL, MOLE_PASSIVE, MOLE_ACTIVE};
00016
00017 class chem_atom;
00018
00019 class chem_element {
00020 explicit chem_element();
00021 chem_element &operator=(const chem_element&);
00022 public:
00023 explicit chem_element(int Z, const char*label) : Z(Z), label(label)
00024 {}
00025 ~chem_element() throw()
00026 {}
00027 const int Z;
00028 const string label;
00029 map<int, count_ptr<chem_atom> > isotopes;
00030
00031
00032
00033 };
00034
00035 typedef map<int, count_ptr<chem_atom> >::iterator isotopes_i;
00036
00037 class chem_atom {
00038 public:
00039
00040
00041
00042
00043 chem_element* el;
00044 int A;
00045 vector<int> ipMl;
00046 realnum mass_amu;
00047 double frac;
00048
00049 bool lgMeanAbundance( void ) const
00050 {
00051 return ( A < 0 );
00052 }
00053
00054
00055 string label( void ) const
00056 {
00057 if( lgMeanAbundance() )
00058 return el->label;
00059 else if( el->Z==1 && A==2 )
00060 {
00061
00062 return "D\0";
00063 }
00064 else
00065 {
00066 char str[4];
00067 sprintf(str,"^%d",A);
00068 return ( str + el->label );
00069 }
00070 }
00071
00072 int compare ( const chem_atom &b ) const
00073 {
00074
00075 if ( el->Z < b.el->Z )
00076 return -1;
00077 else if ( el->Z > b.el->Z )
00078 return 1;
00079
00080 if (mass_amu < b.mass_amu)
00081 return -1;
00082 else if (mass_amu > b.mass_amu)
00083 return 1;
00084 else if (A < b.A )
00085 return -1;
00086 else
00087 return 0;
00088 }
00089 };
00090 inline bool operator< (const chem_atom &a, const chem_atom &b)
00091 {
00092 return a.compare(b) < 0;
00093 }
00094 inline bool operator> (const chem_atom &a, const chem_atom &b)
00095 {
00096 return a.compare(b) > 0;
00097 }
00098 inline bool operator<= (const chem_atom &a, const chem_atom &b)
00099 {
00100 return a.compare(b) <= 0;
00101 }
00102 inline bool operator>= (const chem_atom &a, const chem_atom &b)
00103 {
00104 return a.compare(b) >= 0;
00105 }
00106 inline bool operator== (const chem_atom &a, const chem_atom &b)
00107 {
00108 return a.compare(b) == 0;
00109 }
00110 inline bool operator!= (const chem_atom &a, const chem_atom &b)
00111 {
00112 return !(a == b);
00113 }
00114
00115 typedef vector< count_ptr<chem_atom> > ChemAtomList;
00116 extern ChemAtomList atom_list;
00117 extern ChemAtomList unresolved_atom_list;
00118 extern chem_element *null_element;
00119 extern chem_atom *null_atom;
00120
00121 class element_pointer_value_less
00122 {
00123 public:
00124 bool operator()(const count_ptr<chem_atom>& a,
00125 const count_ptr<chem_atom>& b) const
00126 {
00127 return *a < *b;
00128 }
00129 };
00130
00131
00132 class molecule {
00133 public:
00134 typedef map<const count_ptr<chem_atom>, int,
00135 element_pointer_value_less> nAtomsMap;
00136
00137 string parentLabel;
00138 int parentIndex;
00139 bool isEnabled;
00140
00141
00142 string label;
00143 nAtomsMap nAtom;
00144 int charge;
00145 bool lgExcit;
00146 bool lgGas_Phase;
00147 int n_nuclei(void) const
00148 {
00149 int num = 0;
00150 for (nAtomsMap::const_iterator el = nAtom.begin();
00151 el != nAtom.end(); ++el)
00152 {
00153 num += el->second;
00154 }
00155 return num;
00156 }
00157 bool isMonatomic(void) const
00158 {
00159 if (nAtom.size() == 1 && nAtom.begin()->second == 1)
00160 return true;
00161 return false;
00162 }
00163
00164 realnum form_enthalpy;
00165 realnum mole_mass;
00167
00168 enum mole_state state;
00169 int index, groupnum;
00170
00171 chem_atom *heavyAtom(void)
00172 {
00173 for( nAtomsMap::reverse_iterator it=nAtom.rbegin(); it!=nAtom.rend(); ++it )
00174 {
00175 if (0 != it->second )
00176 {
00177 return it->first.get_ptr();
00178 }
00179 }
00180 return null_atom;
00181 }
00182
00183 int compare(const molecule &mol2) const
00184 {
00185 nAtomsMap::const_reverse_iterator it1, it2;
00186
00187 for( it1 = nAtom.rbegin(), it2 = mol2.nAtom.rbegin();
00188 it1 != nAtom.rend() && it2 != mol2.nAtom.rend(); ++it1, ++it2 )
00189 {
00190 if( *(it1->first) > *(it2->first) )
00191 return 1;
00192 else if( *(it1->first) < *(it2->first) )
00193 return -1;
00194 else if( it1->second > it2->second)
00195 return 1;
00196 else if( it1->second < it2->second)
00197 return -1;
00198 }
00199
00200 if( it1 != nAtom.rend() && it2 == mol2.nAtom.rend() )
00201 return 1;
00202 else if( it1 == nAtom.rend() && it2 != mol2.nAtom.rend() )
00203 return -1;
00204 else
00205 ASSERT( it1 == nAtom.rend() && it2 == mol2.nAtom.rend() );
00206
00207
00208 return ( label.compare(mol2.label) );
00209
00210 }
00211 };
00212
00213
00214 typedef molecule::nAtomsMap::iterator nAtoms_i;
00215 typedef molecule::nAtomsMap::reverse_iterator nAtoms_ri;
00216 typedef molecule::nAtomsMap::const_reverse_iterator nAtoms_cri;
00217
00219 extern void mole_drive(void);
00220
00222 extern void mole_create_react(void);
00223
00224 class mole_reaction;
00225
00226 mole_reaction *mole_findrate_s(const char buf[]);
00227
00228 extern void mole_print_species_reactions( molecule *speciesToPrint );
00229
00230 extern molecule *null_mole;
00231
00232 extern molecule *findspecies(const char buf[]);
00233
00266 class t_mole_global {
00267
00268 public:
00269 void init(void);
00270
00271 void make_species(void);
00272
00274 void zero(void);
00275
00277 bool lgNoMole;
00278
00280 bool lgNoHeavyMole;
00281
00283 bool lgH2Ozer;
00284
00286 bool lgLeidenHack;
00287
00288 bool lgFederman;
00289 bool lgStancil;
00290
00294 bool lgNonEquilChem;
00295
00299 bool lgProtElim;
00300
00304 bool lgNeutrals;
00305
00308 bool lgGrain_mole_deplete;
00309
00310
00311 vector<bool> lgTreatIsotopes;
00312
00314 int num_total, num_calc, num_compacted;
00315
00316 typedef vector<count_ptr<molecule> > MoleculeList;
00317 MoleculeList list;
00318
00319 static void sort(MoleculeList::iterator start,
00320 MoleculeList::iterator end);
00321 static void sort(molecule **start, molecule **end);
00322 };
00323
00324 extern t_mole_global mole_global;
00325
00326 class t_mole_local
00327 {
00328 public:
00329 void set_location( long nelem, long ion, double *dense );
00330 void set_isotope_abundances( void );
00331 double sink_rate_tot(const char chSpecies[]) const;
00332 double sink_rate_tot(const molecule* const sp) const;
00333 double sink_rate(const molecule* const sp, const mole_reaction& rate) const;
00334 double sink_rate(const molecule* const sp, const char buf[]) const;
00335 double source_rate_tot(const char chSpecies[]) const;
00336 double source_rate_tot(const molecule* const sp) const;
00339 double dissoc_rate(const char chSpecies[]) const;
00340 double chem_heat(void) const;
00341 double findrk(const char buf[]) const;
00342 double findrate(const char buf[]) const;
00343
00344 double grain_area, grain_density, grain_saturation;
00345
00347 double elec;
00348
00351 double **source , **sink;
00352
00353 realnum ***xMoleChTrRate;
00354
00355 valarray<class molezone> species;
00356
00357 vector<double> reaction_rks;
00358 vector<double> old_reaction_rks;
00359 long old_zone;
00360 };
00361
00362 extern t_mole_local mole;
00363
00364 class molezone {
00365 public:
00366 molezone()
00367 {
00368 init();
00369 }
00370 void init (void)
00371 {
00372 location = NULL;
00373 levels = NULL;
00374 lines = NULL;
00375 zero();
00376 }
00377 void zero (void)
00378 {
00379 src = 0.;
00380 snk = 0.;
00381 den = 0.;
00382 column = 0.;
00383 nAtomLim = -1;
00384 xFracLim = 0.;
00385 column_old = 0.;
00386 }
00387 double *location;
00390 double src, snk;
00391
00392 qList* levels;
00393 TransitionList* lines;
00394
00395
00396 double den;
00397 realnum column;
00398 int nAtomLim;
00399 realnum xFracLim;
00401
00402 realnum column_old;
00403 };
00404
00405 extern molezone *null_molezone;
00406
00407 extern molezone *findspecieslocal(const char buf[]);
00408
00409 extern void mole_punch(FILE *punit, const char speciesname[], const char args[], bool lgHeader, bool lgData, double depth);
00410
00411 extern void total_molecule_elems(realnum total[LIMELM]);
00412 extern void total_molecule_deut(realnum &total);
00413
00414 extern realnum total_molecules(void);
00415
00416 extern realnum total_molecules_gasphase(void);
00417
00418 extern void mole_make_list(void);
00419 extern void mole_make_groups(void);
00420
00421 void mole_cmp_num_in_out_reactions(void);
00422
00423 bool lgDifferByExcitation( const molecule &mol1, const molecule &mol2 );
00424
00425 extern void mole_update_species_cache(void);
00426
00427 void mole_update_sources(void);
00428
00429 void mole_rk_bigchange(void);
00430
00431 void create_isotopologues_one(
00432 ChemAtomList& atoms,
00433 vector< int >& numAtoms,
00434 string atom_old,
00435 string atom_new,
00436 string embellishments,
00437 vector<string>& newLabels );
00438
00439 bool parse_species_label( const char label[], ChemAtomList &atomsLeftToRight, vector<int> &numAtoms, string &embellishments );
00440 bool parse_species_label( const char mylab[], ChemAtomList &atomsLeftToRight, vector<int> &numAtoms, string &embellishments,
00441 bool &lgExcit, int &charge, bool &lgGas_Phase );
00442
00443 #endif
00444