41 map <string,count_ptr<molecule> >
spectab;
42 map <string,count_ptr<mole_reaction> >
reactab;
43 map <string,count_ptr<chem_element> >
elemtab;
44 map <string,count_ptr<mole_reaction> >
functab;
61 class MoleCmp :
public binary_function<const count_ptr<molecule>,
62 const count_ptr<molecule>,bool>
68 return mol1->
compare(*mol2) < 0;
72 return mol1->
compare(*mol2) < 0;
78 t_mole_global::MoleculeList::iterator end)
80 std::sort(start,end,MoleCmp());
84 std::sort(start,end,MoleCmp());
91 for(
int ielm = 0; ielm <
LIMELM; ielm++ )
100 if( (
unsigned)Z <= lgResolveNelem.size() && lgResolveNelem[Z] )
121 null_molezone->
den = 0.;
144 for(
long nelem=0; nelem<
LIMELM; ++nelem )
164 null_mole->
index = -1;
210 for( ChemNuclideList::iterator atom =
nuclide_list.begin();
213 if ( !(*atom)->lgHasLinkedIon() )
215 long int nelem = (*atom)->el->Z-1;
217 for(
long ion=0; ion<=nelem+1; ion++ )
224 sprintf( temp,
"+" );
226 sprintf( temp,
"+%ld", ion );
227 sprintf( str,
"%s%s", (*atom)->label().c_str(), temp );
231 fixit(
"populate these in a local update");
300 atom->ipMl[sp->
charge] = i;
317 while( getline( ioDATA,line ) )
323 istringstream iss( line );
325 double formation_enthalpy;
327 iss >> formation_enthalpy;
339 vector< int >& numNuclides,
342 string embellishments,
343 vector<string>& newLabels )
347 fixit(
"make sure atom_new and atom_old are isotopes");
348 fixit(
"make sure atom_new is not already present");
351 for(
unsigned position = 0; position < atoms.size(); ++position )
355 if( !newLabel.empty() )
356 newLabels.push_back( newLabel );
365 vector< int >& numNuclides,
368 string embellishments,
371 DEBUG_ENTRY(
"create_isotopologues_one_position()" );
374 if( atoms[position]->label() == atom_old )
376 for(
unsigned i=0; i<position; ++i )
378 str << atoms[i]->label();
379 if( numNuclides[i]>1 )
380 str << numNuclides[i];
383 if( numNuclides[position] > 1 )
386 if( numNuclides[position] > 2 )
387 str << numNuclides[position]-1;
390 if( position+1 == atoms.size() )
393 for(
unsigned i=position+1; i<atoms.size(); ++i )
398 if( atom_new == atoms[i]->label() )
400 str << atoms[i]->label();
401 ASSERT( numNuclides[i] + 1 > 1 );
402 str << numNuclides[i] + 1;
407 str << atoms[i]->label();
408 if( numNuclides[i] > 1 )
409 str << numNuclides[i];
414 str << atoms[i]->label();
415 if( numNuclides[i] > 1 )
416 str << numNuclides[i];
421 str << embellishments;
423 newLabel = str.str();
438 int len = strlen(label)+1;
440 char *mylab = mylab_v.
data();
441 strncpy(mylab,label,len);
442 s = strchr(mylab,
' ');
463 ASSERT( massNumberA < 3 *
LIMELM && ( massNumberA > 0 || massNumberA == -1 ) );
464 ASSERT( mass_amu < 3. * LIMELM && mass_amu > 0. );
465 ASSERT( frac <= 1. + FLT_EPSILON && frac >= 0. );
468 isotope->
A = massNumberA;
471 isotope->
ipMl.resize(el->
Z+1);
473 for (
long int ion = 0; ion < el->
Z+1; ion++)
474 isotope->
ipMl[ion] = -1;
480 el->
isotopes[massNumberA] = isotope;
491 realnum averageMass = 0., fracsum = 0.;
494 fracsum += it->second->frac;
495 averageMass += it->second->mass_amu * it->second->frac;
508 return newspecies( label, type, state, form_enthalpy,
true);
517 bool lgCreateIsotopologues )
523 vector< int > numNuclides;
524 string embellishments;
540 int len = strlen(label)+1;
542 char *mylab = mylab_v.
data();
543 strncpy(mylab,label,len);
544 s = strchr(mylab,
' ');
554 for(
unsigned i = 0; i < nuclidesLeftToRight.size(); ++i )
556 mol->
nNuclide[ nuclidesLeftToRight[i] ] += numNuclides[i];
557 mol->
mole_mass += numNuclides[i] * nuclidesLeftToRight[i]->mass_amu * (
realnum)ATOMIC_MASS_UNIT;
568 fprintf(
ioQQQ,
"No species %s as molecules off\n",label);
580 fprintf(
ioQQQ,
"No species %s as heavy molecules off\n",label);
596 fprintf(
ioQQQ,
"Warning: duplicate species %s - using first one\n",
597 mol->
label.c_str() );
613 if( lgCreateIsotopologues && type==
MOLECULE && mol->
label.compare(
"CO")==0 )
625 it2 != it1->first->el->isotopes.end(); ++it2 )
628 if( it1->first->el->Z-1 ==
ipHYDROGEN && it2->second->A != 2 )
632 if( it2->second->lgMeanAbundance() )
634 vector<string> isoLabs;
635 create_isotopologues( nuclidesLeftToRight, numNuclides, it1->first->label(), it2->second->label(), embellishments, isoLabs );
640 for( vector<string>::iterator newLabel = isoLabs.begin(); newLabel != isoLabs.end(); ++newLabel )
644 if( sp!=NULL && it2->second->A != 2 )
655 bool lgExcit, lgGas_Phase;
657 bool lgOK =
parse_species_label( label, nuclidesLeftToRight, numNuclides, embellishments, lgExcit, charge, lgGas_Phase );
661 bool &lgExcit,
int &charge,
bool &lgGas_Phase )
663 long int i, n, ipAtom;
667 int iend = strlen(label);
670 s = strpbrk(label,
"*");
679 s = strpbrk(label,
"+-");
690 embellishments = s + embellishments;
694 s = strstr(label,
"grn");
698 embellishments = s + embellishments;
710 while (i < iend && label[i] !=
' ' && label[i] !=
'*')
717 thisAtom[ipAtom++] = label[i++];
718 ASSERT( isdigit(label[i]) );
719 thisAtom[ipAtom++] = label[i++];
720 if(isdigit(label[i]))
722 thisAtom[ipAtom++] = label[i++];
726 thisAtom[ipAtom++] = label[i++];
727 if( i<iend && islower(label[i]))
729 thisAtom[ipAtom++] = label[i++];
731 thisAtom[ipAtom] =
'\0';
737 fprintf(stderr,
"Did not recognize atom at %s in \"%s \"[%ld]\n",thisAtom,label,i);
747 if(i < iend && isdigit(label[i]))
751 n = 10*n+(
long int)(label[i]-
'0');
753 }
while (i < iend && isdigit(label[i]));
759 nuclidesLeftToRight.push_back( atom );
760 numNuclides.push_back( n );
793 for (
const char *pb = buf; *pb && *pb !=
' '; ++pb)
801 return &(*p->second);
812 for (
const char *pb = buf; *pb && *pb !=
' '; ++pb)
820 return &(*p->second);
823 fprintf(
ioQQQ,
" PROBLEM specified species, '%s', is not valid or disabled.\n", s.c_str() );
824 fprintf(
ioQQQ,
" The SAVE SPECIES LABELS command will generate a list of species. Species names are case sensitive.\n");
867 double den_times_area, den_grains, adsorbed_density;
871 enum { DEBUG_MOLE =
false };
884 den_times_area = den_grains = 0.0;
885 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
888 den_times_area +=
gv.
bin[nd]->IntArea/4.*
gv.
bin[nd]->cnv_H_pCM3;
889 den_grains +=
gv.
bin[nd]->cnv_GR_pCM3;
892 adsorbed_density = 0.0;
902 double mole_cs = 1e-15;
903 if (4*den_times_area <= mole_cs*adsorbed_density)
984 return ncpt>0 ? sqrt(delta/ncpt) : 0.f;
989 DEBUG_ENTRY(
"t_mole_local::set_isotope_abundances()" );
994 if ( !(*atom)->lgMeanAbundance() )
998 for(
isotopes_i it = (*atom)->el->isotopes.begin(); it != (*atom)->el->isotopes.end(); ++it )
1001 if( it->second->lgMeanAbundance() )
1005 for(
unsigned long ion=0; ion<it->second->ipMl.size(); ++ion )
1007 if( it->second->ipMl[ion] != -1 &&
1008 (
species[ it->second->ipMl[ion] ].location == NULL ) && (*atom)->ipMl[ion] != -1 )
1010 species[ it->second->ipMl[ion] ].den =
species[ (*atom)->ipMl[ion] ].den * it->second->frac;
1021 DEBUG_ENTRY(
"t_mole_local::set_ion_locations()" );
1023 for( ChemNuclideList::iterator atom =
nuclide_list.begin();
1026 if ( !(*atom)->lgHasLinkedIon() )
1028 long nelem = (*atom)->el->Z-1;
1029 for(
long ion=0; ion<nelem+2; ++ion )
1031 long mole_index = (*atom)->ipMl[ion];
1033 if( mole_index != -1 )
1062 for( molecule::nNucsMap::iterator atom=
mole_global.
list[i]->nNuclide.begin();
1065 long int nelem = atom->first->el->Z-1;
1066 if( nelem==0 && atom->first->A==2 )
1091 for( molecule::nNucsMap::iterator atom=
mole_global.
list[i]->nNuclide.begin();
1094 ASSERT( atom->second > 0 );
1095 long int nelem = atom->first->el->Z-1;
1096 if( atom->first->lgMeanAbundance() )
1176 ASSERT( it->second > 0 );
map< string, size_t >::iterator chem_nuclide_i
t_mole_global mole_global
static void sort(MoleculeList::iterator start, MoleculeList::iterator end)
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
const double ENERGY_H2_STAR
molecule::nNucsMap::iterator nNucs_i
realnum getMass(int Anum)
vector< bool > lgTreatIsotopes
realnum total_molecules(void)
map< int, count_ptr< chem_nuclide > >::iterator isotopes_i
void total_molecule_deut(realnum &total)
map< string, count_ptr< molecule > >::iterator molecule_i
bool lgLeiden_Keep_ipMH2s
void mole_make_list(void)
map< string, count_ptr< mole_reaction > > reactab
STATIC molecule * newspecies(const char label[], enum spectype type, enum mole_state state, realnum form_enthalpy)
molezone * ionMole[LIMELM][LIMELM+1]
bool speciesOff(const string &label)
vector< count_ptr< chem_element > > ChemElementList
map< string, size_t > nuclidetab
molezone * findspecieslocal(const char buf[])
void InitDeuteriumIonization()
ChemNuclideList nuclide_list
bool parse_species_label(const char label[], ChemNuclideList &atomsLeftToRight, vector< int > &numAtoms, string &embellishments)
molezone * findspecieslocal_validate(const char buf[])
t_elementnames elementnames
ChemElementList element_list
double xIonDense[LIMELM][LIMELM+1]
element_type * data() const
bool fp_equal(sys_float x, sys_float y, int n=3)
STATIC realnum MeanMassOfElement(const count_ptr< chem_element > &el)
molecule::nNucsMap::reverse_iterator nNucs_ri
const ios_base::openmode mode_r
map< int, count_ptr< chem_nuclide > > isotopes
STATIC bool isactive(const molecule &mol)
map< string, count_ptr< chem_element > > elemtab
STATIC void read_species_file(string filename, bool lgCreateIsotopologues=true)
chem_nuclide * null_nuclide
molecule * findspecies(const char buf[])
void create_isotopologues_one_position(unsigned position, ChemNuclideList &atoms, vector< int > &numAtoms, string atom_old, string atom_new, string embellishments, string &newLabel)
valarray< class molezone > species
count_ptr< chem_nuclide > findnuclide(const char buf[])
STATIC void SetIsotopeFractions(const vector< bool > &lgResolveNelem)
realnum total_molecules_gasphase(void)
map< string, count_ptr< mole_reaction > > functab
STATIC void newelement(const char label[], int ipion)
realnum mole_return_cached_species(const GroupMap &MoleMap)
diatomics h2("h2", 4100.,&hmi.H2_total, Yan_H2_CS)
bool exists(const molecule *m)
void SetDeuteriumFractionation(const realnum &frac)
void SetGasPhaseDeuterium(const realnum &Hdensity)
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
int compare(const molecule &mol2) const
realnum gas_phase[LIMELM]
vector< count_ptr< chem_nuclide > > ChemNuclideList
realnum AtomicWeight[LIMELM]
void set_isotope_abundances(void)
STATIC void start(long, realnum[], realnum[], long, long[], realnum *, int *)
bool isMonatomic(void) const
#define DEBUG_ENTRY(funcname)
STATIC void newisotope(const count_ptr< chem_element > &el, int massNumberA, realnum mass_amu, double frac)
bool lgGrain_mole_deplete
bool lgDifferByExcitation(const molecule &mol1, const molecule &mol2)
void mole_update_species_cache(void)
STATIC bool ispassive(const molecule &mol)
int fprintf(const Output &stream, const char *format,...)
vector< molecule * > groupspecies
sys_float SDIV(sys_float x)
bool lgElemsConserved(void)
void mole_make_groups(void)
map< string, count_ptr< molecule > > spectab
chem_element * null_element
void total_molecule_elems(realnum total[LIMELM])
molecule * findspecies_validate(const char buf[])
void create_isotopologues(ChemNuclideList &atoms, vector< int > &numAtoms, string atom_old, string atom_new, string embellishments, vector< string > &newLabels)