cloudy trunk
Loading...
Searching...
No Matches
mole_species.cpp
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/*CO_Init called from cdInit to initialize co routines */
4/*CO_update_chem_rates update rate coefficients, only temp part - in mole_co_etc.c */
5#include "cdstd.h"
6#include <cctype>
7#include <string.h>
8#include <algorithm>
9#include <stdlib.h>
10#include "cddefines.h"
11#include "co.h"
12#include "colden.h"
13#include "conv.h"
14#include "deuterium.h"
15#include "doppvel.h"
16#include "elementnames.h"
17#include "h2.h"
18#include "iso.h"
19#include "phycon.h"
20#include "physconst.h"
21#include "mole.h"
22#include "mole_priv.h"
23#include "hmi.h"
24#include "radius.h"
25#include "rfield.h"
26#include "rt.h"
27#include "secondaries.h"
28#include "dense.h"
29#include "ionbal.h"
30#include "grainvar.h"
31#include "timesc.h"
32#include "taulines.h"
33#include "trace.h"
34/*lint -e778 constant expression evaluates to 0 in operation '-' */
35
36/* CO_update_chem_rates update rate coefficients, only temp part - in
37 * mole_co_etc.c called in conv_base before any chemistry or
38 * ionization is done */
39
41
42STATIC void read_species_file( string filename, bool lgCreateIsotopologues = true );
43STATIC void newelement(const char label[], int ipion);
44STATIC void newisotope( const count_ptr<chem_element> &el, int massNumberA,
45 realnum mass_amu, double frac );
47// newspecies is overloaded. The first one just calls the second with the additional lgCreateIsotopologues set true
48STATIC molecule *newspecies(const char label[], enum spectype type,
49 enum mole_state state, realnum form_enthalpy);
50STATIC molecule *newspecies(const char label[], enum spectype type,
51 enum mole_state state, realnum form_enthalpy, bool lgCreateIsotopologues);
52STATIC count_ptr<chem_atom> findatom(const char buf[]);
53STATIC bool isactive(const molecule &mol);
54STATIC bool ispassive(const molecule &mol);
55STATIC void ReadIsotopeFractions( const vector<bool>& lgResolveNelem );
56//STATIC void create_isotopologues(count_ptr<molecule> mol, ChemAtomList& atoms, vector<int>& numAtoms);
57
58namespace mole_priv {
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;
64};
69vector< count_ptr <chem_element> > element_list;
73
74#include <functional>
75
76namespace
77{
78 class MoleCmp : public binary_function<const count_ptr<molecule>,
79 const count_ptr<molecule>,bool>
80 {
81 public:
82 bool operator()(const count_ptr<molecule> &mol1,
83 const count_ptr<molecule> &mol2) const
84 {
85 return mol1->compare(*mol2) < 0;
86 }
87 bool operator()(const molecule *mol1, const molecule *mol2) const
88 {
89 return mol1->compare(*mol2) < 0;
90 }
91 };
92}
93
94void t_mole_global::sort(t_mole_global::MoleculeList::iterator start,
95 t_mole_global::MoleculeList::iterator end)
96{
97 std::sort(start,end,MoleCmp());
98}
100{
101 std::sort(start,end,MoleCmp());
102}
103
104STATIC void ReadIsotopeFractions( const vector<bool>& lgResolveNelem )
105{
106 DEBUG_ENTRY( "ReadIsotopeFractions()" );
107
108 char chLine[INPUT_LINE_LENGTH];
109 long i;
110 bool lgEOL;
111
112 static const char chFile[] = "isotope_fractions.dat";
113 FILE *ioDATA = open_data( chFile, "r" );
114 ASSERT( ioDATA != NULL );
115
116 while( read_whole_line( chLine, (int)sizeof(chLine), ioDATA ) != NULL )
117 {
118 if( chLine[0] == '#' )
119 continue;
120
121 i = 1;
122 long Z = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
123 long A = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
124 // file has fractionation as percent; the 0.01 converts to fraction
125 double frac = 0.01 * (double)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
126
127 fixit(); // need real masses here instead of just A (protons + neutrons)
128 if( (unsigned)Z <= lgResolveNelem.size() && lgResolveNelem[Z-1] )
129 newisotope( element_list[Z-1], A, (realnum)A, frac );
130 // always do this to continue history of predicting 13CO
131 else if( Z-1==ipCARBON )
132 newisotope( element_list[Z-1], A, (realnum)A, frac );
133 }
134
135 fclose( ioDATA );
136
137 return;
138}
139
141{
142 DEBUG_ENTRY( "mole_global::make_species()" );
143
144 long int i;
145 molecule *sp;
146
147 null_element = new chem_element(0,"") ;
148 null_atom = new( chem_atom );
149 null_molezone = new( molezone );
150 null_molezone->den = 0.;
151
152 /* set up concordance of elemental species to external Cloudy indices */
153 for (i=0;i<LIMELM;i++)
154 {
155 newelement(elementnames.chElementSym[i], i);
156 }
157
158 // flip this to treat deuterium
159 if( deut.lgElmtOn )
161
162 // read and define isotopes, set default fractionations
164
165 // special code to maintain effect of "set 12C13C" command
166 {
167 count_ptr<chem_atom> atom12C = findatom("^12C");
168 count_ptr<chem_atom> atom13C = findatom("^13C");
169
170 if( co.C12_C13_isotope_ratio_parsed >= 0. )
171 {
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;
176 }
177 else
178 {
179 co.C12_C13_isotope_ratio = atom12C->frac / SDIV( atom13C->frac );
180 }
181 // Make sure only two isotopes are defined: 12C, 13C.
182 // The mean-abundance isotope is defined below.
183 ASSERT( element_list[ipCARBON]->isotopes.size() == 2 );
184 }
185
186
188 {
192 }
193
194 for( long nelem=0; nelem<LIMELM; ++nelem )
195 {
196 realnum mass_amu = MeanMassOfElement( element_list[nelem] );
197 //define generic mean-abundance isotopes
198 newisotope( element_list[nelem], -1, mass_amu, 1.0 );
199 }
200
201 ASSERT( (long) unresolved_atom_list.size() == LIMELM );
202
203 /* set up properties of molecular species -- chemical formulae,
204 array indices, elementary components (parsed from formula),
205 status within CO network, location of stored value external
206 to CO network (must be floating point). */
207
208 /* Sizes of different parts of network are calculated by increments in newspecies */
209 mole_global.num_total = mole_global.num_calc = 0;
210 /* Enthalpies of formation taken from
211 * >>refer mol Le Teuff, Y. E., Millar, T. J., & Markwick, A. J.,2000, A&AS, 146, 157
212 */
213
214 /* Zero density pseudo-species to return when molecule is switched off */
216 null_mole->index = -1;
217
218 read_species_file( "chem_species.dat" );
219
220 if(gv.lgDustOn() && mole_global.lgGrain_mole_deplete )
221 {
222 /* What are formation enthalpies of COgrn, H2Ogrn, OHgrn? For
223 present, take grn as standard state, and neglect adsorption enthalpy.
224
225 -- check e.g. http://www.arxiv.org/abs/astro-ph/0702322 for CO adsorption energy.
226 */
227 if (0)
228 {
229 read_species_file( "chem_species_grn.dat" );
230 }
231 else
232 {
233 sp = newspecies("COgrn ",MOLECULE,MOLE_ACTIVE, -113.8f);
234 sp = newspecies("H2Ogrn",MOLECULE,MOLE_ACTIVE, -238.9f);
235 sp = newspecies("OHgrn ",MOLECULE,MOLE_ACTIVE, 38.4f);
236 //sp = newspecies("Hgrn ",MOLECULE,MOLE_ACTIVE, 216.f);
237 }
238 }
239
240 /* Add passive species to complete network */
241 sp = newspecies("e- ",OTHER,MOLE_PASSIVE, 0.0f);
242 sp->charge = -1; sp->mole_mass = (realnum)ELECTRON_MASS; /* Augment properties for this non-molecular species */
243 sp = newspecies("grn ",OTHER,MOLE_PASSIVE, 0.0f);
244 sp = newspecies("PHOTON",OTHER,MOLE_PASSIVE, 0.0f);
245 sp = newspecies("CRPHOT",OTHER,MOLE_PASSIVE, 0.0f);
246 sp = newspecies("CRP ",OTHER,MOLE_PASSIVE, 0.0f);
247
248 if (!mole_global.lgLeidenHack)
249 sp = newspecies("H- ",MOLECULE,MOLE_ACTIVE, 143.2f);
250 if (hmi.lgLeiden_Keep_ipMH2s)
251 {
252 sp = newspecies("H2* ",MOLECULE,MOLE_ACTIVE,
253 h2.ENERGY_H2_STAR * KJMOL1CM);
254 }
255
256 if( deut.lgElmtOn )
257 {
258 read_species_file( "chem_species_deuterium.dat", false );
259 }
260
261 /* Add species for all other elements and their first ions -- couple to network at least via H- */
262 for( ChemAtomList::iterator atom = unresolved_atom_list.begin();
263 atom != unresolved_atom_list.end(); ++atom)
264 {
265 long int nelem = (*atom)->el->Z-1;
266
267 for( long ion=0; ion<=nelem+1; ion++ )
268 {
270 char temp[CHARS_ION_STAGE+1];
271 if( ion==0 )
272 temp[0] = '\0';
273 else if( ion==1 )
274 sprintf( temp, "+" );
275 else
276 sprintf( temp, "+%ld", ion );
277 sprintf( str, "%s%s", (*atom)->label().c_str(), temp );
278 if (findspecies(str) == null_mole)
279 {
280 sp = newspecies(str,MOLECULE,MOLE_ACTIVE, 0.f);
281 fixit(); // populate these in a local update
282#if 0
283 if( sp != NULL )
284 {
285 sp->levels = NULL;
286 sp->numLevels = 0;
287 }
288#endif
289 }
290 }
291 }
292
293 return;
294}
295
297{
298 DEBUG_ENTRY( "mole_make_list()" );
299
300 /* Create linear list of species and populate it... */
301 mole_global.list.resize(mole_global.num_total);
302
303 /* ...first active species */
304 long int i = 0;
305 for (molecule_i p = mole_priv::spectab.begin(); p != mole_priv::spectab.end(); ++p)
306 {
307 if (isactive(*(p->second)))
308 mole_global.list[i++] = p->second;
309 }
310 ASSERT (i == mole_global.num_calc);
311
312 /* Sort list into a standard ordering */
313 t_mole_global::sort(mole_global.list.begin(),mole_global.list.begin()+mole_global.num_calc);
314
315 for (molecule_i p = mole_priv::spectab.begin(); p != mole_priv::spectab.end(); ++p)
316 {
317 if (ispassive(*(p->second)))
318 mole_global.list[i++] = p->second;
319 }
320 ASSERT (i == mole_global.num_total);
321
322 /* Set molecule indices to order of list just created */
323 for(i=0;i<mole_global.num_total;i++)
324 {
325 mole_global.list[i]->index = i;
326 }
327
328 for(i=0;i<mole_global.num_total;i++)
329 {
330 if( !mole_global.list[i]->parentLabel.empty() )
331 {
332 long parentIndex = findspecies( mole_global.list[i]->parentLabel.c_str() )->index;
333 mole_global.list[i]->parentIndex = parentIndex;
334 }
335 else
336 mole_global.list[i]->parentIndex = -1;
337 }
338
339 /* Register the atomic ladders */
340 for(i=0;i<mole_global.num_total;i++)
341 {
342 molecule *sp = &(*mole_global.list[i]);
343 if (sp->isMonatomic())
344 {
345 ASSERT( (int)sp->nAtom.size() == 1 );
346 count_ptr<chem_atom> atom = sp->nAtom.begin()->first;
347 ASSERT(sp->charge <= atom->el->Z);
348 if(sp->charge >= 0 && sp->lgGas_Phase)
349 {
350 atom->ipMl[sp->charge] = i;
351 }
352 }
353 }
354
355 return;
356
357}
358
359STATIC void read_species_file( string filename, bool lgCreateIsotopologues )
360{
361 DEBUG_ENTRY( "read_sepcies_file()" );
362
363 fstream ioDATA;
364 open_data( ioDATA, filename.c_str(), mode_r );
365 string line;
366
367 while( getline( ioDATA,line ) )
368 {
369 if( line.empty() )
370 break;
371 if( line[0] == '#' )
372 continue;
373 istringstream iss( line );
374 string species;
375 double formation_enthalpy;
376 iss >> species;
377 iss >> formation_enthalpy;
378 ASSERT( iss.eof() );
379 newspecies( species.c_str(), MOLECULE,MOLE_ACTIVE, formation_enthalpy, lgCreateIsotopologues );
380 //fprintf( ioQQQ, "DEBUGGG read_species_file %s\t%f\n", species.c_str(), formation_enthalpy );
381 }
382
383 return;
384}
385/*lint +e778 constant expression evaluates to 0 in operation '-' */
386
389 vector< int >& numAtoms,
390 string atom_old,
391 string atom_new,
392 string embellishments,
393 vector<string>& newLabels )
394{
395 fixit(); // make sure atom_new and atom_old are isotopes
396 fixit(); // make sure atom_new is not already present
397
398 //for( ChemAtomList::iterator it = atoms.begin(); it != atoms.end(); ++it )
399 for( unsigned position = 0; position < atoms.size(); ++position )
400 {
401 stringstream str;
402 if( atoms[position]->label() == atom_old )
403 {
404 for( unsigned i=0; i<position; ++i )
405 {
406 str << atoms[i]->label();
407 if( numAtoms[i]>1 )
408 str << numAtoms[i];
409 }
410
411 if( numAtoms[position] > 1 )
412 {
413 str << atom_old;
414 if( numAtoms[position] > 2 )
415 str << numAtoms[position]-1;
416 }
417
418
419 if( position+1 == atoms.size() )
420 str << atom_new;
421
422 for( unsigned i=position+1; i<atoms.size(); ++i )
423 {
424 if( i==position+1 )
425 {
426 // remove adjacent duplicates
427 if( atom_new == atoms[i]->label() )
428 {
429 str << atoms[i]->label();
430 ASSERT( numAtoms[i] + 1 > 1 );
431 str << numAtoms[i] + 1;
432 }
433 else
434 {
435 str << atom_new;
436 str << atoms[i]->label();
437 if( numAtoms[i] > 1 )
438 str << numAtoms[i];
439 }
440 }
441 else
442 {
443 str << atoms[i]->label();
444 if( numAtoms[i] > 1 )
445 str << numAtoms[i];
446 }
447 }
448
449 // add on charge, grn, and excitation embellishments
450 str << embellishments;
451
452 newLabels.push_back( str.str() );
453 //fprintf( ioQQQ, "DEBUGGG create_isotopologues_one %s\n", newLabels.back().c_str() );
454 }
455 }
456
457 return;
458}
459
460/* Fill element linking structure */
461STATIC void newelement(const char label[], int ipion)
462{
463 char *s;
464
465 DEBUG_ENTRY("newelement()");
466
467 /* Create private workspace for label; copy and strip trailing whitespace */
468 int len = strlen(label)+1;
469 auto_vec<char> mylab_v(new char[len]);
470 char *mylab = mylab_v.data();
471 strncpy(mylab,label,len);
472 s = strchr(mylab,' ');
473 if (s)
474 *s = '\0';
475
476 int exists = (mole_priv::elemtab.find(mylab) != mole_priv::elemtab.end());
477 if (!exists)
478 {
479 count_ptr<chem_element> element(new chem_element(ipion+1,mylab));
480 mole_priv::elemtab[element->label] = element;
481 element_list.push_back(element);
482 }
483 return;
484}
485
486/* Fill isotope lists */
487STATIC void newisotope( const count_ptr<chem_element>& el, int massNumberA,
488 realnum mass_amu, double frac )
489{
490
491 DEBUG_ENTRY("newisotope()");
492
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. );
496
497 count_ptr<chem_atom> isotope(new chem_atom);
498 isotope->A = massNumberA;
499 isotope->mass_amu = mass_amu;
500 isotope->frac = frac;
501 isotope->ipMl.resize(el->Z+1);
502 isotope->el = el.get_ptr();
503 for (long int ion = 0; ion < el->Z+1; ion++)
504 isotope->ipMl[ion] = -1; /* Chemical network species indices not yet defined */
505
506 //int exists = (mole_priv::atomtab.find( isotope->label() ) != mole_priv::atomtab.end());
507 mole_priv::atomtab[ isotope->label() ] = isotope;
508 atom_list.push_back(isotope);
509 if( isotope->lgMeanAbundance() )
510 unresolved_atom_list.push_back(isotope);
511 // register with 'parent' element
512 el->isotopes[massNumberA] = isotope;
513}
514
516{
517 DEBUG_ENTRY("MeanMassOfElement()");
518
519 // if no isotopes have been defined, just use the mean mass stored elsewhere
520 if( el->isotopes.size()==0 )
521 return dense.AtomicWeight[el->Z-1];
522
523 realnum averageMass = 0., fracsum = 0.;
524 for( isotopes_i it = el->isotopes.begin(); it != el->isotopes.end(); ++it )
525 {
526 fracsum += it->second->frac;
527 averageMass += it->second->mass_amu * it->second->frac;
528 }
529 ASSERT( fp_equal( fracsum, realnum(1.) ) );
530
531 return averageMass;
532}
533
535 const char label[],
536 enum spectype type,
537 enum mole_state state,
538 realnum form_enthalpy)
539{
540 return newspecies( label, type, state, form_enthalpy, true);
541}
542
543/* Parse species string to find constituent atoms, charge etc. */
545 const char label[],
546 enum spectype type,
547 enum mole_state state,
548 realnum form_enthalpy,
549 bool lgCreateIsotopologues )/* formation enthalpy at 0K */
550{
551 DEBUG_ENTRY("newspecies()");
552
553 int exists;
554 ChemAtomList atomsLeftToRight;
555 vector< int > numAtoms;
556 string embellishments;
557 char *s;
559
560 mol->parentLabel.clear();
561 mol->isEnabled = true;
562 mol->charge = 0;
563 mol->lgExcit = false;
564 mol->mole_mass = 0.0;
565 mol->state = state;
566 mol->lgGas_Phase = true;
567 mol->form_enthalpy = form_enthalpy;
568 mol->groupnum = -1;
569
570 /* Create private workspace for label; copy and strip trailing whitespace */
571 int len = strlen(label)+1;
572 auto_vec<char> mylab_v(new char[len]);
573 char *mylab = mylab_v.data();
574 strncpy(mylab,label,len);
575 s = strchr(mylab,' ');
576 if (s)
577 *s = '\0';
578 mol->label = mylab;
579
580 if(type == MOLECULE)
581 {
582 if( parse_species_label( mylab, atomsLeftToRight, numAtoms, embellishments, mol->lgExcit, mol->charge, mol->lgGas_Phase ) == false )
583 return NULL;
584
585 for( unsigned i = 0; i < atomsLeftToRight.size(); ++i )
586 {
587 mol->nAtom[ atomsLeftToRight[i] ] += numAtoms[i];
588 mol->mole_mass += numAtoms[i] * atomsLeftToRight[i]->mass_amu * (realnum)ATOMIC_MASS_UNIT;
589 }
590 }
591
592 // we also kill H- if molecules are disabled. This is less than ideal,
593 // but physically motivated by the fact that one of the strongest H- sinks
594 // involves formation of H2 (disabled by "no molecules"), while one the
595 // fastest sources is e- recombination (which would still be allowed).
596 if ( (mol->n_nuclei() > 1 || (mol->isMonatomic() && mol->charge==-1) ) && mole_global.lgNoMole)
597 {
598 if( trace.lgTraceMole )
599 fprintf(ioQQQ,"No species %s as molecules off\n",label);
600 return NULL;
601 }
602
603 if (mol->n_nuclei() > 1 && mole_global.lgNoHeavyMole)
604 {
605 for( nAtoms_ri it=mol->nAtom.rbegin(); it != mol->nAtom.rend(); --it )
606 {
607 if( it->first->el->Z-1 != ipHYDROGEN )
608 {
609 ASSERT( it->second > 0 );
610 if( trace.lgTraceMole )
611 fprintf(ioQQQ,"No species %s as heavy molecules off\n",label);
612 return NULL;
613 }
614 }
615 }
616
617 /* Insert species into hash table */
618 exists = (mole_priv::spectab.find(mol->label) != mole_priv::spectab.end());
619 if( exists )
620 {
621 fprintf( ioQQQ,"Warning: duplicate species %s - using first one\n",
622 mol->label.c_str() );
623 return NULL;
624 }
625 else
626 mole_priv::spectab[mol->label] = mol;
627
628 // all map entries should have strictly positive number of nuclei
629 for( nAtoms_i it=mol->nAtom.begin(); it != mol->nAtom.end(); ++it )
630 ASSERT( it->second > 0 );
631
632 if (state != MOLE_NULL)
633 mole_global.num_total++;
634 if (state == MOLE_ACTIVE)
635 mole_global.num_calc++;
636
637 // this is a special case to always treat 13CO (as has long been done)
638 if( lgCreateIsotopologues && type==MOLECULE && mol->label.compare("CO")==0 )
639 {
640 molecule *sp = newspecies( "^13CO", MOLECULE, mol->state, mol->form_enthalpy, false );
641 sp->parentLabel = mol->label;
642 }
643
644 // create singly-substituted isotopologues
645 if( lgCreateIsotopologues && type==MOLECULE && !mol->isMonatomic() )
646 {
647 for( nAtoms_i it1 = mol->nAtom.begin(); it1 != mol->nAtom.end(); ++it1 )
648 {
649 for( map<int, count_ptr<chem_atom> >::iterator it2 = it1->first->el->isotopes.begin();
650 it2 != it1->first->el->isotopes.end(); ++it2 )
651 {
652 // we don't want to create ^1H isotopologues (only substitute D for H)
653 if( it1->first->el->Z-1 == ipHYDROGEN && it2->second->A != 2 )
654 continue;
655 if( !mole_global.lgTreatIsotopes[it1->first->el->Z-1] )
656 continue;
657 if( it2->second->lgMeanAbundance() )
658 continue;
659 vector<string> isoLabs;
660 create_isotopologues_one( atomsLeftToRight, numAtoms, it1->first->label(), it2->second->label(), embellishments, isoLabs );
661 //fprintf( ioQQQ, " DEBUGGG %10s isotopologues of %10s:", it1->first->label().c_str(), mol->label.c_str() );
662 //for( vector<string>::iterator lab = isoLabs.begin(); lab != isoLabs.end(); ++ lab )
663 // fprintf( ioQQQ, " %10s", lab->c_str() );
664 //fprintf( ioQQQ, "\n" );
665 for( vector<string>::iterator newLabel = isoLabs.begin(); newLabel != isoLabs.end(); ++newLabel )
666 {
667 molecule *sp = newspecies( newLabel->c_str(), MOLECULE, mol->state, mol->form_enthalpy, false );
668 // D is special -- don't set parentLabel
669 if( sp!=NULL && it2->second->A != 2 )
670 sp->parentLabel = mol->label;
671 }
672 }
673 }
674 }
675
676 return &(*mol);
677}
678bool parse_species_label( const char label[], ChemAtomList &atomsLeftToRight, vector<int> &numAtoms, string &embellishments )
679{
680 bool lgExcit, lgGas_Phase;
681 int charge;
682 bool lgOK = parse_species_label( label, atomsLeftToRight, numAtoms, embellishments, lgExcit, charge, lgGas_Phase );
683 return lgOK;
684}
685bool parse_species_label( const char label[], ChemAtomList &atomsLeftToRight, vector<int> &numAtoms, string &embellishments,
686 bool &lgExcit, int &charge, bool &lgGas_Phase )
687{
688 long int i, n, ipAtom;
689 char thisAtom[CHARS_ISOTOPE_SYM];
691 char mylab[CHARS_SPECIES];
692 char *s;
693
694 strncpy( mylab, label, CHARS_SPECIES );
695
696 /* Excitation... */
697 s = strpbrk(mylab,"*");
698 if(s)
699 {
700 lgExcit = true;
701 embellishments = s;
702 *s = '\0';
703 }
704
705 /* ...Charge */
706 s = strpbrk(mylab,"+-");
707 if(s)
708 {
709 if(isdigit(*(s+1)))
710 n = atoi(s+1);
711 else
712 n = 1;
713 if(*s == '+')
714 charge = n;
715 else
716 charge = -n;
717 embellishments = s + embellishments;
718 *s = '\0';
719 }
720 /* ...Grain */
721 s = strstr(mylab,"grn");
722 if(s)
723 {
724 lgGas_Phase = false;
725 embellishments = s + embellishments;
726 *s = '\0';
727 }
728 else
729 {
730 lgGas_Phase = true;
731 }
732 //fprintf( ioQQQ, "DEBUGGG parse_species_label %s\t%s\t%s\n", label, mylab, embellishments.c_str() );
733
734 /* Now analyse chemical formula */
735 i = 0;
736 while (mylab[i] != '\0' && mylab[i] != ' ' && mylab[i] != '*')
737 {
738 /* Select next element in species, matches regexp [A-Z][a-z]? */
739 ipAtom = 0;
740 /* look for isotope prefix */
741 if(mylab[i]=='^')
742 {
743 thisAtom[ipAtom++] = mylab[i++];
744 ASSERT( isdigit(mylab[i]) );
745 thisAtom[ipAtom++] = mylab[i++];
746 if(isdigit(mylab[i]))
747 {
748 thisAtom[ipAtom++] = mylab[i++];
749 }
750 }
751 // should be first character of an element symbol
752 thisAtom[ipAtom++] = mylab[i++];
753 if(islower(mylab[i]))
754 {
755 thisAtom[ipAtom++] = mylab[i++];
756 }
757 thisAtom[ipAtom] = '\0';
758 ASSERT(ipAtom <= CHARS_ISOTOPE_SYM);
759
760 atom = findatom(thisAtom);
761 if(atom.get_ptr() == NULL)
762 {
763 fprintf(stderr,"Did not recognize atom at %s in \"%s \"[%ld]\n",thisAtom,mylab,i);
764 exit(-1);
765 }
766 if(!dense.lgElmtOn[atom->el->Z-1])
767 {
768 if( trace.lgTraceMole )
769 fprintf(ioQQQ,"No species %s as element %s off\n",mylab,atom->el->label.c_str() );
770 return false;
771 }
772
773 if(isdigit(mylab[i])) /* If there is >1 of this atom */
774 {
775 n = 0;
776 do {
777 n = 10*n+(long int)(mylab[i]-'0');
778 i++;
779 // the test of i prevents a warning by g++ 4.8.0: array subscript is above array bounds
780 } while (i < CHARS_SPECIES && isdigit(mylab[i]));
781 }
782 else
783 {
784 n = 1;
785 }
786 atomsLeftToRight.push_back( atom );
787 numAtoms.push_back( n );
788 }
789
790 return true;
791}
792STATIC bool isactive(const molecule &mol)
793{
794 DEBUG_ENTRY("isactive()");
795 return mol.state == MOLE_ACTIVE;
796}
797STATIC bool ispassive(const molecule &mol)
798{
799
800 DEBUG_ENTRY("ispassive()");
801 return mol.state == MOLE_PASSIVE;
802}
803
804bool lgDifferByExcitation( const molecule &mol1, const molecule &mol2 )
805{
806 if( mol1.label == mol2.label + "*" )
807 return true;
808 else if( mol2.label == mol1.label + "*" )
809 return true;
810 else
811 return false;
812}
813
814molecule *findspecies(const char buf[])
815{
816 DEBUG_ENTRY("findspecies()");
817
818 // strip string of the first space and anything after it
819 string s;
820 for (const char *pb = buf; *pb && *pb != ' '; ++pb)
821 {
822 s += *pb;
823 }
824
825 const molecule_i p = mole_priv::spectab.find(s);
826
827 if(p != mole_priv::spectab.end())
828 return &(*p->second);
829 else
830 return null_mole;
831}
832
833molezone *findspecieslocal(const char buf[])
834{
835 DEBUG_ENTRY("findspecieslocal()");
836
837 // strip string of the first space and anything after it
838 string s;
839 for (const char *pb = buf; *pb && *pb != ' '; ++pb)
840 {
841 s += *pb;
842 }
843
844 const molecule_i p = mole_priv::spectab.find(s);
845
846 if(p != mole_priv::spectab.end())
847 return &mole.species[ p->second->index ];
848 else
849 return null_molezone;
850}
851
853{
854 chem_atom_i p;
855
856 DEBUG_ENTRY("findatom()");
857
858 p = mole_priv::atomtab.find(buf);
859
860 if(p != mole_priv::atomtab.end())
861 return p->second;
862 else
863 return count_ptr<chem_atom>(NULL);
864}
865
867{
868 int i;
869 double den_times_area, den_grains, adsorbed_density;
870
871 DEBUG_ENTRY("mole_update_species_cache()");
872
873 enum { DEBUG_MOLE = false };
874
875 /* For rates that are dependent on grain physics. This includes grain density,
876 cross sectional area, and dust temperature of each constituent. Note that
877
878 gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3
879
880 is the integrated projected grain surface area per cm^3 of gas for each grain size bin */
881
882 /* >>chng 06 feb 28, turn off this rate when no grain molecules */
883 /* >>chng 06 dec 05 rjrw: do this in newreact rather than rate */
884 if( gv.lgDustOn() )
885 {
886 den_times_area = den_grains = 0.0;
887 for( size_t nd=0; nd < gv.bin.size(); nd++ )
888 {
889 /* >>chng 06 mar 04, update expression for projected grain surface area, PvH */
890 den_times_area += gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
891 den_grains += gv.bin[nd]->cnv_GR_pCM3;
892 }
893
894 adsorbed_density = 0.0;
895 for (i=0;i<mole_global.num_total;i++)
896 {
897 if( !mole_global.list[i]->lgGas_Phase && mole_global.list[i]->parentLabel.empty() )
898 adsorbed_density += mole.species[i].den;
899 }
900
901 mole.grain_area = den_times_area;
902 mole.grain_density = den_grains;
903
904 double mole_cs = 1e-15;
905 if (4*den_times_area <= mole_cs*adsorbed_density)
906 mole.grain_saturation = 1.0;
907 else
908 mole.grain_saturation = ((mole_cs*adsorbed_density)/(4.*den_times_area));
909 }
910 else
911 {
912 mole.grain_area = 0.0;
913 mole.grain_density = 0.0;
914 mole.grain_saturation = 1.0;
915 }
916
917 if (DEBUG_MOLE)
918 fprintf(ioQQQ,"Dust: %f %f %f\n",
919 mole.grain_area,mole.grain_density,mole.grain_saturation);
920
921 for (i=0;i<mole_global.num_total;i++)
922 {
923 if(mole.species[i].location != NULL)
924 {
925 ASSERT( mole_global.list[i]->parentLabel.empty() );
926 mole.species[i].den = *(mole.species[i].location);
927 if (DEBUG_MOLE)
928 fprintf(ioQQQ,"Updating %s: %15.8g\n",mole_global.list[i]->label.c_str(),mole.species[i].den);
929 }
930 }
931
932 mole.set_isotope_abundances();
933}
934
936// Finding the total atom density from MoleMap.molElems w
937{
938 int i;
939
940 /* These two updates should together maintain the abundance invariant */
941
942 // Assert invariant
944
945 // Update total of non-ladder species
946 dense.updateXMolecules();
947 total_molecule_deut( deut.xMolecules );
948
949 /* charge on molecules */
950 mole.elec = 0.;
951 for(i=0;i<mole_global.num_calc;i++)
952 {
953 if (mole.species[i].location == NULL && mole_global.list[i]->parentLabel.empty())
954 mole.elec += mole.species[i].den*mole_global.list[i]->charge;
955 }
956
957 // Update ionization ladder species
958 realnum delta = 0.0;
959 long ncpt = 0;
960
961 for (i=0;i<mole_global.num_total;i++)
962 {
963 if(mole.species[i].location && mole_global.list[i]->state == MOLE_ACTIVE)
964 {
965 realnum new_pop = mole.species[i].den;
966
967 if( mole_global.list[i]->isMonatomic() )
968 {
969 realnum old_pop = *(mole.species[i].location);
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]);
973 delta += frac*frac;
974 ++ncpt;
975 }
976
977 //if( iteration >= 3 && nzone >= 100 )
978 // fprintf( ioQQQ, "DEBUGGG mole_return_ %i\t%s\t%.12e\t%.12e\t%.12e\t%.12e\t%li\n",
979 // i, mole_global.list[i]->label.c_str(), new_pop, old_pop, frac, delta, ncpt );
980 *(mole.species[i].location) = new_pop;
981 }
982 }
983
984 // Assert invariant
986 return ncpt>0 ? sqrt(delta/ncpt) : 0.f;
987}
988
990{
991 DEBUG_ENTRY("t_mole_local::set_isotope_abundances()");
992
993 // loop over unresolved elements
994 for(ChemAtomList::iterator atom = unresolved_atom_list.begin(); atom != unresolved_atom_list.end(); ++atom)
995 {
996 // loop over all isotopes of each element
997 for( isotopes_i it = (*atom)->el->isotopes.begin(); it != (*atom)->el->isotopes.end(); ++it )
998 {
999 // skip mean-abundance "isotopes"
1000 if( it->second->lgMeanAbundance() )
1001 continue;
1002
1003 // loop over all ions of the isotope
1004 for( unsigned long ion=0; ion<it->second->ipMl.size(); ++ion )
1005 {
1006 if( it->second->ipMl[ion] != -1 &&
1007 (species[ it->second->ipMl[ion] ].location == NULL ) && (*atom)->ipMl[ion] != -1 )
1008 {
1009 species[ it->second->ipMl[ion] ].den = species[ (*atom)->ipMl[ion] ].den * it->second->frac;
1010 }
1011 }
1012 }
1013 }
1014
1015 return;
1016}
1017
1018void t_mole_local::set_location( long nelem, long ion, double *density )
1019{
1020 DEBUG_ENTRY( "t_mole_local::set_location()" );
1021
1022 ASSERT( nelem < LIMELM );
1023 ASSERT( ion < nelem + 2 );
1024 long mole_index = unresolved_atom_list[nelem]->ipMl[ion];
1025 // element not enabled if index is -1
1026 if( mole_index == -1 )
1027 return;
1028 ASSERT( mole_index < mole_global.num_total );
1029 species[mole_index].location = density;
1030
1031 return;
1032}
1033
1035{
1036 DEBUG_ENTRY( "total_molecule_deut()" );
1037
1038 double total = 0.;
1039
1040 if( !deut.lgElmtOn )
1041 return;
1042
1043 for (long int i=0;i<mole_global.num_calc;++i)
1044 {
1045 if (mole.species[i].location == NULL && mole_global.list[i]->parentLabel.empty() )
1046 {
1047 for( molecule::nAtomsMap::iterator atom=mole_global.list[i]->nAtom.begin();
1048 atom != mole_global.list[i]->nAtom.end(); ++atom)
1049 {
1050 long int nelem = atom->first->el->Z-1;
1051 if( nelem==0 && atom->first->A==2 )
1052 {
1053 total += mole.species[i].den*atom->second;
1054 }
1055 }
1056 }
1057 }
1058
1059 total_f = (realnum)total;
1060
1061 return;
1062}
1064{
1065 DEBUG_ENTRY( "total_molecule_elems()" );
1066
1067 /* now set total density of each element locked in molecular species */
1068 for ( long int nelem=ipHYDROGEN;nelem<LIMELM; ++nelem )
1069 {
1070 total[nelem] = 0.;
1071 }
1072 for (long int i=0;i<mole_global.num_calc;++i)
1073 {
1074 if (mole.species[i].location == NULL && mole_global.list[i]->parentLabel.empty() )
1075 {
1076 for( molecule::nAtomsMap::iterator atom=mole_global.list[i]->nAtom.begin();
1077 atom != mole_global.list[i]->nAtom.end(); ++atom)
1078 {
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;
1083 }
1084 }
1085 }
1086}
1087void total_network_elems(double total[LIMELM])
1088{
1089 DEBUG_ENTRY( "total_network_elems()" );
1090
1091 /* now set total density of each element locked in molecular species */
1092 for ( long int nelem=ipHYDROGEN;nelem<LIMELM; ++nelem )
1093 {
1094 total[nelem] = 0.;
1095 }
1096 for (long int i=0;i<mole_global.num_calc;++i)
1097 {
1098 if (mole_global.list[i]->parentLabel.empty())
1099 {
1100 for( molecule::nAtomsMap::iterator atom=mole_global.list[i]->nAtom.begin();
1101 atom != mole_global.list[i]->nAtom.end(); ++atom)
1102 {
1103 long int nelem = atom->first->el->Z-1;
1104 total[nelem] += (realnum) mole.species[i].den*atom->second;
1105 }
1106 }
1107 }
1108}
1110{
1111 long int i;
1112 realnum total;
1113
1114 DEBUG_ENTRY( "total_molecules()" );
1115
1116 total = 0.;
1117 for (i=0;i<mole_global.num_calc;++i)
1118 {
1119 if (mole.species[i].location == NULL && mole_global.list[i]->parentLabel.empty())
1120 total += (realnum) mole.species[i].den;
1121 }
1122 return total;
1123}
1125{
1126 long int i;
1127 realnum total;
1128
1129 DEBUG_ENTRY( "total_molecules_gasphase()" );
1130
1131 total = 0.;
1132 for (i=0;i<mole_global.num_calc;++i)
1133 {
1134 if (mole_global.list[i]->lgGas_Phase && mole.species[i].location == NULL && mole_global.list[i]->parentLabel.empty())
1135 total += (realnum) mole.species[i].den;
1136 }
1137 return total;
1138}
1139/*lint +e778 const express eval to 0 */
1140/*lint +e725 expect positive indentation */
1141
1143{
1144 long int i, j;
1145 /* Neutrals and positive ions will be treated as single species inside
1146 molecular equilibrium solver, to facilitate coupling with ionization
1147 solvers */
1148 DEBUG_ENTRY ("mole_make_groups()");
1149 if (mole_global.num_calc == 0)
1150 {
1151 groupspecies = NULL;
1152 mole_global.num_compacted = 0;
1153 return;
1154 }
1155 groupspecies = (molecule **) MALLOC(mole_global.num_calc*sizeof(molecule *));
1156 for (i=0,j=0;i<mole_global.num_calc;i++)
1157 {
1158 if( mole_global.list[i]->parentLabel.empty() && ( !mole_global.list[i]->isMonatomic() || mole_global.list[i]->charge <= 0 || ! mole_global.list[i]->lgGas_Phase ) )
1159 {
1160 /* Compound molecules and negative ions are represented individually */
1161 mole_global.list[i]->groupnum = j;
1162 groupspecies[j++] = &(*mole_global.list[i]);
1163 }
1164 else
1165 {
1166 /* All positive ions are collapsed into single macrostate (i.e. H+ -> H) */
1167 /* Need to increase constant if higher ions are included */
1168 ASSERT (mole_global.list[i]->charge < LIMELM+1);
1169 ASSERT (mole_global.list[i]->groupnum == -1);
1170 }
1171 }
1172 mole_global.num_compacted = j;
1174 mole_global.num_compacted*sizeof(molecule *));
1175
1176 for (i=0;i<mole_global.num_calc;i++)
1177 {
1178 if (mole_global.list[i]->groupnum == -1)
1179 {
1180 if( mole_global.list[i]->isMonatomic() && mole_global.list[i]->parentLabel.empty() )
1181 {
1182 for( nAtoms_i it = mole_global.list[i]->nAtom.begin(); it != mole_global.list[i]->nAtom.end(); ++it )
1183 {
1184 ASSERT( it->second > 0 );
1185 mole_global.list[i]->groupnum = mole_global.list[it->first->ipMl[0]]->groupnum;
1186 break;
1187 }
1188 }
1189 else
1190 {
1191 ASSERT( !mole_global.list[i]->parentLabel.empty() );
1192 mole_global.list[i]->groupnum = mole_global.list[ mole_global.list[i]->parentIndex ]->groupnum;
1193 }
1194 }
1195
1196 ASSERT( mole_global.list[i]->groupnum != -1);
1197 }
1198
1199 return;
1200}
1201
t_atoms atoms
Definition atoms.cpp:5
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
@ CHARS_ISOTOPE_SYM
Definition cddefines.h:275
#define STATIC
Definition cddefines.h:97
#define REALLOC
Definition cddefines.h:519
const int ipCARBON
Definition cddefines.h:310
@ CHARS_SPECIES
Definition cddefines.h:274
#define MALLOC(exp)
Definition cddefines.h:501
struct t_species species
Definition cddefines.h:1224
const int LIMELM
Definition cddefines.h:258
const int INPUT_LINE_LENGTH
Definition cddefines.h:254
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition service.cpp:381
float realnum
Definition cddefines.h:103
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void fixit(void)
Definition service.cpp:991
element_type * data() const
Definition cddefines.h:1175
T * get_ptr() const
Definition count_ptr.h:49
bool lgGas_Phase
Definition mole.h:146
int charge
Definition mole.h:144
string parentLabel
Definition mole.h:137
nAtomsMap nAtom
Definition mole.h:143
realnum mole_mass
Definition mole.h:165
string label
Definition mole.h:142
enum mole_state state
Definition mole.h:168
bool isMonatomic(void) const
Definition mole.h:157
int index
Definition mole.h:169
iterator begin(size_type i1)
void make_species(void)
static void sort(MoleculeList::iterator start, MoleculeList::iterator end)
vector< bool > lgTreatIsotopes
Definition mole.h:311
void set_location(long nelem, long ion, double *dense)
void set_isotope_abundances(void)
valarray< class molezone > species
Definition mole.h:398
t_co co
Definition co.cpp:5
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
const ios_base::openmode mode_r
Definition cpu.h:212
bool lgElemsConserved(void)
Definition dense.cpp:99
t_dense dense
Definition dense.cpp:24
void SetGasPhaseDeuterium(const realnum &Hdensity)
Definition deuterium.cpp:69
void InitDeuteriumIonization()
Definition deuterium.cpp:19
void SetDeuteriumFractionation(const realnum &frac)
Definition deuterium.cpp:61
t_deuterium deut
Definition deuterium.cpp:8
t_elementnames elementnames
@ CHARS_ION_STAGE
Definition elementnames.h:8
GrainVar gv
Definition grainvar.cpp:5
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
t_hmi hmi
Definition hmi.cpp:5
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
chem_element * null_element
ChemAtomList atom_list
t_mole_global mole_global
Definition mole.cpp:6
chem_atom * null_atom
mole_state
Definition mole.h:15
@ MOLE_NULL
Definition mole.h:15
@ MOLE_ACTIVE
Definition mole.h:15
@ MOLE_PASSIVE
Definition mole.h:15
ChemAtomList unresolved_atom_list
molecule::nAtomsMap::iterator nAtoms_i
Definition mole.h:214
molecule * findspecies(const char buf[])
molezone * null_molezone
map< int, count_ptr< chem_atom > >::iterator isotopes_i
Definition mole.h:35
molecule * null_mole
vector< count_ptr< chem_atom > > ChemAtomList
Definition mole.h:115
map< string, count_ptr< chem_atom > >::iterator chem_atom_i
Definition mole_priv.h:41
molecule ** groupspecies
map< string, count_ptr< molecule > >::iterator molecule_i
Definition mole_priv.h:39
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)
void mole_make_list()
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])
spectype
@ MOLECULE
@ OTHER
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
Definition physconst.h:91
UNUSED const double KJMOL1CM
Definition physconst.h:170
UNUSED const double ATOMIC_MASS_UNIT
Definition physconst.h:88
t_state state
Definition state.cpp:19
t_trace trace
Definition trace.cpp:5