mole.h

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */

#ifndef MOLE_H_
#define MOLE_H_

/* mole.h */

#include "module.h"
#include "container_classes.h"

class TransitionList;
class qList;
class species;
class molezone;

static const double SMALLABUND = 1e-24;

enum mole_state {MOLE_NULL, MOLE_PASSIVE, MOLE_ACTIVE};

class chem_nuclide;

class chem_element {
        explicit chem_element(); // Do not implement
        chem_element &operator=(const chem_element&); // Do not implement
public:
        explicit chem_element(int Z, const char*label) : Z(Z), label(label)
        {}
        ~chem_element() throw()
        {}
        const int Z;
        const string label;
        map<int, shared_ptr<chem_nuclide> > isotopes;
        //(first -> Atomic A; second -> chem_nuclide )
        //(first -> -1 for bogus isotope, i.e. where no
        //      isotopes have been explicitly defined)
};

typedef map<int, shared_ptr<chem_nuclide> >::iterator isotopes_i;

class chem_nuclide {
public:
        // Link back to basic element for convenience -- not a shared_ptr
        // as this would lead to a reference cycle (and hence the
        // destructors never being called).  Many-to-one relation suggests
        // that the weak link should be this way around.
        chem_element* wel;
        int A;               /* mass number */
        unsigned long index; /* index in list for temporary use */
        vector<int> ipMl;    /* Atom and ion species in molecule arrays */
        realnum mass_amu;    /* mass of isotope in AMU */
        double frac;         /* fraction of element in this isotope */

        chem_element* el() const
        {
                return wel;
        }
        bool lgMeanAbundance( void ) const
        {
                return ( A < 0 );
        }

        // May have independent logic in the future...
        bool lgHasLinkedIon( void ) const
        {
                return lgMeanAbundance();
        }

        /* Chemical symbols for elements */
        string label( void ) const
        {
                if( lgMeanAbundance() )
                        return el()->label;
                else if( el()->Z==1 && A==2 )
                {
                        // Deuterium is a special case
                        return "D";
                }
                else
                {
                        char str[20];
                        sprintf(str,"^%d",A);
                        return ( str + el()->label );
                }
        }
        
        int compare ( const chem_nuclide &b ) const
        {
                // sort by proton number first
                if ( el()->Z < b.el()->Z )
                        return -1;
                if ( el()->Z > b.el()->Z )
                        return 1;

                if ( A < b.A )
                        return -1;
                if ( A > b.A )
                        return 1;

                return 0;
        }
};
inline bool operator< (const chem_nuclide &a, const chem_nuclide &b)
{
        return a.compare(b) < 0;
}
inline bool operator> (const chem_nuclide &a, const chem_nuclide &b)
{
        return a.compare(b) > 0;
}
inline bool operator<= (const chem_nuclide &a, const chem_nuclide &b)
{
        return a.compare(b) <= 0;
}
inline bool operator>= (const chem_nuclide &a, const chem_nuclide &b)
{
        return a.compare(b) >= 0;
}
inline bool operator== (const chem_nuclide &a, const chem_nuclide &b)
{
        return a.compare(b) == 0;
}
inline bool operator!= (const chem_nuclide &a, const chem_nuclide &b)
{
        return !(a == b);
}

typedef vector< shared_ptr<chem_nuclide> > ChemNuclideList;
extern ChemNuclideList nuclide_list;
typedef vector< shared_ptr<chem_element> > ChemElementList;
extern ChemElementList element_list;
extern chem_nuclide *null_nuclide;

class element_pointer_value_less
{
public:
        bool operator()(const shared_ptr<chem_nuclide>& a,
                                                const shared_ptr<chem_nuclide>& b) const
        {
                return *a < *b;
        }
};

/* Structure containing molecule data, initially only CO */
class molecule {
public:
        typedef map<const shared_ptr<chem_nuclide>, int,
                element_pointer_value_less> nNucsMap;

        string parentLabel;
        bool isIsotopicTotalSpecies() const
        {
                return parentLabel.empty();
        }
        int parentIndex;
        bool isEnabled;   /* Is it enabled? */

        /* Species physical data */
        string    label;        
        nNucsMap nNuclide;      
        int     charge;         
        bool    lgExcit;        
        bool    lgGas_Phase;    
        int     n_nuclei(void) const       
        {
                int num = 0;
                for (nNucsMap::const_iterator el = nNuclide.begin(); 
                          el != nNuclide.end(); ++el)
                {
                        num += el->second;
                }
                return num;
        }
        int nElement(int nelem)
        {
                int num = 0;
                for (nNucsMap::const_iterator el = nNuclide.begin(); 
                          el != nNuclide.end(); ++el)
                {
                        if (el->first->el()->Z-1 == nelem)
                                num += el->second;
                }
                return num;
        }

        bool isMonatomic(void) const       
        {
                if (nNuclide.size() == 1 && nNuclide.begin()->second == 1)
                        return true;
                return false;
        }

        bool isMolecule() const
        {
                if( nNuclide.size() > 1 || ( nNuclide.size() == 1 && nNuclide.begin()->second > 1 ) )
                        return true;
                return false;
        }

        realnum form_enthalpy;  
        realnum mole_mass;      
        /* Parameters as computational object */
        enum mole_state state;
        int index, groupnum;
        int n_react;            
        class molezone *local(void) const;

        chem_nuclide *heavyAtom(void) //const
        {
                for( nNucsMap::reverse_iterator it=nNuclide.rbegin(); it!=nNuclide.rend(); ++it )
                {
                        if (0 != it->second )
                        {
                                return it->first.get();
                        }
                }
                return null_nuclide;
        }

        int compare(const molecule &mol2) const
        {
                nNucsMap::const_reverse_iterator it1, it2;

                for( it1 = nNuclide.rbegin(), it2 = mol2.nNuclide.rbegin();
                        it1 != nNuclide.rend() && it2 != mol2.nNuclide.rend(); ++it1, ++it2 )
                {
                        if( *(it1->first) > *(it2->first) )
                                return 1;
                        else if( *(it1->first) < *(it2->first) )
                                return -1;
                        else if( it1->second > it2->second)
                                return 1;
                        else if( it1->second < it2->second)
                                return -1;
                }

                if( it1 != nNuclide.rend() && it2 == mol2.nNuclide.rend() )
                        return 1;
                else if( it1 == nNuclide.rend() && it2 != mol2.nNuclide.rend() )
                        return -1;
                else
                        ASSERT( it1 == nNuclide.rend() && it2 == mol2.nNuclide.rend() );

                // sort by label if falls through to here       
                return ( label.compare(mol2.label) );
                        
        }
}; 

/* iterators on nNuclide */     
typedef molecule::nNucsMap::iterator nNucs_i;
typedef molecule::nNucsMap::reverse_iterator nNucs_ri;
typedef molecule::nNucsMap::const_reverse_iterator nNucs_cri;

extern void mole_drive(void);

extern void mole_create_react(void);

class mole_reaction;

mole_reaction *mole_findrate_s(const char buf[]);

extern molecule *null_mole;

inline bool exists (const molecule* m)
{
        return m != null_mole;
}

extern molecule *findspecies(const char buf[]);
extern molecule *findspecies_validate(const char buf[]);

extern molezone *findspecieslocal(const char buf[]);
extern molezone *findspecieslocal_validate(const char buf[]);

extern shared_ptr<chem_nuclide> findnuclide(const char buf[]);

class Properties;


class t_mole_global : public module {

public:
        const char *chName() const
        {
                return "mole_global";
        }
        void zero();
        void comment(t_warnings&) {}
        void init(void);

        void make_species(void);
                
        bool lgNoMole;

        bool lgNoHeavyMole;

        bool lgH2Ozer;

        bool lgLeidenHack;

        bool lgFederman;

        bool lgStancil;

        bool lgNonEquilChem;

        bool lgProtElim;

        bool lgNeutrals;

        bool lgGrain_mole_deplete;

        // flag saying whether to model isotopes (and isotopologues) of a given element
        vector<bool> lgTreatIsotopes;

        int num_total, num_calc, num_compacted;

        typedef vector<shared_ptr<molecule> > MoleculeList;
        MoleculeList list;

        map<string,bool> offReactions;
        map<string,Properties > speciesProperties;

        static void sort(MoleculeList::iterator start,
                                 MoleculeList::iterator end);
        static void sort(molecule **start, molecule **end);

};

extern t_mole_global mole_global;

class molezone {
public:
        molezone()
        {
                init();
        }
        void init (void)
        {
                location = NULL;
                levels = NULL;
                lines = NULL;
                dbase = NULL;
                zero();
        }
        void zero (void)
        {
                src = 0.;
                snk = 0.;
                den = 0.;
                column = 0.;
                atomLim = null_nuclide;
                xFracLim = 0.;
                column_old = 0.;
                index = -1;
        }
        int index;
        const molecule* global() const
        {
                if (index == -1)
                        return null_mole;
                return mole_global.list[index].get();
        }
        
        double *location;      
        double src, snk;

        species* dbase;
        qList* levels;
        TransitionList* lines;

        /* Current zone data */
        double  den;       
        realnum column;    
        chem_nuclide* atomLim; /* atom which is closest to limiting molecule density */
        realnum xFracLim;  
        /* Historical solution data */
        realnum column_old;   
};

extern molezone *null_molezone;

class t_mole_local : public module
{
public:
        void alloc();
        void zero();
        void comment(t_warnings&) {}
        const char *chName() const
        {
                return "mole";
        }
        void set_ion_locations( );
        void set_isotope_abundances( void );
        double sink_rate_tot(const char chSpecies[]) const;
        double sink_rate_tot(const molecule* const sp) const;
        double sink_rate(const molecule* const sp, const mole_reaction& rate) const;
        double sink_rate(const molecule* const sp, const char buf[]) const;
        double source_rate_tot(const char chSpecies[]) const;
        double source_rate_tot(const molecule* const sp) const;
        double dissoc_rate(const char chSpecies[]) const;
        double chem_heat(void) const;
        double findrk(const char buf[]) const;
        double findrate(const char buf[]) const;

        double grain_area, grain_density, grain_saturation;

        double elec;

        multi_arr<double,2> source, sink;
        
        multi_arr<realnum,3> xMoleChTrRate;

        valarray<class molezone> species;

        vector<double> reaction_rks;
        vector<double> old_reaction_rks;
        long old_zone;
};

extern t_mole_local mole;

inline bool exists (const molezone* m)
{
        return m != null_molezone;
}

inline molezone *molecule::local(void) const
{
        if ( index != -1)
                return &mole.species[index];
        else
                return null_molezone;
}

extern void total_molecule_elems(realnum total[LIMELM]);
extern void total_molecule_deut(realnum &total);

extern realnum total_molecules(void);

extern realnum total_molecules_gasphase(void);

realnum species_gasphase_density( const string &chSpecies );

extern void mole_make_list(void);
extern void mole_make_groups(void);

void mole_cmp_num_in_out_reactions(void);

bool lgDifferByExcitation( const molecule &mol1, const molecule &mol2 );

extern void mole_update_species_cache(void);

void mole_update_sources(void);

void mole_rk_bigchange(void);

void create_isotopologues(
        ChemNuclideList& atoms,
        vector< int >& numAtoms,
        string atom_old,
        string atom_new,
        string embellishments,
        vector<string>& newLabels );

void create_isotopologues_one_position(
        unsigned position,
        ChemNuclideList& atoms,
        vector< int >& numAtoms,
        string atom_old,
        string atom_new,
        string embellishments,
        string& newLabel );

bool parse_species_label( const char label[], ChemNuclideList &atomsLeftToRight, vector<int> &numAtoms, string &embellishments );
bool parse_species_label( const char mylab[], ChemNuclideList &atomsLeftToRight, vector<int> &numAtoms, string &embellishments,
        bool &lgExcit, int &charge, bool &lgGas_Phase );

/*HMRATE compile molecular rates using Hollenbach and McKee fits */
/* #define HMRATE(a,b,c) ( ((b) == 0 && (c) == 0) ? (a) : \
 *      ( ((c) == 0) ? (a)*pow(phycon.te/300.,(b)) : \
 *      ( ((c)/phycon.te > 50.) ? 0. : ( ((b) == 0) ?  (a)*exp(-(c)/phycon.te) : \
 *                                       (a)*pow(phycon.te/300.,(b))*exp(-(c)/phycon.te) ) ) ) ) */

inline double hmrate4( double a, double b, double c, double te )
{
        if( b == 0. && c == 0. )
                return a;
        else if( c == 0. )
                return a*pow(te/300.,b);
        else if( b == 0. )
                return a*sexp(c/te);
        else
                return a*pow(te/300.,b)*sexp(c/te);
}

class Parser;
void ParseChemistry(Parser &p);

bool isSpecies( const string &chSpecies );

bool isMolecule( const string &chSpecies );

void test_isMolecule();

void getMolecules( vector<string> &allMolecules );

#endif /* MOLE_H_ */