/home66/gary/public_html/cloudy/c08_branch/source/punch_average.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*punch_average routine to read in list of species to output as averages */
#include "cddefines.h"
#include "cddrive.h"
#include "input.h"
#include "elementnames.h"
#include "punch.h"

/* return value is number of lines, -1 if file could not be opened */
void punch_average( 
        /* the file we will write to */
        long int ipPun, 
        /* this is the job we shall do - READ and PUNCH */
        const char *chJob,
        char *chHeader)
{
        long int i;
        bool lgEOF , lgEOL;
        long int nLine;

        char chLine[FILENAME_PATH_LENGTH_2] ;

        DEBUG_ENTRY( "punch_average()" );

        if( strcmp(chJob,"READ") == 0 )
        {
                char chCap[INPUT_LINE_LENGTH];
                char chTemp[INPUT_LINE_LENGTH];

                /* this is index for first line we will read.  use this to count 
                 * total number of species */
                nLine = input.nRead;
                /* very first time this routine is called, chJob is "READ" and we read
                * in lines from the input stream.  The species labels and other information
                * are store in the punch structure.  These are output in a later call
                * to this routine with argument PUNCH  */

                /* number of lines we will save */
                punch.nAverageList[ipPun] = 0;

                /* get the next line, and check for eof */
                input_readarray(chLine,&lgEOF);
                if( lgEOF )
                {
                        fprintf( ioQQQ, 
                                " Punch average hit EOF while reading list; use END to end list.\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                /* convert line to caps */
                strcpy( chCap, chLine );
                caps(chCap);

                /* keep reading until we hit END */
                while( strncmp(chCap, "END" ,3 ) != 0 )
                {
                        /* count number of species we will save */
                        ++punch.nAverageList[ipPun];

                        /* get next line and check for eof */
                        input_readarray(chLine,&lgEOF);
                        if( lgEOF )
                        {
                                fprintf( ioQQQ, " punch averages hit EOF while reading species list; use END to end list.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }

                        /* convert it to caps */
                        strcpy( chCap, chLine );
                        caps(chCap);
                }
/*#             define PADEBUG*/
#               ifdef PADEBUG
                fprintf(ioQQQ , "DEBUG punch_average %li species read in.\n", 
                punch.nAverageList[ipPun] );
#               endif

                /* now create space that will be needed to hold these arrays */

                if( punch.ipPnunit[ipPun] == NULL )
                {
                        /* nAverageIonList is set of ions for averages */
                        punch.nAverageIonList[ipPun] = (int *)MALLOC(
                                (size_t)(punch.nAverageList[ipPun]*sizeof(int)) );

                        /* nAverage2ndPar is set of second parameters for averages */
                        punch.nAverage2ndPar[ipPun] = (int *)MALLOC((size_t)
                                (punch.nAverageList[ipPun]*sizeof(int)) );

                        /* chAverageType is label for type of average */
                        punch.chAverageType[ipPun] = (char **)MALLOC((size_t)
                                (punch.nAverageList[ipPun]*sizeof(char *)) );

                        /* chAverageSpeciesLabel is label for species */
                        punch.chAverageSpeciesLabel[ipPun] = (char **)MALLOC((size_t)
                                (punch.nAverageList[ipPun]*sizeof(char *)) );
                        for( i=0; i<punch.nAverageList[ipPun]; ++i )
                        {
                                /* create space for labels themselves */
                                punch.chAverageType[ipPun][i] = (char *)MALLOC((size_t)
                                        (5*sizeof(char)) );

                                /* chAverageSpeciesLabel is label for species */
                                punch.chAverageSpeciesLabel[ipPun][i] = (char *)MALLOC((size_t)
                                        (5*sizeof(char)) );
                        }
                }

                /* reset array input read to first of the species lines */
                input.nRead = nLine;

#               ifdef PADEBUG
                fprintf(ioQQQ , "DEBUG punch_average %li species read in.\n", 
                        punch.nAverageList[ipPun] );
#               endif

                /* reread the input lines and save the data */
                input_readarray(chLine,&lgEOF);
                if( lgEOF )
                {
                        fprintf( ioQQQ, 
                                " Punch average hit EOF while reading list; use END to end list.\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                /* convert line to caps */
                strcpy( chCap, chLine );
                caps(chCap);

                /* use this to count number of species, and will assert equal to
                 * total malloced above */
                nLine = 0;
                /* keep reading until we hit END */
                while( strncmp(chCap, "END" ,3 ) != 0 )
                {
                        /* count number of species we will save */
                        ++nLine;
                        if( nMatch("TEMP" , chCap ))
                        {
                                /* temperature */
                                strcpy( punch.chAverageType[ipPun][nLine-1] , "TEMP" );
                        }
                        else if( nMatch("COLU" , chCap ))
                        {
                                /* column density */
                                strcpy( punch.chAverageType[ipPun][nLine-1] , "COLU" );
                        }
                        else if( nMatch("IONI" , chCap ))
                        {
                                /* ionization fraction */
                                strcpy( punch.chAverageType[ipPun][nLine-1] , "IONI" );
                        }
                        else
                        {
                                fprintf(ioQQQ,"PROBLEM one of the jobs TEMPerature, COLUmn density, or IONIzation, must appear.\n");
                                cdEXIT(EXIT_FAILURE);
                        }

                        /* get element name, a string we shall pass on to the routine
                         * that computes final quantities */
                        if( (i = GetElem( chCap )) < 0 )
                        {
                                /* no name found */
                                fprintf(ioQQQ, "punch average did not an element on this line, sorry %s\n",
                                        chCap );
                                cdEXIT(EXIT_FAILURE);
                        }
                        strcpy( punch.chAverageSpeciesLabel[ipPun][nLine-1] , elementnames.chElementNameShort[i]);

                        /* now get ionization stage */
                        i = 5;
                        punch.nAverageIonList[ipPun][nLine-1] = (int)
                                FFmtRead( chCap , &i,INPUT_LINE_LENGTH,&lgEOL );
                        if( lgEOF )
                        {
                                /* error - needed that ionization stage */
                                NoNumb( chCap );
                        }

                        /* look for volume keyword, otherwise will be radius 
                         * only used for some options */
                        if( nMatch( "VOLU" , chCap ) )
                        {
                                /* volume */
                                punch.nAverage2ndPar[ipPun][nLine-1] = 1;
                        }
                        else
                        {
                                /* radius */
                                punch.nAverage2ndPar[ipPun][nLine-1] = 0;
                        }

                        /* get next line and check for eof */
                        input_readarray(chLine,&lgEOF);
                        if( lgEOF )
                        {
                                fprintf( ioQQQ, " punch averages hit EOF while reading species list; use END to end list.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }

                        /* convert it to caps */
                        strcpy( chCap, chLine );
                        caps(chCap);
                }

                /* these must equal or we have a major logical error */
                ASSERT( nLine == punch.nAverageList[ipPun]);

#               ifdef PADEBUG
                for( i=0; i<nLine ; ++i )
                {
                        fprintf(ioQQQ, "PDDEBUG %s %s %i %i\n",
                                punch.chAverageType[ipPun][i],
                                punch.chAverageSpeciesLabel[ipPun][i] ,
                                punch.nAverageIonList[ipPun][i] ,
                                punch.nAverage2ndPar[ipPun][i] );
                }
#               endif

                /* punch headers */
                sprintf(chHeader, "#averages");
                for( i=0; i<nLine ; ++i )
                {
                        sprintf(chTemp, "\t %s %s %i %i",
                                punch.chAverageType[ipPun][i],
                                punch.chAverageSpeciesLabel[ipPun][i] ,
                                punch.nAverageIonList[ipPun][i] ,
                                punch.nAverage2ndPar[ipPun][i] );
                        strcat( chHeader, chTemp );
                }
                strcat( chHeader, "\n");
        }
        else if( strcmp(chJob,"PUNCH") == 0 )
        {
                /* do the output */
                for( i=0; i<punch.nAverageList[ipPun] ; ++i )
                {
                        double result;
                        char chWeight[7];
                        if( punch.nAverage2ndPar[ipPun][i] == 0 )
                                strcpy( chWeight , "RADIUS");
                        else
                                strcpy( chWeight , "VOLUME");

                        if( strncmp( punch.chAverageType[ipPun][i] , "TEMP" , 4 ) == 0)
                        {
                                /* temperature */
                                if( cdTemp( 
                                        punch.chAverageSpeciesLabel[ipPun][i] ,
                                        punch.nAverageIonList[ipPun][i] ,
                                        &result , 
                                        chWeight ) )
                                {
                                        fprintf( ioQQQ, " punch average temperature could not identify the species.\n" );
                                        cdEXIT(EXIT_FAILURE);
                                }
                                /* will report log of temperature */
                                result = log10( result );
                        }
                        else if( strncmp( punch.chAverageType[ipPun][i] , "IONI" , 4 ) == 0)
                        {
                                /* ionization fraction 
                                 * H2 is a special case, HYDRO 0 requests
                                 * the H2 fraction, n(H2)/n(H) */
                                if( strncmp( "HYDR" , 
                                        punch.chAverageSpeciesLabel[ipPun][i] , 4)==0 &&
                                        punch.nAverageIonList[ipPun][i]== 0 )
                                        strncpy( punch.chAverageSpeciesLabel[ipPun][i],
                                        "H2  " , 4 );
                                if( cdIonFrac( 
                                        punch.chAverageSpeciesLabel[ipPun][i] ,
                                        punch.nAverageIonList[ipPun][i] ,
                                        &result , 
                                        chWeight ,
                                        false 
                                        ) )
                                {
                                        fprintf( ioQQQ, " punch average ionization fraction could not identify the species.\n" );
                                        cdEXIT(EXIT_FAILURE);
                                }
                                /* will report log of ionization fraction */
                                result = log10( result );
                        }
                        else if( strncmp( punch.chAverageType[ipPun][i] , "COLU" , 4 ) == 0)
                        {
                                /* column density */
                                if( cdColm( 
                                        punch.chAverageSpeciesLabel[ipPun][i] ,
                                        punch.nAverageIonList[ipPun][i] ,
                                        &result ) )
                                {
                                        fprintf( ioQQQ, " punch average column density fraction could not identify the species.\n" );
                                        cdEXIT(EXIT_FAILURE);
                                }
                                /* will report log of column density */
                                result = log10( result );
                        }
                        else
                                TotalInsanity();

                        fprintf(punch.ipPnunit[ipPun], "\t %.3f", result );
                }
                fprintf(punch.ipPnunit[ipPun], "\n");
        }
        else
        {
                /* total insanity */
                TotalInsanity();
        }
        /* return */
        return;
}

---