parse_atom_h2.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*ParseAtomH2 parse information from the atom command line */
#include "cddefines.h" 
#include "hmi.h" 
#include "h2.h" 
#include "h2_priv.h" 
#include "parser.h" 
#include "thirdparty.h"
#include "deuterium.h"

/*ParseAtomH2 parse information from the atom command line */
void ParseAtomH2(Parser &p )
{
        long int j;

        DEBUG_ENTRY( "ParseAtomH2()" );

        fixit(); // this must be generalized!!!
                // easiest way is to create and populate diatoms before parsing,
                // so we can simply iterator of diatoms and match strings.
                // probably will want to put label on command line in quotes
                // and change command to, for example,
                // diatom "H2"
                // do that work before entering this routine.
        diatomics *diatom = NULL;

        if( p.nMatch(" H2 " ) )
        {
                diatom = &h2;
                /* this command has a 2 in the H2 label - must not parse the two by
                 * accident.  Get the first number off the line image, and confirm that
                 * it is a 2 */
                j = (long int)p.FFmtRead();
                if( j != 2 )
                {
                        fprintf( ioQQQ, " Something is wrong with the order of the numbers on this line.\n" );
                        fprintf( ioQQQ, " The first number I encounter should be a 2.\n Sorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }
        else if( p.nMatch(" HD " ) )
        {
                diatom = &hd;
                // turn deuterium on
                deut.lgElmtOn = true;
        }

        /* the mere calling of this routine turns the large H2 molecule on */
        diatom->lgEnabled = true;

        if( p.nMatch("LEVE") )
        {
                /* number of electronic levels */

                /* lgREAD_DATA is false at start of calculation, set true when 
                 * space allocated for the H lines.  Once done we must ignore all 
                 * future changes in the number of levels */
                if( !diatom->lgREAD_DATA )
                {
                        diatom->n_elec_states = (long int)p.FFmtRead();
                        if( p.lgEOL() )
                        {
                                if( p.nMatch("LARG") )
                                {
                                        /* LARGE is option to use the most number of electronic levels */
                                        diatom->n_elec_states = N_ELEC;
                                }
                                else
                                {
                                        p.NoNumb("number of electronic levels");
                                }
                        }

                        /* do not allow fewer than 3 - that includes Lyman & Werner bands */
                        if( diatom->n_elec_states < 3 )
                        {
                                fprintf( ioQQQ, " This would be too few electronic levels - resetting to 3.\n" );
                                diatom->n_elec_states = 3;
                        }
                        /* N_ELEC is the greatest number of elec lev possible */
                        else if( diatom->n_elec_states > N_ELEC )
                        {
                                fprintf( ioQQQ, 
                                        " This would be too many levels, the limit is %i.\n" , 
                                        N_ELEC);
                                cdEXIT(EXIT_FAILURE);
                        }
                }
        }

        else if( p.nMatch("LIMI") )
        {
                /* the limit to the H2 / Htot ratio - 
                 * if smaller than this, do not compute large H2 mole */
                diatom->H2_to_H_limit = p.FFmtRead();
                if( p.lgEOL() )
                {
                        /* did not find a number, either mistake or key " off" */
                        if( p.nMatch( " OFF"  ) )
                        {
                                /* turn off limit */
                                diatom->H2_to_H_limit = -1.;
                        }
                        else
                        {
                                p.NoNumb( "limit to the H2 / Htot ratio" );
                        }
                }
                else
                {
                        /* got a number, check if negative and so a log */
                        /* a number <= 0 is the log of the ratio */
                        if( diatom->H2_to_H_limit <= 0. )
                                diatom->H2_to_H_limit = pow(10., diatom->H2_to_H_limit);
                }
        }
        else if( p.nMatch("GBAR" ) )
        {
                /* option to either use, or not use, gbar approximation for low X 
                 * levels with no collision data - by default it is on */
                if( p.nMatch(" OFF" ) )
                {
                        diatom->lgColl_gbar = false;
                }
                else if( p.nMatch(" ON " ) )
                {
                        diatom->lgColl_gbar = true;
                }
                else
                {
                        fprintf( ioQQQ, 
                                " The gbar approximation must be off (\" OFF\") or on (\" ON \").\n");
                        cdEXIT(EXIT_FAILURE);
                }
        }
        /* option to turn collisional effects off or on */
        else if( p.nMatch("COLL" ) )
        {
                /* option to turn collisional dissociation off or on */
                if( p.nMatch("DISS" ) )
                {
                        /* option to turn collisions off */
                        if( p.nMatch(" ON " ) )
                        {
                                /* this is the default, leave collisions off */
                                diatom->lgColl_dissoc_coll = true;
                        }
                        else
                        {
                                /* default (and only reason for this command) is to turn off collisions */
                                diatom->lgColl_dissoc_coll = false;
                        }
                }
                /* option to turn collisional dissociation off or on 
                 * >>chng 06 mar 01, had been simply if - so all collisions were turned off
                 * when dissociation collisions turned off - 
                 * due to bucket else at end */
                else if( p.nMatch("ORTH" ) && p.nMatch("PARA" ) )
                {
                        /* option to turn ortho - para collisions with particles off */
                        if( p.nMatch(" ON " ) )
                        {
                                /* this is the default, leave collisions off */
                                diatom->lgH2_ortho_para_coll_on = true;
                        }
                        else
                        {
                                /* default (and only reason for this command) is to turn off 
                                 * ortho-para collisions */
                                diatom->lgH2_ortho_para_coll_on = false;
                        }
                }

                /* option to turn collisional effects off or on */
                else if( p.nMatch("GRAI" ) )
                {
                        /* option to turn collisions off */
                        if( p.nMatch(" ON" ) )
                        {
                                /* this is the default, leave collisions off */
                                diatom->lgH2_grain_deexcitation = true;
                        }
                        else
                        {
                                /* default (and only reason for this command) is to turn off collisions */
                                diatom->lgH2_grain_deexcitation = false;
                        }
                }
                else if( p.nMatch(" HE " ) )
                {
                        /* atom H2 He collisions ORNL (the default), Le BOURlot, and OFF
                         * which data set for He collisions,
                         * Teck Lee et al. ApJ to be submitted */
                        if( p.nMatch(" NEW" ) || p.nMatch("ORNL" ) )
                        {
                                /* use the new coefficients */
                                diatom->lgH2_He_ORNL = true;
                                diatom->coll_source[1].filename = "coll_rates_He_ORNL.dat";
                        }
                        else if( p.nMatch(" OLD" ) || p.nMatch("BOUR" ) )
                        {
                                /* use the coefficients from
                                 *>>refer       H2      collision       Le Bourlot, J., Pineau des Forets, 
                                 *>>refercon    G., & Flower, D.R. 1999, MNRAS, 305, 802*/
                                diatom->lgH2_He_ORNL = false;
                                diatom->coll_source[1].filename = "coll_rates_He_LeBourlot.dat";
                        }
                        else
                        {
                                fprintf( ioQQQ, 
                                        " I did not find a keyword on this ATOM H2 HE command - I know about the keys ORNL and Le BOURlot\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                }

                /*>>chng 08 feb 27, GS*/
                else if( p.nMatch("ORH2" ) )
                {
                        /* atom H2 H2ortho collisions ORNL (the default), Le BOURlot, and OFF
                         * which data set for H2 collisions,
                         * Teck Lee et al. ApJ to be submitted */
                        if( p.nMatch("ORNL" ) )
                        {
                                /* use the new coefficients */
                                diatom->lgH2_ORH2_ORNL = true;
                                diatom->coll_source[2].filename = "coll_rates_H2ortho_ORNL.dat";
                        }
                        else if( p.nMatch("BOUR" ) )
                        {
                                /* use the coefficients from
                                 *>>refer       H2      collision       Le Bourlot, J., Pineau des Forets, 
                                 *>>refercon    G., & Flower, D.R. 1999, MNRAS, 305, 802*/
                                diatom->lgH2_ORH2_ORNL = false;
                                diatom->coll_source[2].filename = "coll_rates_H2ortho_LeBourlot.dat";
                        }
                        else
                        {
                                fprintf( ioQQQ, 
                                        " I did not find a keyword on this ATOM H2 ohH2 command - I know about the keys ORNL and Le BOURlot\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                }

                else if( p.nMatch("PAH2" ) )
                {
                        /* atom H2 H2ortho collisions ORNL (the default), Le BOURlot, and OFF
                         * which data set for He collisions,
                         * Teck Lee et al. ApJ to be submitted */
                        if(  p.nMatch("ORNL" ) )
                        {
                                /* use the new coefficients */
                                diatom->lgH2_PAH2_ORNL = true;
                                diatom->coll_source[3].filename = "coll_rates_H2para_ORNL.dat";
                        }
                        else if( p.nMatch("BOUR" ) )
                        {
                                /* use the coefficients from
                                 *>>refer       H2      collision       Le Bourlot, J., Pineau des Forets, 
                                 *>>refercon    G., & Flower, D.R. 1999, MNRAS, 305, 802*/
                                diatom->lgH2_PAH2_ORNL = false;
                                diatom->coll_source[3].filename = "coll_rates_H2para_LeBourlot.dat";
                        }
                        else
                        {
                                fprintf( ioQQQ, 
                                        " I did not find a keyword on this ATOM H2 paH2 command - I know about the keys ORNL and Le BOURlot\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                }

                else
                {
                        /* option to turn all collisions off */
                        if( p.nMatch(" ON " ) )
                        {
                                /* this is the default, leave collisions on */
                                diatom->lgColl_deexec_Calc = true;
                        }
                        else
                        {
                                /* default (and only reason for this command) is to turn off collisions */
                                diatom->lgColl_deexec_Calc = false;
                        }
                }
        }

        /* set number of levels in matrix, but not trace matrix option */
        else if( p.nMatch("MATR" ) && !p.nMatch("TRAC" ) )
        {
                /* matrix option sets the number of levels that will
                 * be included in the matrix solution */
                long numLevels = (long)p.FFmtRead();
                if( p.nMatch(" ALL") )
                {
                        /* " all" means do all of X, but space has not yet been allocated,
                         * so we do know know how many levels are within X - set special
                         * flag that will be used then this is known */
                        numLevels = -1;
                }
                else if( p.lgEOL() && !(p.nMatch(" OFF") || p.nMatch("NONE") ) )
                {
                        /* this branch hit eol but OFF or NONE is not on line - this is a mistake */
                        fprintf( ioQQQ, 
                                " The total number of levels used in the matrix solver must be entered, or keywords ALL or NONE entered.\n Sorry.\n");
                        cdEXIT(EXIT_FAILURE);
                }

                diatom->set_numLevelsMatrix( numLevels );
                /* cannot check less than total number of levels within x since not yet set 
                 * We do not certify that matrix limits are greater than 1 -
                 * zero or <0 limits just turns if off, as did the off option */
        }
        else if( p.nMatch(" LTE" ) )
        {
                /* LTE option causes code to assume LTE for level populations  */
                diatom->lgLTE = true;
        }

        else if( p.nMatch("TRAC" ) )
        {
                /* these are used to set trace levels of output 
                diatom->n_trace_final = 1;
                diatom->n_trace_iterations = 2;
                diatom->n_trace_full = 3;
                diatom->n_trace_matrix = 4*/

                /* turns on trace printout - there are multiple levels */
                if( p.nMatch("FINA" ) )
                {
                        /* FINAL gives only final information when solver exits */
                        diatom->nTRACE = diatom->n_trace_final;
                }
                else if( p.nMatch("ITER" ) )
                {
                        /* follow iterations within each call */
                        diatom->nTRACE = diatom->n_trace_iterations;
                }
                else if( p.nMatch("FULL" ) )
                {
                        /* full details of solution - this is also the default*/
                        diatom->nTRACE = diatom->n_trace_full;
                }
                else if( p.nMatch("MATR" ) )
                {
                        /* print the matrices used for X */
                        diatom->nTRACE = diatom->n_trace_matrix;
                }
                else
                {
                        /* full details of solution is also the default*/
                        diatom->nTRACE = diatom->n_trace_full;
                }
        }
        else if( p.nMatch("NOIS" ) )
        {
                unsigned int iseed;
                /* check on effects of uncertainties in collision rates */
                diatom->lgH2_NOISE = true;
                diatom->lgH2_NOISECOSMIC = true;

                /* optional mean - default is 0 */
                diatom->xMeanNoise = p.FFmtRead();
                if( p.lgEOL() )
                        diatom->xMeanNoise = 0.;

                /* this is the standard deviation for the mole, with default */
                diatom->xSTDNoise = p.FFmtRead();
                if( p.lgEOL() )
                        diatom->xSTDNoise = 0.5;

                /* this may be a seed for the random number generator.  if no seed is
                 * set then use system time, and always get different sequence */
                iseed = (unsigned int)p.FFmtRead();
                /* returned 0 if eol hit */
                if( iseed > 0 )
                {
                        /* user set seed */
                        init_genrand( iseed );
                }
                else
                {
                        init_genrand( (unsigned)time( NULL ) );
                }
        }

        else if( p.nMatch("THER" ) )
        {
                /* change the treatment of the heating - cooling effects of H2,
                 * options are simple (use TH85 expressions) and full (use large molecule)*/
                if( p.nMatch("SIMP" ) )
                {
                        hmi.lgH2_Thermal_BigH2 = false;
                }
                else if( p.nMatch("FULL" ) )
                {
                        /* this is the default - use big atom */
                        hmi.lgH2_Thermal_BigH2 = true;
                }
        }

        else if( p.nMatch("CHEM" ) )
        {
                /* atom h2 chemistry simple command
                 * change the treatment of the chemistry - formation and destruction,
                 * options are simple (use TH85 expressions) and full (use large molecule)*/
                if( p.nMatch("SIMP" ) )
                {
                        hmi.lgH2_Chemistry_BigH2 = false;
                }
                else if( p.nMatch("FULL" ) )
                {
                        /* this is the default - use big atom */
                        hmi.lgH2_Chemistry_BigH2 = true;
                }
        }

        /* there is no final branch - if we do not find a keyword, simply
         * turn on the H2 molecule */
        return;
}