parse_drive.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 */
/*ParseDrive parse the drive command - drive calls to various subs */
/*DrvCaseBHS allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */
/*DrvHyas allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */
/*dgaunt drive gaunt factor routines by letting user query values */
#include "cddefines.h"
#include "trace.h"
#include "hydro_bauman.h"
#include "hydroeinsta.h"
#include "thirdparty.h"
#include "phycon.h"
#include "atmdat.h"  
#include "taulines.h"
#include "thermal.h"
#include "thirdparty.h"
#include "abund.h"
#include "rt.h"
#include "continuum.h"
#include "parser.h"
#include "physconst.h"
 
/*dgaunt drive gaunt factor routines by letting user query values */
STATIC void dgaunt(void);
 
/*DrvHyas allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */
STATIC void DrvHyas(void);
 
/* drive escape probability routines */
STATIC void DrvEscP( void );
 
/*DrvCaseBHS allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */
STATIC void DrvCaseBHS(void);
 
void ParseDrive(Parser &p )
{
        bool lgEOL;
        long int n, 
          i;
        double fac, 
          zed;
 
        DEBUG_ENTRY( "ParseDrive()" );
 
        /* NB evolve all following names to style DrvSomething */
 
        /* option to drive cloudy, which one? */
        if( p.nMatch("FFMT") )
        {
                /* free format parser */
                char chInput[INPUT_LINE_LENGTH];
                fprintf( ioQQQ, " FFmtRead ParseDrive entered.  Enter number.\n" );
                lgEOL = false;
                while( !lgEOL )
                {
                        if( read_whole_line( chInput , (int)sizeof(chInput) , ioStdin ) == NULL )
                        {
                                fprintf( ioQQQ, " ParseDrive.dat error getting magic number\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        i = 1;
                        fac = FFmtRead(chInput,&i,sizeof(chInput),&lgEOL);
                        if( lgEOL )
                        {
                                fprintf( ioQQQ, " FFmtRead hit the EOL with no value, return=%10.2e\n", 
                                  fac );
                                break;
                        }
                        else if( fac == 0. )
                        {
                                break;
                        }
                        else
                        {
                                fprintf( ioQQQ, " FFmtRead returned with value%11.4e\n", 
                                  fac );
                        }
                        fprintf( ioQQQ, " Enter 0 to stop, or another value.\n" );
                }
                fprintf( ioQQQ, " FFmtRead ParseDrive exits.\n" );
        }
 
        else if( p.nMatch("CASE") )
        {
                /* option to interpolate on Hummer and Storey case b hydrogenic emission routines */
                DrvCaseBHS( );
        }
 
        else if( p.nMatch("CDLI") )
        {
                /* drive cdLine to check that it finds all the right lines, routine is in lines.c */
                trace.lgDrv_cdLine = true;
        }
 
        else if( p.nMatch(" E1 ") )
        {
                // option to drive exponential integral routines
                // first, special case given in Abramovitc & Stegan
                double tau = 1.275;
                for( i=0; i<50; ++i )
                {
                        fprintf(ioQQQ,"tau\t%.3e\t exp-tau\t%.5e\t e1 tau\t%.5e  \t e2 "
                                "\t%.5e \te2n %.5e \t e3\t%.5e \t e4\t%.5e \n",
                                tau , sexp(tau) , ee1(tau) , e2(tau ), expn(2, tau) ,
                                expn(3 , tau ), expn(4 , tau ) );
                        tau = pow( 10. , ((double)i/4. - 9.) );
                }
                cdEXIT(EXIT_SUCCESS);
        }
 
        else if( p.nMatch("ESCA") )
        {
                /* option to drive escape probability routines */
                DrvEscP( );
        }
 
        else if( p.nMatch("HYAS") )
        {
                /* option to drive Jason's hydrogen transition probabilities */
                DrvHyas();
        }
 
        else if( p.nMatch("GAUN") )
        {
                /* drive gaunt factor routine */
                dgaunt();
        }
 
        else if( p.nMatch("POIN") )
        {
                /* later on, check cell pointers, centers, widths */
                fprintf( ioQQQ, " Define entire model first, later I will ask for energies.\n" );
                trace.lgPtrace = true;
        }
 
        else if( p.nMatch("PUMP") )
        {
                char chInput[INPUT_LINE_LENGTH];
                lgEOL = false;
                fprintf( ioQQQ, " Continuum pump ParseDrive entered - Enter log tau\n" );
                while( !lgEOL )
                {
                        if( read_whole_line( chInput , (int)sizeof(chInput) , ioStdin ) == NULL )
                        {
                                fprintf( ioQQQ, " ParseDrive.dat error getting magic number\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        /* get tau */
                        i = 1;
                        fac = FFmtRead(chInput,&i,sizeof(chInput),&lgEOL);
                        if( lgEOL )
                                break;
                        fac = pow(10.,fac);
                        fprintf( ioQQQ, " Tau =%11.4e\n", fac );
                        (*TauDummy).Emis().TauIn() = (realnum)fac;
                        // not sure what this is supposed to do, but TauDummy.Hi->nelem() and therefore the doppler width are ill-defined here!
                        // just put 1 for hydrogen
                        fac = DrvContPump(*TauDummy,1.f); 
                        fprintf( ioQQQ, " ContPump =%11.4e\n", fac );
                        fprintf( ioQQQ, " Enter null to stop, or another value.\n" );
                }
                fprintf( ioQQQ, " ContPump ParseDrive exits.\n" );
        }
 
        else if( p.nMatch("STAR") )
        {
                char chInput[INPUT_LINE_LENGTH];
                /* get starburst abundances */
                for( i=0; i < 40; i++ )
                {
                        zed = ((double)i+1.)/4. + 0.01;
                        sprintf( chInput, "starburst, zed=%10.4f", zed );
                        p.setline(chInput);
                        abund_starburst(p);
                        fprintf( ioQQQ, "%8.1e", zed );
                        for(n=0; n < LIMELM; n++)
                                fprintf( ioQQQ, "%8.1e", abund.solar[n] );
                        fprintf( ioQQQ, "\n" );
                }
        }
 
        else if( p.nMatch("VOIGT") )
        {
                /* create tab-delimited table giving Voigt function */
                FILE *ioVOIGT = fopen("voigt.dat" , "w");
                fprintf(ioVOIGT,"x \\ a");
                const realnum DampLogMin = -4., DampLogMax = 4.01;
                for( realnum damplog=DampLogMin; damplog<DampLogMax; ++damplog)
                        fprintf(ioVOIGT,"\tlog a=%.3e",pow(10.,damplog));
                fprintf(ioVOIGT , "\n");
 
                for( realnum x=-2.; x<5.;x+=0.05)
                {
                        realnum xlin = (realnum)pow(10.,x);
                        fprintf(ioVOIGT,"%.3e",xlin);
                        for( realnum damplog=DampLogMin; damplog<DampLogMax; ++damplog)
                        {
                                realnum xval[1];
                                xval[0] = xlin;
                                realnum damp = (realnum)pow(10. , damplog);
                                realnum yval[1];
                                VoigtH(damp,xval,yval,1);
                                fprintf(ioVOIGT , "\t%.3e",yval[0]);
                        }
                        fprintf(ioVOIGT , "\n");
                }
                fclose(ioVOIGT);
                cdEXIT(EXIT_SUCCESS);
        }
 
        else
        {
                fprintf( ioQQQ, 
                        " Unrecognized key; keys are CASE, CDLIne, E1 , ESCApe, FFMTread, GAUNt, "
                        "HYAS, POINt, PUMP, STAR, and VOIGt.  Sorry.\n" );
                cdEXIT(EXIT_FAILURE);
        }
        return;
}
 
/*DrvEscP user queries escape probability routines, which return values */
STATIC void DrvEscP( void )
{
        char chCard[INPUT_LINE_LENGTH];
        bool lgEOL;
        long i;
        double tau;
 
        DEBUG_ENTRY( "DrvEscP()" );
 
        /* this routine is enterd with the command escape probability, and
         * drives the escape probability routine to compare answers */
        fprintf( ioQQQ, " Enter the log of the one-sided optical depth; line with no number to stop.\n" );
 
        lgEOL = false;
        while( !lgEOL )
        {
                if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
                {
                        break;
                }
 
                i = 1;
                tau = FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
                if( lgEOL )
                {
                        break;
                }
 
                tau = pow(10.,tau);
                fprintf( ioQQQ, "tau was %e\n", tau );
                fprintf( ioQQQ, " ESCINC=%11.3e\n", esc_PRD_1side(tau,1e-4) );
                fprintf( ioQQQ, " ESCCOM=%11.3e\n", esc_CRDwing_1side(tau,1e-4 ) );
                fprintf( ioQQQ, " ESCA0K2=%11.3e\n", esca0k2(tau) );
 
        }
        return;
}
 
/*DrvCaseBHS allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */
STATIC void DrvCaseBHS(void)
{
        char chCard[INPUT_LINE_LENGTH];
        bool lgEOL,
                lgHit;
        long int i, 
          n1, 
          nelem ,
          n2;
        double Temperature,
                Density;
 
        DEBUG_ENTRY( "DrvCaseBHS()" );
 
        /* this routine is entered with the command DRIVE CASEB, and
         * interpolates on the Hummer & Storey case b data set */
 
        /* read in some external data files, but only if this is first call */
        fprintf(ioQQQ," I will get needed H data files. This will take a second.\n");
        atmdat_readin();
 
        {
                /* following should be set true to print input lines */
                /*@-redef@*/
                enum {DEBUG_LOC=false};
                /*@+redef@*/
                if( DEBUG_LOC )
                {
                        double xLyman , alpha;
                        long int ipHi;
                        nelem = 2;
                        Temperature = 2e4;
                        Density = 1e2;
                        for( ipHi=3; ipHi<25; ++ipHi )
                        {
                                double photons = (1./POW2(ipHi-1.)-1./POW2((double)ipHi) ) /(1.-1./ipHi/ipHi );
                                xLyman = atmdat_HS_caseB(1,ipHi, nelem,Temperature , Density , 'A' );
                                alpha = atmdat_HS_caseB(ipHi-1,ipHi, nelem,Temperature , Density , 'A' );
                                fprintf(ioQQQ,"%li\t%.3e\t%.3e\n", ipHi, xLyman/alpha*photons, photons );
                        }
                        cdEXIT(EXIT_SUCCESS);
                }
        }
 
        /* first get the charge, only H and He at present */
        lgHit = false;
        nelem = 0;
        while( !lgHit )
        {
                fprintf( ioQQQ, " Enter atomic number of species, either 1(H) or 2(He).\n" );
                if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
                {
                        fprintf( ioQQQ, " error getting species \n" );
                        cdEXIT(EXIT_FAILURE);
                }
 
                i = 1;
                nelem = (long int)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
                if( lgEOL || nelem< 1 || nelem > 2 )
                {
                        fprintf( ioQQQ, " must be either 1 or 2!\n" );
                }
                else
                {
                        lgHit = true;
                }
        }
 
        fprintf(ioQQQ," In the following temperatures <10 are log, >=10 linear.\n");
        fprintf(ioQQQ," The density is always a log.\n");
        fprintf(ioQQQ," The order of the quantum numbers do not matter.\n");
        fprintf(ioQQQ," The smallest must not be smaller than 2,\n");
        fprintf(ioQQQ," and the largest must not be larger than 25.\n");
        fprintf(ioQQQ," Units of emissivity are erg cm^3 s^-1\n\n");
        fprintf(ioQQQ," The limits of the HS tables are 2 <= n <= 25.\n");
 
        lgHit = true;
        /* this is always true */
        while( lgHit )
        {
                fprintf( ioQQQ, " Enter 4 numbers, temperature, density, 2 quantum numbers, null line stop.\n" );
                if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
                {
                        fprintf( ioQQQ, " Thanks for interpolating on the Hummer & Storey data set!\n" );
                        cdEXIT(EXIT_FAILURE);
                }
 
                i = 1;
                Temperature = FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
                if( lgEOL )
                {
                        fprintf( ioQQQ, " error getting temperature!\n" );
                        break;
                }
 
                /* log if less than 10 */
                if( Temperature < 10. )
                {
                        Temperature = pow(10., Temperature );
                }
                fprintf(ioQQQ," Temperature is %g\n", Temperature );
 
                Density = FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
                if( lgEOL )
                {
                        fprintf( ioQQQ, " error getting density!\n" );
                        break;
                }
                Density = pow(10., Density );
                fprintf(ioQQQ," Density is %g\n", Density );
 
                /* these quantum numbers can be in any order */
                n1 = (long)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
                if( lgEOL )
                {
                        fprintf( ioQQQ, " error getting quantum number!\n" );
                        break;
                }
 
                n2 = (long)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
                if( lgEOL )
                {
                        fprintf( ioQQQ, " error getting quantum number!\n" );
                        break;
                }
 
                if( MAX2( n1 , n2 ) > 25 )
                {
                        fprintf( ioQQQ," The limits of the HS tables are 2 <= n <= 25.  Sorry.\n");
                        break;
                }
 
                fprintf( ioQQQ, 
                        " 4pJ(%ld,%ld)/n_e n_p=%11.3e\n", 
                         n1, n2, 
                         atmdat_HS_caseB(n1,n2, nelem,Temperature , Density , 'B' ) );
 
                /* this is check that we were in bounds */
#if     0
                {
                        long j;
                        double tempTable[33] = {
                                        11870.,12490.,12820.,
                                        11060.,17740.,12560.,
                                        16390.,16700.,11360.,
                                        10240.,20740.,12030.,
                                        14450.,19510.,12550.,
                                        16470.,16560.,12220.,
                                        15820.,12960.,10190.,
                                        12960.,14060.,12560.,
                                        11030.,10770.,13360.,
                                        10780.,11410.,11690.,
                                        12500.,13190.,21120. };
                        double edenTable[33] = {
                                        10.,270.,80.,10.,70.,
                                        110.,200.,10.,40.,90.,
                                        340.,80.,60.,340.,30.,
                                        120.,10.,50.,450.,30.,
                                        180.,20.,170.,60.,20.,
                                        40.,30.,20.,100.,130.,
                                        10.,10.,110. };
 
 
                        for( j=0; j<33; j++ )
                        {
                                double halpha, hbeta, hgamma;
 
                                halpha = atmdat_HS_caseB(2,3, 1,tempTable[j] , edenTable[j] , 'B' );
                                hbeta =  atmdat_HS_caseB(2,4, 1,tempTable[j] , edenTable[j] , 'B' );
                                hgamma = atmdat_HS_caseB(2,5, 1,tempTable[j] , edenTable[j] , 'B' );
 
                                fprintf( ioQQQ, "%e\t%e\t%e\t%e\n",
                                        tempTable[j], 
                                        edenTable[j],
                                        halpha/hbeta,
                                        hgamma/hbeta );
                        }
                }
#endif
        }
 
        fprintf( ioQQQ, " Thanks for interpolating on the Hummer & Storey data set!\n" );
        cdEXIT(EXIT_FAILURE);
 
}
 
/*DrvHyas allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */
STATIC void DrvHyas(void)
{
        char chCard[INPUT_LINE_LENGTH];
        bool lgEOL;
        long int i, nHi, lHi, nLo, lLo;
 
        DEBUG_ENTRY( "DrvHyas()" );
 
        /* this routine is entered with the command DRIVE HYAS, and
         * drives Jason's hydrogen einstein A routines */
 
        nHi = 1;
        /* nHi never lt 1 */
        while( nHi != 0 )
        {
                fprintf( ioQQQ, " Enter four quantum numbers (n, l, n', l'), null line to stop.\n" );
                if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
                {
                        fprintf( ioQQQ, " error getting drvhyas \n" );
                        cdEXIT(EXIT_FAILURE);
                }
 
                i = 1;
                nHi = (long int)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
                if( lgEOL )
                        break;
 
                lHi = (long int)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
                if( lgEOL )
                {
                        fprintf( ioQQQ, " must be four numbers!\n" );
                        break;
                }
 
                if( lHi >= nHi )
                {
                        fprintf( ioQQQ, " l must be less than n!\n" );
                        break;
                }
 
                nLo = (long int)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
                if( lgEOL )
                {
                        fprintf( ioQQQ, " must be four numbers!\n" );
                        break;
                }
 
                lLo = (long int)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
                if( lgEOL )
                {
                        fprintf( ioQQQ, " must be four numbers!\n" );
                        break;
                }
 
                if( lLo >= nLo )
                {
                        fprintf( ioQQQ, " l must be less than n!\n" );
                        break;
                }
 
                if( nLo > nHi )
                {
                        long nTemp, lTemp;
 
                        /* swap hi and lo */
                        nTemp = nLo;
                        lTemp = lLo;
                        nLo = nHi;
                        lLo = lHi;
                        nHi = nTemp;
                        lHi = lTemp;
                }
 
                fprintf( ioQQQ, " A(%3ld,%3ld->%3ld,%3ld)=%11.3e\n", 
                        nHi, lHi, nLo, lLo,
                        H_Einstein_A( nHi, lHi, nLo, lLo, 1 ) );
 
        }
        fprintf( ioQQQ, " Driver exits, enter next line.\n" );
 
        return;
}
 
/*dgaunt drive gaunt factor routines by letting user query values */
STATIC void dgaunt(void)
{
        char chCard[INPUT_LINE_LENGTH];
        bool lgEOL;
        int inputflag;
        long int i, 
          ierror;
        realnum enerlin[1];
        double SaveTemp;
        double z,mygaunt=0.;
        double loggamma2, logu;
 
        DEBUG_ENTRY( "dgaunt()" );
 
        SaveTemp = phycon.te;
 
        /* this routine is entered with the command DRIVE GAUNT, and
         * drives the gaunt factor routine to check range
         * */
        fprintf( ioQQQ, " Enter 0 to input temp, energy, and net charge, or 1 for u and gamma**2.\n" );
        if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
        {
                fprintf( ioQQQ, " dgaunt error getting magic number\n" );
                cdEXIT(EXIT_FAILURE);
        }
        i = 1;
        inputflag = (int)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
 
        if( inputflag == 0 )
        {
                fprintf( ioQQQ, " Enter the temperature (log if <=10), energy (Ryd), and net charge. Null line to stop.\n" );
                /* >>chng 96 july 07, got rid of statement labels replacing with do while
                 * */
                ierror = 0;
                while( ierror == 0 )
                {
                        if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
                        {
                                fprintf( ioQQQ, " dgaunt error getting magic number\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        i = 1;
                        phycon.alogte = FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
                        /* the line may be trash but ierror will pick it up  */
                        if( lgEOL  )
                        {
                                fprintf( ioQQQ, " Gaunt driver exits, enter next line.\n" );
                                break;
                        }
                        /* numbers less than or equal to 10 are the log of the temperature */
                        double TeNew;
                        if( phycon.alogte > 10. )
                        {
                                TeNew = phycon.alogte;
                        }
                        else
                        {
                                TeNew = pow(10.,phycon.alogte);
                        }
                        TempChange(TeNew , false);
 
                        enerlin[0] = (realnum)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
                        if( lgEOL || enerlin[0] == 0. )
                        {
                                fprintf( ioQQQ, " Sorry, but there should be two more numbers, energy and charge.\n" );
                        }
 
                        z = (double)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
                        if( lgEOL || z == 0. )
                        {
                                fprintf( ioQQQ, " Sorry, but there should be a third number, charge.\n" );
                        }
 
                        /* This is non-thermally averaged gaunt factors.        */
                        mygaunt = cont_gaunt_calc( (double)phycon.te, z, enerlin[0] );
 
                        fprintf( ioQQQ, " Using my  routine, Gff= \t" );
                        fprintf( ioQQQ, "%11.3e\n", mygaunt );
 
                }
        }
        else
        {
                /* this routine is entered with the command DRIVE GAUNT, and
                 * drives the gaunt factor routine to check range
                 * */
                fprintf( ioQQQ, " Enter log u and log gamma2. Null line to stop.\n" );
                ierror = 0;
                while( ierror == 0 )
                {
                        if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
                        {
                                fprintf( ioQQQ, " dgaunt error getting magic number\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        i = 1;
                        logu = FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
                        /* the line may be trash but ierror will pick it up  */
                        if( lgEOL  )
                        {
                                fprintf( ioQQQ, " Gaunt driver exits, enter next line.\n" );
                                break;
                        }
 
                        loggamma2 = FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
                        if( lgEOL  )
                        {
                                fprintf( ioQQQ, " Sorry, but there should be another numbers, log gamma2.\n" );
                        }
 
                        /* This is my attempt to calculate non-thermally averaged gaunt factors.        */
                        mygaunt = cont_gaunt_calc( TE1RYD/pow(10.,loggamma2), 1., pow(10.,logu-loggamma2) );
 
                        TempChange(TE1RYD/pow(10.,loggamma2) , false);
 
                        fprintf( ioQQQ, " Using my  routine, Gaunt factor is" );
                        fprintf( ioQQQ, "%11.3e\n", mygaunt );
                }
        }
 
        TempChange(SaveTemp , false);
        return;
}