00001
00002
00003
00004
00005 #include "cdstd.h"
00006 #include <cctype>
00007 #include <string.h>
00008 #include <algorithm>
00009 #include <stdlib.h>
00010 #include "cddefines.h"
00011 #include "colden.h"
00012 #include "conv.h"
00013 #include "deuterium.h"
00014 #include "doppvel.h"
00015 #include "elementnames.h"
00016 #include "h2.h"
00017 #include "iso.h"
00018 #include "phycon.h"
00019 #include "physconst.h"
00020 #include "mole.h"
00021 #include "mole_priv.h"
00022 #include "hmi.h"
00023 #include "radius.h"
00024 #include "rfield.h"
00025 #include "rt.h"
00026 #include "secondaries.h"
00027 #include "dense.h"
00028 #include "ionbal.h"
00029 #include "grainvar.h"
00030 #include "timesc.h"
00031 #include "taulines.h"
00032 #include "trace.h"
00033
00034
00035
00036
00037
00038
00039 enum spectype {MOLECULE, OTHER};
00040
00041 STATIC void read_species_file( string filename, bool lgCreateIsotopologues = true );
00042 STATIC void newelement(const char label[], int ipion);
00043 STATIC void newisotope( const count_ptr<chem_element> &el, int massNumberA,
00044 realnum mass_amu, double frac );
00045 STATIC realnum MeanMassOfElement( const count_ptr<chem_element> &el );
00046
00047 STATIC molecule *newspecies(const char label[], enum spectype type,
00048 enum mole_state state, realnum form_enthalpy);
00049 STATIC molecule *newspecies(const char label[], enum spectype type,
00050 enum mole_state state, realnum form_enthalpy, bool lgCreateIsotopologues);
00051 STATIC count_ptr<chem_atom> findatom(const char buf[]);
00052 STATIC bool isactive(const molecule &mol);
00053 STATIC bool ispassive(const molecule &mol);
00054 STATIC void ReadIsotopeFractions( const vector<bool>& lgResolveNelem );
00055
00056
00057 namespace mole_priv {
00058 map <string,count_ptr<molecule> > spectab;
00059 map <string,count_ptr<mole_reaction> > reactab;
00060 map <string,count_ptr<chem_element> > elemtab;
00061 map <string,count_ptr<chem_atom> > atomtab;
00062 map <string,count_ptr<mole_reaction> > functab;
00063 };
00064 molecule *null_mole;
00065 molezone *null_molezone;
00066 chem_element *null_element;
00067 chem_atom *null_atom;
00068 vector< count_ptr <chem_element> > element_list;
00069 ChemAtomList unresolved_atom_list;
00070 ChemAtomList atom_list;
00071 molecule **groupspecies;
00072
00073 #include <functional>
00074
00075 namespace
00076 {
00077 class MoleCmp : public binary_function<const count_ptr<molecule>,
00078 const count_ptr<molecule>,bool>
00079 {
00080 public:
00081 bool operator()(const count_ptr<molecule> &mol1,
00082 const count_ptr<molecule> &mol2) const
00083 {
00084 return mol1->compare(*mol2) < 0;
00085 }
00086 bool operator()(const molecule *mol1, const molecule *mol2) const
00087 {
00088 return mol1->compare(*mol2) < 0;
00089 }
00090 };
00091 }
00092
00093 void t_mole_global::sort(t_mole_global::MoleculeList::iterator start,
00094 t_mole_global::MoleculeList::iterator end)
00095 {
00096 std::sort(start,end,MoleCmp());
00097 }
00098 void t_mole_global::sort(molecule **start, molecule **end)
00099 {
00100 std::sort(start,end,MoleCmp());
00101 }
00102
00103 STATIC void ReadIsotopeFractions( const vector<bool>& lgResolveNelem )
00104 {
00105 DEBUG_ENTRY( "ReadIsotopeFractions()" );
00106
00107 char chLine[INPUT_LINE_LENGTH];
00108 long i;
00109 bool lgEOL;
00110
00111 static const char chFile[] = "isotope_fractions.dat";
00112 FILE *ioDATA = open_data( chFile, "r" );
00113 ASSERT( ioDATA != NULL );
00114
00115 while( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) != NULL )
00116 {
00117 if( chLine[0] == '#' )
00118 continue;
00119
00120 i = 1;
00121 long Z = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00122 long A = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00123
00124 double frac = 0.01 * (double)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00125
00126 fixit();
00127 if( (unsigned)Z <= lgResolveNelem.size() && lgResolveNelem[Z-1] )
00128 newisotope( element_list[Z-1], A, (realnum)A, frac );
00129
00130 else if( Z-1==ipCARBON )
00131 newisotope( element_list[Z-1], A, (realnum)A, frac );
00132 }
00133
00134 fclose( ioDATA );
00135
00136 return;
00137 }
00138
00139 void t_mole_global::make_species(void)
00140 {
00141 DEBUG_ENTRY( "mole_global::make_species()" );
00142
00143 long int i;
00144 molecule *sp;
00145
00146 null_element = new chem_element(0,"") ;
00147 null_atom = new( chem_atom );
00148 null_molezone = new( molezone );
00149 null_molezone->den = 0.;
00150
00151
00152 for (i=0;i<LIMELM;i++)
00153 {
00154 newelement(elementnames.chElementSym[i], i);
00155 }
00156
00157
00158 if( deut.lgElmtOn )
00159 lgTreatIsotopes[ipHYDROGEN] = true;
00160
00161
00162 ReadIsotopeFractions( lgTreatIsotopes );
00163
00164 if( lgTreatIsotopes[ipHYDROGEN] )
00165 {
00166 SetDeuteriumFractionation( element_list[ipHYDROGEN]->isotopes[2]->frac );
00167 SetGasPhaseDeuterium( dense.gas_phase[ipHYDROGEN] );
00168 InitDeuteriumIonization();
00169 }
00170
00171 for( long nelem=0; nelem<LIMELM; ++nelem )
00172 {
00173 realnum mass_amu = MeanMassOfElement( element_list[nelem] );
00174
00175 newisotope( element_list[nelem], -1, mass_amu, 1.0 );
00176 }
00177
00178 ASSERT( (long) unresolved_atom_list.size() == LIMELM );
00179
00180
00181
00182
00183
00184
00185
00186 mole_global.num_total = mole_global.num_calc = 0;
00187
00188
00189
00190
00191
00192 null_mole = newspecies(" ",OTHER,MOLE_NULL, 0.);
00193 null_mole->index = -1;
00194
00195 read_species_file( "chem_species.dat" );
00196
00197 if(gv.lgDustOn() && mole_global.lgGrain_mole_deplete )
00198 {
00199
00200
00201
00202
00203
00204 if (0)
00205 {
00206 read_species_file( "chem_species_grn.dat" );
00207 }
00208 else
00209 {
00210 sp = newspecies("COgrn ",MOLECULE,MOLE_ACTIVE, -113.8f);
00211 sp = newspecies("H2Ogrn",MOLECULE,MOLE_ACTIVE, -238.9f);
00212 sp = newspecies("OHgrn ",MOLECULE,MOLE_ACTIVE, 38.4f);
00213
00214 }
00215 }
00216
00217
00218 sp = newspecies("e- ",OTHER,MOLE_PASSIVE, 0.0f);
00219 sp->charge = -1; sp->mole_mass = (realnum)ELECTRON_MASS;
00220 sp = newspecies("grn ",OTHER,MOLE_PASSIVE, 0.0f);
00221 sp = newspecies("PHOTON",OTHER,MOLE_PASSIVE, 0.0f);
00222 sp = newspecies("CRPHOT",OTHER,MOLE_PASSIVE, 0.0f);
00223 sp = newspecies("CRP ",OTHER,MOLE_PASSIVE, 0.0f);
00224
00225 if (!mole_global.lgLeidenHack)
00226 sp = newspecies("H- ",MOLECULE,MOLE_ACTIVE, 143.2f);
00227 if (hmi.lgLeiden_Keep_ipMH2s)
00228 sp = newspecies("H2* ",MOLECULE,MOLE_ACTIVE, 0.0f); fixit();
00229
00230 if( deut.lgElmtOn )
00231 {
00232 read_species_file( "chem_species_deuterium.dat", false );
00233 }
00234
00235
00236 for( ChemAtomList::iterator atom = unresolved_atom_list.begin();
00237 atom != unresolved_atom_list.end(); ++atom)
00238 {
00239 long int nelem = (*atom)->el->Z-1;
00240
00241 for( long ion=0; ion<=nelem+1; ion++ )
00242 {
00243 char str[CHARS_ISOTOPE_SYM+CHARS_ION_STAGE+1];
00244 char temp[CHARS_ION_STAGE+1];
00245 if( ion==0 )
00246 temp[0] = '\0';
00247 else if( ion==1 )
00248 sprintf( temp, "+" );
00249 else
00250 sprintf( temp, "+%ld", ion );
00251 sprintf( str, "%s%s", (*atom)->label().c_str(), temp );
00252 if (findspecies(str) == null_mole)
00253 {
00254 sp = newspecies(str,MOLECULE,MOLE_ACTIVE, 0.f);
00255 fixit();
00256 #if 0
00257 if( sp != NULL )
00258 {
00259 sp->levels = NULL;
00260 sp->numLevels = 0;
00261 }
00262 #endif
00263 }
00264 }
00265 }
00266
00267 return;
00268 }
00269
00270 void mole_make_list()
00271 {
00272 DEBUG_ENTRY( "mole_make_list()" );
00273
00274
00275 mole_global.list.resize(mole_global.num_total);
00276
00277
00278 long int i = 0;
00279 for (molecule_i p = mole_priv::spectab.begin(); p != mole_priv::spectab.end(); ++p)
00280 {
00281 if (isactive(*(p->second)))
00282 mole_global.list[i++] = p->second;
00283 }
00284 ASSERT (i == mole_global.num_calc);
00285
00286
00287 t_mole_global::sort(mole_global.list.begin(),mole_global.list.begin()+mole_global.num_calc);
00288
00289 for (molecule_i p = mole_priv::spectab.begin(); p != mole_priv::spectab.end(); ++p)
00290 {
00291 if (ispassive(*(p->second)))
00292 mole_global.list[i++] = p->second;
00293 }
00294 ASSERT (i == mole_global.num_total);
00295
00296
00297 for(i=0;i<mole_global.num_total;i++)
00298 {
00299 mole_global.list[i]->index = i;
00300 }
00301
00302 for(i=0;i<mole_global.num_total;i++)
00303 {
00304 if( !mole_global.list[i]->parentLabel.empty() )
00305 {
00306 long parentIndex = findspecies( mole_global.list[i]->parentLabel.c_str() )->index;
00307 mole_global.list[i]->parentIndex = parentIndex;
00308 }
00309 else
00310 mole_global.list[i]->parentIndex = -1;
00311 }
00312
00313
00314 for(i=0;i<mole_global.num_total;i++)
00315 {
00316 molecule *sp = &(*mole_global.list[i]);
00317 if (sp->isMonatomic())
00318 {
00319 ASSERT( (int)sp->nAtom.size() == 1 );
00320 count_ptr<chem_atom> atom = sp->nAtom.begin()->first;
00321 ASSERT(sp->charge <= atom->el->Z);
00322 if(sp->charge >= 0 && sp->lgGas_Phase)
00323 {
00324 atom->ipMl[sp->charge] = i;
00325 }
00326 }
00327 }
00328
00329 return;
00330
00331 }
00332
00333 STATIC void read_species_file( string filename, bool lgCreateIsotopologues )
00334 {
00335 DEBUG_ENTRY( "read_sepcies_file()" );
00336
00337 fstream ioDATA;
00338 open_data( ioDATA, filename.c_str(), mode_r );
00339 string line;
00340
00341 while( getline( ioDATA,line ) )
00342 {
00343 if( line.empty() )
00344 break;
00345 if( line[0] == '#' )
00346 continue;
00347 istringstream iss( line );
00348 string species;
00349 double formation_enthalpy;
00350 iss >> species;
00351 iss >> formation_enthalpy;
00352 ASSERT( iss.eof() );
00353 newspecies( species.c_str(), MOLECULE,MOLE_ACTIVE, formation_enthalpy, lgCreateIsotopologues );
00354
00355 }
00356
00357 return;
00358 }
00359
00360
00361 void create_isotopologues_one(
00362 ChemAtomList& atoms,
00363 vector< int >& numAtoms,
00364 string atom_old,
00365 string atom_new,
00366 string embellishments,
00367 vector<string>& newLabels )
00368 {
00369 fixit();
00370 fixit();
00371
00372
00373 for( unsigned position = 0; position < atoms.size(); ++position )
00374 {
00375 stringstream str;
00376 if( atoms[position]->label() == atom_old )
00377 {
00378 for( unsigned i=0; i<position; ++i )
00379 {
00380 str << atoms[i]->label();
00381 if( numAtoms[i]>1 )
00382 str << numAtoms[i];
00383 }
00384
00385 if( numAtoms[position] > 1 )
00386 {
00387 str << atom_old;
00388 if( numAtoms[position] > 2 )
00389 str << numAtoms[position]-1;
00390 }
00391
00392
00393 if( position+1 == atoms.size() )
00394 str << atom_new;
00395
00396 for( unsigned i=position+1; i<atoms.size(); ++i )
00397 {
00398 if( i==position+1 )
00399 {
00400
00401 if( atom_new == atoms[i]->label() )
00402 {
00403 str << atoms[i]->label();
00404 ASSERT( numAtoms[i] + 1 > 1 );
00405 str << numAtoms[i] + 1;
00406 }
00407 else
00408 {
00409 str << atom_new;
00410 str << atoms[i]->label();
00411 if( numAtoms[i] > 1 )
00412 str << numAtoms[i];
00413 }
00414 }
00415 else
00416 {
00417 str << atoms[i]->label();
00418 if( numAtoms[i] > 1 )
00419 str << numAtoms[i];
00420 }
00421 }
00422
00423
00424 str << embellishments;
00425
00426 newLabels.push_back( str.str() );
00427
00428 }
00429 }
00430
00431 return;
00432 }
00433
00434
00435 STATIC void newelement(const char label[], int ipion)
00436 {
00437 char *s;
00438
00439 DEBUG_ENTRY("newelement()");
00440
00441
00442 int len = strlen(label)+1;
00443 auto_vec<char> mylab_v(new char[len]);
00444 char *mylab = mylab_v.data();
00445 strncpy(mylab,label,len);
00446 s = strchr(mylab,' ');
00447 if (s)
00448 *s = '\0';
00449
00450 int exists = (mole_priv::elemtab.find(mylab) != mole_priv::elemtab.end());
00451 if (!exists)
00452 {
00453 count_ptr<chem_element> element(new chem_element(ipion+1,mylab));
00454 mole_priv::elemtab[element->label] = element;
00455 element_list.push_back(element);
00456 }
00457 return;
00458 }
00459
00460
00461 STATIC void newisotope( const count_ptr<chem_element>& el, int massNumberA,
00462 realnum mass_amu, double frac )
00463 {
00464
00465 DEBUG_ENTRY("newisotope()");
00466
00467 ASSERT( massNumberA < 3 * LIMELM && ( massNumberA > 0 || massNumberA == -1 ) );
00468 ASSERT( mass_amu < 3. * LIMELM && mass_amu > 0. );
00469 ASSERT( frac <= 1. + FLT_EPSILON && frac >= 0. );
00470
00471 count_ptr<chem_atom> isotope(new chem_atom);
00472 isotope->A = massNumberA;
00473 isotope->mass_amu = mass_amu;
00474 isotope->frac = frac;
00475 isotope->ipMl.resize(el->Z+1);
00476 isotope->el = el.get_ptr();
00477 for (long int ion = 0; ion < el->Z+1; ion++)
00478 isotope->ipMl[ion] = -1;
00479
00480
00481 mole_priv::atomtab[ isotope->label() ] = isotope;
00482 atom_list.push_back(isotope);
00483 if( isotope->lgMeanAbundance() )
00484 unresolved_atom_list.push_back(isotope);
00485
00486 el->isotopes[massNumberA] = isotope;
00487 }
00488
00489 STATIC realnum MeanMassOfElement( const count_ptr<chem_element>& el )
00490 {
00491 DEBUG_ENTRY("MeanMassOfElement()");
00492
00493
00494 if( el->isotopes.size()==0 )
00495 return dense.AtomicWeight[el->Z-1];
00496
00497 realnum averageMass = 0., fracsum = 0.;
00498 for( isotopes_i it = el->isotopes.begin(); it != el->isotopes.end(); ++it )
00499 {
00500 fracsum += it->second->frac;
00501 averageMass += it->second->mass_amu * it->second->frac;
00502 }
00503 ASSERT( fp_equal( fracsum, 1.f ) );
00504
00505 return averageMass;
00506 }
00507
00508 STATIC molecule *newspecies(
00509 const char label[],
00510 enum spectype type,
00511 enum mole_state state,
00512 realnum form_enthalpy)
00513 {
00514 return newspecies( label, type, state, form_enthalpy, true);
00515 }
00516
00517
00518 STATIC molecule *newspecies(
00519 const char label[],
00520 enum spectype type,
00521 enum mole_state state,
00522 realnum form_enthalpy,
00523 bool lgCreateIsotopologues )
00524 {
00525 DEBUG_ENTRY("newspecies()");
00526
00527 int exists;
00528 ChemAtomList atomsLeftToRight;
00529 vector< int > numAtoms;
00530 string embellishments;
00531 char *s;
00532 count_ptr<molecule> mol(new molecule);
00533
00534 mol->parentLabel.clear();
00535 mol->isEnabled = true;
00536 mol->charge = 0;
00537 mol->lgExcit = false;
00538 mol->mole_mass = 0.0;
00539 mol->state = state;
00540 mol->lgGas_Phase = true;
00541 mol->form_enthalpy = form_enthalpy;
00542 mol->groupnum = -1;
00543
00544
00545 int len = strlen(label)+1;
00546 auto_vec<char> mylab_v(new char[len]);
00547 char *mylab = mylab_v.data();
00548 strncpy(mylab,label,len);
00549 s = strchr(mylab,' ');
00550 if (s)
00551 *s = '\0';
00552 mol->label = mylab;
00553
00554 if(type == MOLECULE)
00555 {
00556 if( parse_species_label( mylab, atomsLeftToRight, numAtoms, embellishments, mol->lgExcit, mol->charge, mol->lgGas_Phase ) == false )
00557 return NULL;
00558
00559 for( unsigned i = 0; i < atomsLeftToRight.size(); ++i )
00560 {
00561 mol->nAtom[ atomsLeftToRight[i] ] += numAtoms[i];
00562 mol->mole_mass += numAtoms[i] * atomsLeftToRight[i]->mass_amu * (realnum)ATOMIC_MASS_UNIT;
00563 }
00564 }
00565
00566
00567
00568
00569
00570 if ( (mol->n_nuclei() > 1 || (mol->isMonatomic() && mol->charge==-1) ) && mole_global.lgNoMole)
00571 {
00572 if( trace.lgTraceMole )
00573 fprintf(ioQQQ,"No species %s as molecules off\n",label);
00574 return NULL;
00575 }
00576
00577 if (mol->n_nuclei() > 1 && mole_global.lgNoHeavyMole)
00578 {
00579 for( nAtoms_ri it=mol->nAtom.rbegin(); it != mol->nAtom.rend(); --it )
00580 {
00581 if( it->first->el->Z-1 != ipHYDROGEN )
00582 {
00583 ASSERT( it->second > 0 );
00584 if( trace.lgTraceMole )
00585 fprintf(ioQQQ,"No species %s as heavy molecules off\n",label);
00586 return NULL;
00587 }
00588 }
00589 }
00590
00591
00592 exists = (mole_priv::spectab.find(mol->label) != mole_priv::spectab.end());
00593 if( exists )
00594 {
00595 fprintf( ioQQQ,"Warning: duplicate species %s - using first one\n",
00596 mol->label.c_str() );
00597 return NULL;
00598 }
00599 else
00600 mole_priv::spectab[mol->label] = mol;
00601
00602
00603 for( nAtoms_i it=mol->nAtom.begin(); it != mol->nAtom.end(); ++it )
00604 ASSERT( it->second > 0 );
00605
00606 if (state != MOLE_NULL)
00607 mole_global.num_total++;
00608 if (state == MOLE_ACTIVE)
00609 mole_global.num_calc++;
00610
00611
00612 if( lgCreateIsotopologues && type==MOLECULE && mol->label.compare("CO")==0 )
00613 {
00614 molecule *sp = newspecies( "^13CO", MOLECULE, mol->state, mol->form_enthalpy, false );
00615 sp->parentLabel = mol->label;
00616 }
00617
00618
00619 if( lgCreateIsotopologues && type==MOLECULE && !mol->isMonatomic() )
00620 {
00621 for( nAtoms_i it1 = mol->nAtom.begin(); it1 != mol->nAtom.end(); ++it1 )
00622 {
00623 for( map<int, count_ptr<chem_atom> >::iterator it2 = it1->first->el->isotopes.begin();
00624 it2 != it1->first->el->isotopes.end(); ++it2 )
00625 {
00626
00627 if( it1->first->el->Z-1 == ipHYDROGEN && it2->second->A != 2 )
00628 continue;
00629 if( !mole_global.lgTreatIsotopes[it1->first->el->Z-1] )
00630 continue;
00631 if( it2->second->lgMeanAbundance() )
00632 continue;
00633 vector<string> isoLabs;
00634 create_isotopologues_one( atomsLeftToRight, numAtoms, it1->first->label(), it2->second->label(), embellishments, isoLabs );
00635
00636
00637
00638
00639 for( vector<string>::iterator newLabel = isoLabs.begin(); newLabel != isoLabs.end(); ++newLabel )
00640 {
00641 molecule *sp = newspecies( newLabel->c_str(), MOLECULE, mol->state, mol->form_enthalpy, false );
00642
00643 if( sp!=NULL && it2->second->A != 2 )
00644 sp->parentLabel = mol->label;
00645 }
00646 }
00647 }
00648 }
00649
00650 return &(*mol);
00651 }
00652 bool parse_species_label( const char label[], ChemAtomList &atomsLeftToRight, vector<int> &numAtoms, string &embellishments )
00653 {
00654 bool lgExcit, lgGas_Phase;
00655 int charge;
00656 bool lgOK = parse_species_label( label, atomsLeftToRight, numAtoms, embellishments, lgExcit, charge, lgGas_Phase );
00657 return lgOK;
00658 }
00659 bool parse_species_label( const char label[], ChemAtomList &atomsLeftToRight, vector<int> &numAtoms, string &embellishments,
00660 bool &lgExcit, int &charge, bool &lgGas_Phase )
00661 {
00662 long int i, n, ipAtom;
00663 char thisAtom[CHARS_ISOTOPE_SYM];
00664 count_ptr<chem_atom> atom;
00665 char mylab[CHARS_SPECIES];
00666 char *s;
00667
00668 strncpy( mylab, label, CHARS_SPECIES );
00669
00670
00671 s = strpbrk(mylab,"*");
00672 if(s)
00673 {
00674 lgExcit = true;
00675 embellishments = s;
00676 *s = '\0';
00677 }
00678
00679
00680 s = strpbrk(mylab,"+-");
00681 if(s)
00682 {
00683 if(isdigit(*(s+1)))
00684 n = atoi(s+1);
00685 else
00686 n = 1;
00687 if(*s == '+')
00688 charge = n;
00689 else
00690 charge = -n;
00691 embellishments = s + embellishments;
00692 *s = '\0';
00693 }
00694
00695 s = strstr(mylab,"grn");
00696 if(s)
00697 {
00698 lgGas_Phase = false;
00699 embellishments = s + embellishments;
00700 *s = '\0';
00701 }
00702 else
00703 {
00704 lgGas_Phase = true;
00705 }
00706
00707
00708
00709 i = 0;
00710 while (mylab[i] != '\0' && mylab[i] != ' ' && mylab[i] != '*')
00711 {
00712
00713 ipAtom = 0;
00714
00715 if(mylab[i]=='^')
00716 {
00717 thisAtom[ipAtom++] = mylab[i++];
00718 ASSERT( isdigit(mylab[i]) );
00719 thisAtom[ipAtom++] = mylab[i++];
00720 if(isdigit(mylab[i]))
00721 {
00722 thisAtom[ipAtom++] = mylab[i++];
00723 }
00724 }
00725
00726 thisAtom[ipAtom++] = mylab[i++];
00727 if(islower(mylab[i]))
00728 {
00729 thisAtom[ipAtom++] = mylab[i++];
00730 }
00731 thisAtom[ipAtom] = '\0';
00732 ASSERT(ipAtom <= CHARS_ISOTOPE_SYM);
00733
00734 atom = findatom(thisAtom);
00735 if(atom.get_ptr() == NULL)
00736 {
00737 fprintf(stderr,"Did not recognize atom at %s in \"%s \"[%ld]\n",thisAtom,mylab,i);
00738 exit(-1);
00739 }
00740 if(!dense.lgElmtOn[atom->el->Z-1])
00741 {
00742 if( trace.lgTraceMole )
00743 fprintf(ioQQQ,"No species %s as element %s off\n",mylab,atom->el->label.c_str() );
00744 return false;
00745 }
00746
00747 if(isdigit(mylab[i]))
00748 {
00749 n = 0;
00750 do {
00751 n = 10*n+(long int)(mylab[i]-'0');
00752 i++;
00753 } while (isdigit(mylab[i]));
00754 }
00755 else
00756 {
00757 n = 1;
00758 }
00759 atomsLeftToRight.push_back( atom );
00760 numAtoms.push_back( n );
00761 }
00762
00763 return true;
00764 }
00765 STATIC bool isactive(const molecule &mol)
00766 {
00767 DEBUG_ENTRY("isactive()");
00768 return mol.state == MOLE_ACTIVE;
00769 }
00770 STATIC bool ispassive(const molecule &mol)
00771 {
00772
00773 DEBUG_ENTRY("ispassive()");
00774 return mol.state == MOLE_PASSIVE;
00775 }
00776
00777 bool lgDifferByExcitation( const molecule &mol1, const molecule &mol2 )
00778 {
00779 if( mol1.label == mol2.label + "*" )
00780 return true;
00781 else if( mol2.label == mol1.label + "*" )
00782 return true;
00783 else
00784 return false;
00785 }
00786
00787 molecule *findspecies(const char buf[])
00788 {
00789 DEBUG_ENTRY("findspecies()");
00790
00791
00792 string s;
00793 for (const char *pb = buf; *pb && *pb != ' '; ++pb)
00794 {
00795 s += *pb;
00796 }
00797
00798 const molecule_i p = mole_priv::spectab.find(s);
00799
00800 if(p != mole_priv::spectab.end())
00801 return &(*p->second);
00802 else
00803 return null_mole;
00804 }
00805
00806 molezone *findspecieslocal(const char buf[])
00807 {
00808 DEBUG_ENTRY("findspecieslocal()");
00809
00810
00811 string s;
00812 for (const char *pb = buf; *pb && *pb != ' '; ++pb)
00813 {
00814 s += *pb;
00815 }
00816
00817 const molecule_i p = mole_priv::spectab.find(s);
00818
00819 if(p != mole_priv::spectab.end())
00820 return &mole.species[ p->second->index ];
00821 else
00822 return null_molezone;
00823 }
00824
00825 STATIC count_ptr<chem_atom> findatom(const char buf[])
00826 {
00827 chem_atom_i p;
00828
00829 DEBUG_ENTRY("findatom()");
00830
00831 p = mole_priv::atomtab.find(buf);
00832
00833 if(p != mole_priv::atomtab.end())
00834 return p->second;
00835 else
00836 return count_ptr<chem_atom>(NULL);
00837 }
00838
00839 void mole_update_species_cache(void)
00840 {
00841 int i;
00842 double den_times_area, den_grains, adsorbed_density;
00843
00844 DEBUG_ENTRY("mole_update_species_cache()");
00845
00846 enum { DEBUG_MOLE = false };
00847
00848 for (i=0;i<mole_global.num_total;i++)
00849 {
00850 if(mole.species[i].location != NULL)
00851 {
00852 ASSERT( mole_global.list[i]->parentLabel.empty() );
00853 mole.species[i].den = *(mole.species[i].location);
00854 if (DEBUG_MOLE)
00855 fprintf(ioQQQ,"%s: %f\n",mole_global.list[i]->label.c_str(),mole.species[i].den);
00856 }
00857 }
00858
00859 mole.set_isotope_abundances();
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870 if( gv.lgDustOn() )
00871 {
00872 den_times_area = den_grains = 0.0;
00873 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00874 {
00875
00876 den_times_area += gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
00877 den_grains += gv.bin[nd]->cnv_GR_pCM3;
00878 }
00879
00880 adsorbed_density = 0.0;
00881 for (i=0;i<mole_global.num_total;i++)
00882 {
00883 if( !mole_global.list[i]->lgGas_Phase && mole_global.list[i]->parentLabel.empty() )
00884 adsorbed_density += mole.species[i].den;
00885 }
00886
00887 mole.grain_area = (realnum) den_times_area;
00888 mole.grain_density = (realnum) den_grains;
00889
00890 double mole_cs = 1e-15;
00891 if (4*den_times_area <= mole_cs*adsorbed_density)
00892 mole.grain_saturation = 1.0;
00893 else
00894 mole.grain_saturation = (realnum)((mole_cs*adsorbed_density)/(4.*den_times_area));
00895 }
00896 else
00897 {
00898 mole.grain_area = 0.0;
00899 mole.grain_density = 0.0;
00900 mole.grain_saturation = 1.0;
00901 }
00902 if (DEBUG_MOLE)
00903 fprintf(ioQQQ,"Dust: %f %f %f\n",
00904 mole.grain_area,mole.grain_density,mole.grain_saturation);
00905
00906 }
00907
00908 realnum mole_return_cached_species(void)
00909 {
00910 int i;
00911
00912
00913
00914
00915 ASSERT(lgElemsConserved());
00916
00917
00918 dense.updateXMolecules();
00919 total_molecule_deut( deut.xMolecules );
00920
00921
00922 mole.elec = 0.;
00923 for(i=0;i<mole_global.num_calc;i++)
00924 {
00925 if (mole.species[i].location == NULL && mole_global.list[i]->parentLabel.empty())
00926 mole.elec += mole.species[i].den*mole_global.list[i]->charge;
00927 }
00928
00929
00930 realnum delta = 0.0;
00931 long ncpt = 0;
00932
00933 for (i=0;i<mole_global.num_total;i++)
00934 {
00935 if(mole.species[i].location && mole_global.list[i]->state == MOLE_ACTIVE)
00936 {
00937 realnum old_pop = *(mole.species[i].location),
00938 new_pop = mole.species[i].den;
00939 realnum frac = (new_pop-old_pop)/SDIV(new_pop+old_pop);
00940
00941 if( mole_global.list[i]->isMonatomic() )
00942 {
00943 delta += frac*frac;
00944 ++ncpt;
00945 }
00946
00947
00948
00949
00950 *(mole.species[i].location) = new_pop;
00951 }
00952 }
00953
00954
00955 ASSERT(lgElemsConserved());
00956 return ncpt>0 ? sqrt(delta/ncpt) : 0.f;
00957 }
00958
00959 void t_mole_local::set_isotope_abundances( void )
00960 {
00961 DEBUG_ENTRY("t_mole_local::set_isotope_abundances()");
00962
00963
00964 for(ChemAtomList::iterator atom = unresolved_atom_list.begin(); atom != unresolved_atom_list.end(); ++atom)
00965 {
00966
00967 for( isotopes_i it = (*atom)->el->isotopes.begin(); it != (*atom)->el->isotopes.end(); ++it )
00968 {
00969
00970 if( it->second->lgMeanAbundance() )
00971 continue;
00972
00973
00974 for( unsigned long ion=0; ion<it->second->ipMl.size(); ++ion )
00975 {
00976 if( it->second->ipMl[ion] != -1 &&
00977 (species[ it->second->ipMl[ion] ].location == NULL ) && (*atom)->ipMl[ion] != -1 )
00978 {
00979 species[ it->second->ipMl[ion] ].den = species[ (*atom)->ipMl[ion] ].den * it->second->frac;
00980 }
00981 }
00982 }
00983 }
00984
00985 return;
00986 }
00987
00988 void t_mole_local::set_location( long nelem, long ion, double *density )
00989 {
00990 DEBUG_ENTRY( "t_mole_local::set_location()" );
00991
00992 ASSERT( nelem < LIMELM );
00993 ASSERT( ion < nelem + 2 );
00994 long mole_index = unresolved_atom_list[nelem]->ipMl[ion];
00995
00996 if( mole_index == -1 )
00997 return;
00998 ASSERT( mole_index < mole_global.num_total );
00999 species[mole_index].location = density;
01000
01001 return;
01002 }
01003
01004 void total_molecule_deut( realnum &total_f )
01005 {
01006 DEBUG_ENTRY( "total_molecule_deut()" );
01007
01008 double total = 0.;
01009
01010 if( !deut.lgElmtOn )
01011 return;
01012
01013 for (long int i=0;i<mole_global.num_calc;++i)
01014 {
01015 if (mole.species[i].location == NULL && mole_global.list[i]->parentLabel.empty() )
01016 {
01017 for( molecule::nAtomsMap::iterator atom=mole_global.list[i]->nAtom.begin();
01018 atom != mole_global.list[i]->nAtom.end(); ++atom)
01019 {
01020 long int nelem = atom->first->el->Z-1;
01021 if( nelem==0 && atom->first->A==2 )
01022 {
01023 total += mole.species[i].den*atom->second;
01024 }
01025 }
01026 }
01027 }
01028
01029 total_f = (realnum)total;
01030
01031 return;
01032 }
01033 void total_molecule_elems(realnum total[LIMELM])
01034 {
01035 DEBUG_ENTRY( "total_molecule_elems()" );
01036
01037
01038 for ( long int nelem=ipHYDROGEN;nelem<LIMELM; ++nelem )
01039 {
01040 total[nelem] = 0.;
01041 }
01042 for (long int i=0;i<mole_global.num_calc;++i)
01043 {
01044 if (mole.species[i].location == NULL && mole_global.list[i]->parentLabel.empty() )
01045 {
01046 for( molecule::nAtomsMap::iterator atom=mole_global.list[i]->nAtom.begin();
01047 atom != mole_global.list[i]->nAtom.end(); ++atom)
01048 {
01049 ASSERT( atom->second > 0 );
01050 long int nelem = atom->first->el->Z-1;
01051 if( atom->first->lgMeanAbundance() )
01052 total[nelem] += (realnum) mole.species[i].den*atom->second;
01053 }
01054 }
01055 }
01056 }
01057 void total_network_elems(double total[LIMELM])
01058 {
01059 DEBUG_ENTRY( "total_network_elems()" );
01060
01061
01062 for ( long int nelem=ipHYDROGEN;nelem<LIMELM; ++nelem )
01063 {
01064 total[nelem] = 0.;
01065 }
01066 for (long int i=0;i<mole_global.num_calc;++i)
01067 {
01068 if (mole_global.list[i]->parentLabel.empty())
01069 {
01070 for( molecule::nAtomsMap::iterator atom=mole_global.list[i]->nAtom.begin();
01071 atom != mole_global.list[i]->nAtom.end(); ++atom)
01072 {
01073 long int nelem = atom->first->el->Z-1;
01074 total[nelem] += (realnum) mole.species[i].den*atom->second;
01075 }
01076 }
01077 }
01078 }
01079 realnum total_molecules(void)
01080 {
01081 long int i;
01082 realnum total;
01083
01084 DEBUG_ENTRY( "total_molecules()" );
01085
01086 total = 0.;
01087 for (i=0;i<mole_global.num_calc;++i)
01088 {
01089 if (mole.species[i].location == NULL && mole_global.list[i]->parentLabel.empty())
01090 total += (realnum) mole.species[i].den;
01091 }
01092 return total;
01093 }
01094 realnum total_molecules_gasphase(void)
01095 {
01096 long int i;
01097 realnum total;
01098
01099 DEBUG_ENTRY( "total_molecules_gasphase()" );
01100
01101 total = 0.;
01102 for (i=0;i<mole_global.num_calc;++i)
01103 {
01104 if (mole_global.list[i]->lgGas_Phase && mole.species[i].location == NULL && mole_global.list[i]->parentLabel.empty())
01105 total += (realnum) mole.species[i].den;
01106 }
01107 return total;
01108 }
01109
01110
01111
01112 void mole_make_groups(void)
01113 {
01114 long int i, j;
01115
01116
01117
01118 DEBUG_ENTRY ("mole_make_groups()");
01119 if (mole_global.num_calc == 0)
01120 {
01121 groupspecies = NULL;
01122 mole_global.num_compacted = 0;
01123 return;
01124 }
01125 groupspecies = (molecule **) MALLOC(mole_global.num_calc*sizeof(molecule *));
01126 for (i=0,j=0;i<mole_global.num_calc;i++)
01127 {
01128 if( mole_global.list[i]->parentLabel.empty() && ( !mole_global.list[i]->isMonatomic() || mole_global.list[i]->charge <= 0 || ! mole_global.list[i]->lgGas_Phase ) )
01129 {
01130
01131 mole_global.list[i]->groupnum = j;
01132 groupspecies[j++] = &(*mole_global.list[i]);
01133 }
01134 else
01135 {
01136
01137
01138 ASSERT (mole_global.list[i]->charge < LIMELM+1);
01139 ASSERT (mole_global.list[i]->groupnum == -1);
01140 }
01141 }
01142 mole_global.num_compacted = j;
01143 groupspecies = (molecule **) REALLOC((void *)groupspecies,
01144 mole_global.num_compacted*sizeof(molecule *));
01145
01146 for (i=0;i<mole_global.num_calc;i++)
01147 {
01148 if (mole_global.list[i]->groupnum == -1)
01149 {
01150 if( mole_global.list[i]->isMonatomic() && mole_global.list[i]->parentLabel.empty() )
01151 {
01152 for( nAtoms_i it = mole_global.list[i]->nAtom.begin(); it != mole_global.list[i]->nAtom.end(); ++it )
01153 {
01154 ASSERT( it->second > 0 );
01155 mole_global.list[i]->groupnum = mole_global.list[it->first->ipMl[0]]->groupnum;
01156 break;
01157 }
01158 }
01159 else
01160 {
01161 ASSERT( !mole_global.list[i]->parentLabel.empty() );
01162 mole_global.list[i]->groupnum = mole_global.list[ mole_global.list[i]->parentIndex ]->groupnum;
01163 }
01164 }
01165
01166 ASSERT( mole_global.list[i]->groupnum != -1);
01167 }
01168
01169 return;
01170 }
01171