45 realnum mass_amu,
double frac );
59 map <string,count_ptr<molecule> >
spectab;
60 map <string,count_ptr<mole_reaction> >
reactab;
61 map <string,count_ptr<chem_element> >
elemtab;
62 map <string,count_ptr<chem_atom> >
atomtab;
63 map <string,count_ptr<mole_reaction> >
functab;
78 class MoleCmp :
public binary_function<const count_ptr<molecule>,
79 const count_ptr<molecule>,bool>
85 return mol1->compare(*mol2) < 0;
89 return mol1->compare(*mol2) < 0;
95 t_mole_global::MoleculeList::iterator end)
97 std::sort(
start,end,MoleCmp());
101 std::sort(
start,end,MoleCmp());
112 static const char chFile[] =
"isotope_fractions.dat";
116 while(
read_whole_line( chLine, (
int)
sizeof(chLine), ioDATA ) != NULL )
118 if( chLine[0] ==
'#' )
128 if( (
unsigned)Z <= lgResolveNelem.size() && lgResolveNelem[Z-1] )
170 if(
co.C12_C13_isotope_ratio_parsed >= 0. )
172 atom12C->frac =
co.C12_C13_isotope_ratio_parsed /
173 (
co.C12_C13_isotope_ratio_parsed + 1.);
174 atom13C->frac = 1. / (
co.C12_C13_isotope_ratio_parsed + 1.);
175 co.C12_C13_isotope_ratio =
co.C12_C13_isotope_ratio_parsed;
179 co.C12_C13_isotope_ratio = atom12C->frac /
SDIV( atom13C->frac );
194 for(
long nelem=0; nelem<
LIMELM; ++nelem )
250 if (
hmi.lgLeiden_Keep_ipMH2s)
265 long int nelem = (*atom)->el->Z-1;
267 for(
long ion=0; ion<=nelem+1; ion++ )
274 sprintf( temp,
"+" );
276 sprintf( temp,
"+%ld", ion );
277 sprintf( str,
"%s%s", (*atom)->label().c_str(), temp );
350 atom->ipMl[sp->
charge] = i;
367 while( getline( ioDATA,line ) )
373 istringstream iss( line );
375 double formation_enthalpy;
377 iss >> formation_enthalpy;
389 vector< int >& numAtoms,
392 string embellishments,
393 vector<string>& newLabels )
399 for(
unsigned position = 0; position <
atoms.size(); ++position )
402 if(
atoms[position]->label() == atom_old )
404 for(
unsigned i=0; i<position; ++i )
406 str <<
atoms[i]->label();
411 if( numAtoms[position] > 1 )
414 if( numAtoms[position] > 2 )
415 str << numAtoms[position]-1;
419 if( position+1 ==
atoms.size() )
422 for(
unsigned i=position+1; i<
atoms.size(); ++i )
427 if( atom_new ==
atoms[i]->label() )
429 str <<
atoms[i]->label();
430 ASSERT( numAtoms[i] + 1 > 1 );
431 str << numAtoms[i] + 1;
436 str <<
atoms[i]->label();
437 if( numAtoms[i] > 1 )
443 str <<
atoms[i]->label();
444 if( numAtoms[i] > 1 )
450 str << embellishments;
452 newLabels.push_back( str.str() );
468 int len = strlen(label)+1;
470 char *mylab = mylab_v.
data();
471 strncpy(mylab,label,len);
472 s = strchr(mylab,
' ');
488 realnum mass_amu,
double frac )
493 ASSERT( massNumberA < 3 *
LIMELM && ( massNumberA > 0 || massNumberA == -1 ) );
494 ASSERT( mass_amu < 3. * LIMELM && mass_amu > 0. );
495 ASSERT( frac <= 1. + FLT_EPSILON && frac >= 0. );
498 isotope->A = massNumberA;
499 isotope->mass_amu = mass_amu;
500 isotope->frac = frac;
501 isotope->ipMl.resize(el->Z+1);
503 for (
long int ion = 0; ion < el->Z+1; ion++)
504 isotope->ipMl[ion] = -1;
509 if( isotope->lgMeanAbundance() )
512 el->isotopes[massNumberA] = isotope;
520 if( el->isotopes.size()==0 )
521 return dense.AtomicWeight[el->Z-1];
523 realnum averageMass = 0., fracsum = 0.;
524 for(
isotopes_i it = el->isotopes.begin(); it != el->isotopes.end(); ++it )
526 fracsum += it->second->frac;
527 averageMass += it->second->mass_amu * it->second->frac;
549 bool lgCreateIsotopologues )
555 vector< int > numAtoms;
556 string embellishments;
560 mol->parentLabel.clear();
561 mol->isEnabled =
true;
563 mol->lgExcit =
false;
564 mol->mole_mass = 0.0;
566 mol->lgGas_Phase =
true;
567 mol->form_enthalpy = form_enthalpy;
571 int len = strlen(label)+1;
573 char *mylab = mylab_v.
data();
574 strncpy(mylab,label,len);
575 s = strchr(mylab,
' ');
582 if(
parse_species_label( mylab, atomsLeftToRight, numAtoms, embellishments, mol->lgExcit, mol->charge, mol->lgGas_Phase ) ==
false )
585 for(
unsigned i = 0; i < atomsLeftToRight.size(); ++i )
587 mol->nAtom[ atomsLeftToRight[i] ] += numAtoms[i];
596 if ( (mol->n_nuclei() > 1 || (mol->isMonatomic() && mol->charge==-1) ) &&
mole_global.lgNoMole)
598 if(
trace.lgTraceMole )
599 fprintf(
ioQQQ,
"No species %s as molecules off\n",label);
603 if (mol->n_nuclei() > 1 &&
mole_global.lgNoHeavyMole)
605 for(
nAtoms_ri it=mol->nAtom.rbegin(); it != mol->nAtom.rend(); --it )
610 if(
trace.lgTraceMole )
611 fprintf(
ioQQQ,
"No species %s as heavy molecules off\n",label);
621 fprintf(
ioQQQ,
"Warning: duplicate species %s - using first one\n",
622 mol->label.c_str() );
629 for(
nAtoms_i it=mol->nAtom.begin(); it != mol->nAtom.end(); ++it )
638 if( lgCreateIsotopologues && type==
MOLECULE && mol->label.compare(
"CO")==0 )
645 if( lgCreateIsotopologues && type==
MOLECULE && !mol->isMonatomic() )
647 for(
nAtoms_i it1 = mol->nAtom.begin(); it1 != mol->nAtom.end(); ++it1 )
650 it2 != it1->first->el->isotopes.end(); ++it2 )
653 if( it1->first->el->Z-1 ==
ipHYDROGEN && it2->second->A != 2 )
655 if( !
mole_global.lgTreatIsotopes[it1->first->el->Z-1] )
657 if( it2->second->lgMeanAbundance() )
659 vector<string> isoLabs;
660 create_isotopologues_one( atomsLeftToRight, numAtoms, it1->first->label(), it2->second->label(), embellishments, isoLabs );
665 for( vector<string>::iterator newLabel = isoLabs.begin(); newLabel != isoLabs.end(); ++newLabel )
669 if( sp!=NULL && it2->second->A != 2 )
680 bool lgExcit, lgGas_Phase;
682 bool lgOK =
parse_species_label( label, atomsLeftToRight, numAtoms, embellishments, lgExcit, charge, lgGas_Phase );
686 bool &lgExcit,
int &charge,
bool &lgGas_Phase )
688 long int i, n, ipAtom;
697 s = strpbrk(mylab,
"*");
706 s = strpbrk(mylab,
"+-");
717 embellishments = s + embellishments;
721 s = strstr(mylab,
"grn");
725 embellishments = s + embellishments;
736 while (mylab[i] !=
'\0' && mylab[i] !=
' ' && mylab[i] !=
'*')
743 thisAtom[ipAtom++] = mylab[i++];
744 ASSERT( isdigit(mylab[i]) );
745 thisAtom[ipAtom++] = mylab[i++];
746 if(isdigit(mylab[i]))
748 thisAtom[ipAtom++] = mylab[i++];
752 thisAtom[ipAtom++] = mylab[i++];
753 if(islower(mylab[i]))
755 thisAtom[ipAtom++] = mylab[i++];
757 thisAtom[ipAtom] =
'\0';
763 fprintf(stderr,
"Did not recognize atom at %s in \"%s \"[%ld]\n",thisAtom,mylab,i);
766 if(!
dense.lgElmtOn[atom->el->Z-1])
768 if(
trace.lgTraceMole )
769 fprintf(
ioQQQ,
"No species %s as element %s off\n",mylab,atom->el->label.c_str() );
773 if(isdigit(mylab[i]))
777 n = 10*n+(
long int)(mylab[i]-
'0');
786 atomsLeftToRight.push_back( atom );
787 numAtoms.push_back( n );
820 for (
const char *pb = buf; *pb && *pb !=
' '; ++pb)
828 return &(*p->second);
839 for (
const char *pb = buf; *pb && *pb !=
' '; ++pb)
847 return &
mole.species[ p->second->index ];
869 double den_times_area, den_grains, adsorbed_density;
873 enum { DEBUG_MOLE =
false };
886 den_times_area = den_grains = 0.0;
887 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
890 den_times_area +=
gv.bin[nd]->IntArea/4.*
gv.bin[nd]->cnv_H_pCM3;
891 den_grains +=
gv.bin[nd]->cnv_GR_pCM3;
894 adsorbed_density = 0.0;
898 adsorbed_density +=
mole.species[i].den;
901 mole.grain_area = den_times_area;
902 mole.grain_density = den_grains;
904 double mole_cs = 1e-15;
905 if (4*den_times_area <= mole_cs*adsorbed_density)
906 mole.grain_saturation = 1.0;
908 mole.grain_saturation = ((mole_cs*adsorbed_density)/(4.*den_times_area));
912 mole.grain_area = 0.0;
913 mole.grain_density = 0.0;
914 mole.grain_saturation = 1.0;
918 fprintf(
ioQQQ,
"Dust: %f %f %f\n",
919 mole.grain_area,
mole.grain_density,
mole.grain_saturation);
923 if(
mole.species[i].location != NULL)
926 mole.species[i].den = *(
mole.species[i].location);
932 mole.set_isotope_abundances();
946 dense.updateXMolecules();
953 if (
mole.species[i].location == NULL &&
mole_global.list[i]->parentLabel.empty())
970 long nelem =
mole_global.list[i]->nAtom.begin()->first->el->Z-1;
971 realnum frac = (new_pop-old_pop)/
SDIV(new_pop+old_pop+1e-8*
972 dense.gas_phase[nelem]);
980 *(
mole.species[i].location) = new_pop;
986 return ncpt>0 ? sqrt(delta/ncpt) : 0.f;
991 DEBUG_ENTRY(
"t_mole_local::set_isotope_abundances()");
997 for(
isotopes_i it = (*atom)->el->isotopes.begin(); it != (*atom)->el->isotopes.end(); ++it )
1000 if( it->second->lgMeanAbundance() )
1004 for(
unsigned long ion=0; ion<it->second->ipMl.size(); ++ion )
1006 if( it->second->ipMl[ion] != -1 &&
1007 (
species[ it->second->ipMl[ion] ].location == NULL ) && (*atom)->ipMl[ion] != -1 )
1009 species[ it->second->ipMl[ion] ].den =
species[ (*atom)->ipMl[ion] ].den * it->second->frac;
1023 ASSERT( ion < nelem + 2 );
1026 if( mole_index == -1 )
1029 species[mole_index].location = density;
1040 if( !
deut.lgElmtOn )
1045 if (
mole.species[i].location == NULL &&
mole_global.list[i]->parentLabel.empty() )
1047 for( molecule::nAtomsMap::iterator atom=
mole_global.list[i]->nAtom.begin();
1050 long int nelem = atom->first->el->Z-1;
1051 if( nelem==0 && atom->first->A==2 )
1053 total +=
mole.species[i].den*atom->second;
1074 if (
mole.species[i].location == NULL &&
mole_global.list[i]->parentLabel.empty() )
1076 for( molecule::nAtomsMap::iterator atom=
mole_global.list[i]->nAtom.begin();
1079 ASSERT( atom->second > 0 );
1080 long int nelem = atom->first->el->Z-1;
1081 if( atom->first->lgMeanAbundance() )
1082 total[nelem] += (
realnum)
mole.species[i].den*atom->second;
1100 for( molecule::nAtomsMap::iterator atom=
mole_global.list[i]->nAtom.begin();
1103 long int nelem = atom->first->el->Z-1;
1104 total[nelem] += (
realnum)
mole.species[i].den*atom->second;
1119 if (
mole.species[i].location == NULL &&
mole_global.list[i]->parentLabel.empty())
1184 ASSERT( it->second > 0 );
const int INPUT_LINE_LENGTH
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
bool fp_equal(sys_float x, sys_float y, int n=3)
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
sys_float SDIV(sys_float x)
#define DEBUG_ENTRY(funcname)
element_type * data() const
bool isMonatomic(void) const
iterator begin(size_type i1)
static void sort(MoleculeList::iterator start, MoleculeList::iterator end)
vector< bool > lgTreatIsotopes
void set_location(long nelem, long ion, double *dense)
void set_isotope_abundances(void)
valarray< class molezone > species
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
const ios_base::openmode mode_r
bool lgElemsConserved(void)
void SetGasPhaseDeuterium(const realnum &Hdensity)
void InitDeuteriumIonization()
void SetDeuteriumFractionation(const realnum &frac)
t_elementnames elementnames
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
t_mole_global mole_global
molecule::nAtomsMap::reverse_iterator nAtoms_ri
chem_element * null_element
t_mole_global mole_global
ChemAtomList unresolved_atom_list
molecule::nAtomsMap::iterator nAtoms_i
molecule * findspecies(const char buf[])
map< int, count_ptr< chem_atom > >::iterator isotopes_i
vector< count_ptr< chem_atom > > ChemAtomList
map< string, count_ptr< chem_atom > >::iterator chem_atom_i
map< string, count_ptr< molecule > >::iterator molecule_i
STATIC realnum MeanMassOfElement(const count_ptr< chem_element > &el)
STATIC bool isactive(const molecule &mol)
STATIC void newisotope(const count_ptr< chem_element > &el, int massNumberA, realnum mass_amu, double frac)
realnum total_molecules_gasphase(void)
STATIC void ReadIsotopeFractions(const vector< bool > &lgResolveNelem)
void mole_make_groups(void)
realnum total_molecules(void)
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[])
STATIC void newelement(const char label[], int ipion)
vector< count_ptr< chem_element > > element_list
STATIC bool ispassive(const molecule &mol)
STATIC void read_species_file(string filename, bool lgCreateIsotopologues=true)
realnum mole_return_cached_species(const GroupMap &)
void mole_update_species_cache(void)
bool lgDifferByExcitation(const molecule &mol1, const molecule &mol2)
STATIC count_ptr< chem_atom > findatom(const char buf[])
void total_network_elems(double total[LIMELM])
bool parse_species_label(const char label[], ChemAtomList &atomsLeftToRight, vector< int > &numAtoms, string &embellishments)
molecule * findspecies(const char buf[])
STATIC molecule * newspecies(const char label[], enum spectype type, enum mole_state state, realnum form_enthalpy)
void total_molecule_elems(realnum total[LIMELM])
void total_molecule_deut(realnum &total_f)
map< string, count_ptr< chem_atom > > atomtab
map< string, count_ptr< mole_reaction > > functab
map< string, count_ptr< mole_reaction > > reactab
map< string, count_ptr< molecule > > spectab
map< string, count_ptr< chem_element > > elemtab
STATIC void start(long, realnum[], realnum[], long, long[], realnum *, int *)
UNUSED const double ELECTRON_MASS
UNUSED const double KJMOL1CM
UNUSED const double ATOMIC_MASS_UNIT