save_average.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 */
/*save_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 "save.h"
#include "parser.h"
 
/* return value is number of lines, -1 if file could not be opened */
void parse_save_average( 
        Parser &p,
        /* the file we will write to */
        long int ipPun, 
        char *chHeader)
{
        long int i;
        long int nLine;
 
        DEBUG_ENTRY( "parse_save_average()" );
 
        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 save structure.  These are output in a later call
         * to this routine with argument PUNCH  */
        
        /* number of lines we will save */
        save.nAverageList[ipPun] = 0;
        
        /* get the next line, and check for eof */
        p.getline();
        if( p.m_lgEOF )
        {
                fprintf( ioQQQ, 
                                        " Punch average hit EOF while reading list; use END to end list.\n" );
                cdEXIT(EXIT_FAILURE);
        }
                
        /* keep reading until we hit END */
        while( p.strcmp( "END" ) != 0 )
        {
                /* count number of species we will save */
                ++save.nAverageList[ipPun];
                
                /* get next line and check for eof */
                p.getline();
                if( p.m_lgEOF )
                {
                        fprintf( ioQQQ, " save averages hit EOF while reading species list; use END to end list.\n" );
                        cdEXIT(EXIT_FAILURE);
                }               
        }
/*#             define PADEBUG*/
#               ifdef PADEBUG
        fprintf(ioQQQ , "DEBUG save_average %li species read in.\n", 
                          save.nAverageList[ipPun] );
#               endif
        
        /* now create space that will be needed to hold these arrays */
        
        if( save.ipPnunit[ipPun] == NULL )
        {
                /* make sure the memory from a previous call is freed */
                save.SaveAverageFree(ipPun);
 
                /* nAverageIonList is set of ions for averages */
                save.nAverageIonList[ipPun].resize(save.nAverageList[ipPun]);
                
                /* nAverage2ndPar is set of second parameters for averages */
                save.nAverage2ndPar[ipPun].resize(save.nAverageList[ipPun]);
 
                /* chAverageType is label for type of average */
                save.chAverageType[ipPun].resize(save.nAverageList[ipPun]);
                
                /* chAverageSpeciesLabel is label for species */
                save.chAverageSpeciesLabel[ipPun].resize(save.nAverageList[ipPun]);
 
                for( i=0; i < save.nAverageList[ipPun]; ++i )
                {
                        /* create space for labels themselves */
                        save.chAverageType[ipPun][i] = new char[5];
                        /* chAverageSpeciesLabel is label for species */
                        save.chAverageSpeciesLabel[ipPun][i] = new char[5];
                }
        }
        
        /* reset array input read to first of the species lines */
        // Evil input rewind CodeSmell -- grrw
        input.nRead = nLine;
        
#               ifdef PADEBUG
        fprintf(ioQQQ , "DEBUG save_average %li species read in.\n", 
                          save.nAverageList[ipPun] );
#               endif
        
        /* reread the input lines and save the data */
        p.getline();
        if( p.m_lgEOF )
        {
                fprintf( ioQQQ, 
                                        " Punch average hit EOF while reading list; use END to end list.\n" );
                cdEXIT(EXIT_FAILURE);
        }
        
        /* use this to count number of species, and will assert equal to
         * total malloced above */
        nLine = 0;
        /* keep reading until we hit END */
        while( p.strcmp("END" ) != 0 )
        {
                /* count number of species we will save */
                ++nLine;
                if( p.nMatch("TEMP" ))
                {
                        /* temperature */
                                strcpy( save.chAverageType[ipPun][nLine-1] , "TEMP" );
                }
                else if( p.nMatch("COLU" ))
                {
                        /* column density */
                        strcpy( save.chAverageType[ipPun][nLine-1] , "COLU" );
                }
                else if( p.nMatch("IONI" ))
                {
                        /* ionization fraction */
                        strcpy( save.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 = p.GetElem( )) < 0 )
                {
                        /* no name found */
                        fprintf(ioQQQ, "save average did not an element on this line, sorry %s\n",
                                          chCap );
                        cdEXIT(EXIT_FAILURE);
                }
                strcpy( save.chAverageSpeciesLabel[ipPun][nLine-1] , elementnames.chElementNameShort[i]);
                
                /* now get ionization stage */
                save.nAverageIonList[ipPun][nLine-1] = (int) p.FFmtRead();
                if( p.lgEOL() )
                {
                        /* error - needed that ionization stage */
                        p.NoNumb("ionization stage" );
                }
                
                /* look for volume keyword, otherwise will be radius 
                 * only used for some options */
                if( p.nMatch( "VOLU" ) )
                {
                        /* volume */
                        save.nAverage2ndPar[ipPun][nLine-1] = 1;
                }
                else
                {
                        /* radius */
                        save.nAverage2ndPar[ipPun][nLine-1] = 0;
                }
                
                /* get next line and check for eof */
                p.getline();
                if( p.m_lgEOF )
                {
                        fprintf( ioQQQ, " save averages hit EOF while reading species list; use END to end list.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }
        
        /* these must equal or we have a major logical error */
        ASSERT( nLine == save.nAverageList[ipPun]);
        
#               ifdef PADEBUG
        for( i=0; i<nLine ; ++i )
        {
                fprintf(ioQQQ, "PDDEBUG %s %s %i %i\n",
                                  save.chAverageType[ipPun][i],
                                  save.chAverageSpeciesLabel[ipPun][i] ,
                                  save.nAverageIonList[ipPun][i] ,
                                  save.nAverage2ndPar[ipPun][i] );
        }
#               endif
        
        /* save headers */
        sprintf(chHeader, "#averages");
        for( i=0; i<nLine ; ++i )
        {
                sprintf(chTemp, "\t %s %s %i %i",
                                  save.chAverageType[ipPun][i],
                                  save.chAverageSpeciesLabel[ipPun][i] ,
                                save.nAverageIonList[ipPun][i] ,
                                  save.nAverage2ndPar[ipPun][i] );
                strcat( chHeader, chTemp );
        }
        strcat( chHeader, "\n");
}
 
void save_average( 
        /* the file we will write to */
        long int ipPun)
{
        long int i;
 
        DEBUG_ENTRY( "save_average()" );
 
                /* do the output */
                for( i=0; i<save.nAverageList[ipPun] ; ++i )
                {
                        double result;
                        char chWeight[7];
                        if( save.nAverage2ndPar[ipPun][i] == 0 )
                                strcpy( chWeight , "RADIUS");
                        else
                                strcpy( chWeight , "VOLUME");
 
                        if( strncmp( save.chAverageType[ipPun][i] , "TEMP" , 4 ) == 0)
                        {
                                /* temperature */
                                if( cdTemp( 
                                        save.chAverageSpeciesLabel[ipPun][i] ,
                                        save.nAverageIonList[ipPun][i] ,
                                        &result , 
                                        chWeight ) )
                                {
                                        fprintf( ioQQQ, " save average temperature could not identify the species.\n" );
                                        cdEXIT(EXIT_FAILURE);
                                }
                                /* will report log of temperature */
                                result = log10( result );
                        }
                        else if( strncmp( save.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" , 
                                        save.chAverageSpeciesLabel[ipPun][i] , 4)==0 &&
                                        save.nAverageIonList[ipPun][i]== 0 )
                                        strncpy( save.chAverageSpeciesLabel[ipPun][i],
                                        "H2  " , 4 );
                                if( cdIonFrac( 
                                        save.chAverageSpeciesLabel[ipPun][i] ,
                                        save.nAverageIonList[ipPun][i] ,
                                        &result , 
                                        chWeight ,
                                        false 
                                        ) )
                                {
                                        fprintf( ioQQQ, " save average ionization fraction could not identify the species.\n" );
                                        cdEXIT(EXIT_FAILURE);
                                }
                                /* will report log of ionization fraction */
                                result = log10( result );
                        }
                        else if( strncmp( save.chAverageType[ipPun][i] , "COLU" , 4 ) == 0)
                        {
                                /* column density */
                                if( cdColm( 
                                        save.chAverageSpeciesLabel[ipPun][i] ,
                                        save.nAverageIonList[ipPun][i] ,
                                        &result ) )
                                {
                                        fprintf( ioQQQ, " save 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(save.ipPnunit[ipPun], "\t %.3f", result );
                }
                fprintf(save.ipPnunit[ipPun], "\n");
}