save_species.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 */
/*SaveSpecies generate output for the save species command */
#include "cddefines.h"
#include "opacity.h"
#include "taulines.h"
#include "radius.h"
#include "phycon.h"
#include "save.h"
#include "mole.h"
 
/* save results for one particular species */
STATIC void SaveSpeciesOne( const size_t species_index, const char chKey[],
                FILE *ioPUN, long int ipPun, size_t maxLevels );
 
/*SaveSpecies generate output for the save species command */
void SaveSpecies(
                FILE* ioPUN,
                long int ipPun)
{
        DEBUG_ENTRY( "SaveSpecies()" );
 
        if( strcmp( save.chSaveArgs[ipPun], "LABE" )==0 )
        {
                if( save.lgPunHeader[ipPun] )
                {
                        /* save list of species labels */
                        fprintf( ioPUN, "#Species labels\n" );
                        save.lgPunHeader[ipPun] = false;
                        for( size_t i=0; i<mole_global.list.size(); ++i )
                        {
                                molecule *spg = &(*mole_global.list[i]);
                                fprintf( ioPUN, "%s\n", spg->label.c_str() );
                        }
                }
                return;
        }
 
        else if( strcmp(save.chSaveArgs[ipPun],"LEVL") == 0 )
        {
                /* number of levels active in this zone */
                if( save.lgPunHeader[ipPun] )
                {
                        /* save list of species labels */
                        fprintf( ioPUN, "#Species\tnumber of levels\n" );
                        save.lgPunHeader[ipPun] = false;
                }
                for( size_t i=0; i<mole_global.list.size(); ++i )
                {
                        molecule *spg = &(*mole_global.list[i]);
                        molezone *sp = &mole.species[i];
                        fprintf( ioPUN, "%s", spg->label.c_str() );
                        if( sp->levels == NULL )
                                fprintf( ioPUN, "\t%4lu\n", 0UL );
                        else
                                fprintf( ioPUN, "\t%4lu\n", (unsigned long)sp->levels->size() );
                }
                return;
        }
 
        /* remaining options are column densities, populations, and energies
         * first branch; save results for all species if "" */
        if( strcmp( save.chSaveSpecies[ipPun], "" ) == 0 )
        {
                // max number of levels, for header print
                size_t mostLevels = 0;
                for( size_t i=0; i<mole_global.list.size(); ++i )
                {
                        molezone *sp = &mole.species[i];
                        if( sp->levels != NULL )
                                mostLevels = MAX2(mostLevels, sp->levels->size() );
                }
                ASSERT( mostLevels > 1 );
                ASSERT( mostLevels < 10000 );
 
                // loop over species
                for( size_t i=0; i<mole_global.list.size(); ++i )
                        SaveSpeciesOne( i, save.chSaveArgs[ipPun], ioPUN, ipPun, mostLevels );
        }
        else
        {
                const molecule *saveSpeciesGlobal = findspecies(save.chSaveSpecies[ipPun]);
                const molezone *saveSpecies = findspecieslocal(save.chSaveSpecies[ipPun]);
 
                if( saveSpecies == null_molezone )
                {
                        fprintf( ioQQQ,"Could not find species %s, so SAVE SPECIES LABELS to get a list of all species."
                                        "\nSorry.\n", save.chSaveSpecies[ipPun] );
                        cdEXIT(EXIT_FAILURE);
                }
        
                size_t numLevels = 0;
                if( saveSpecies->levels != NULL )
                        numLevels = saveSpecies->levels->size();        
                SaveSpeciesOne( saveSpeciesGlobal->index, save.chSaveArgs[ipPun], ioPUN, ipPun, numLevels );
        }
 
        return;
}
 
/* print 0.000e+00 as simply 0 */
STATIC void PrintShortZero( FILE *ioPUN , double arg )
{
        DEBUG_ENTRY( "PrintShortZero()" );
        if( arg==0. )
                fprintf(ioPUN,"\t0");
        else
                fprintf(ioPUN,"\t%.3e", arg);
 
}
 
/* save results for one particular species */
STATIC void SaveSpeciesOne( const size_t species_index, const char chKey[],
                FILE *ioPUN, long int ipPun, size_t maxLevels )
{
        DEBUG_ENTRY( "SaveSpeciesOne()" );
 
        molecule *spg = &(*mole_global.list[species_index]);
        molezone *sp = &mole.species[species_index];
 
        if( spg == null_mole || sp == null_molezone )
                return;
        
        // one time print of energy levels
        if( strcmp( chKey, "ENER" )==0 )
        {
                if( save.lgPunHeader[ipPun] )
                {
                        save.lgPunHeader[ipPun] = false;
 
                        fprintf( ioPUN, "#species energies");
                        for( size_t i = 0; i < maxLevels; ++i )
                        {
                                fprintf( ioPUN, "\t%lu", (unsigned long)i );
                        }
                        fprintf( ioPUN, "\n");
                }
 
                fprintf( ioPUN, "%s", spg->label.c_str() );
                if( sp->levels == NULL || sp->levels->size() == 0 )
                {
                        fprintf( ioPUN, "\t%.6e", 0. );
                }
                else
                {
                        for( qList::const_iterator st = sp->levels->begin(); st != sp->levels->end(); ++st )
                        {
                                ASSERT( (*st).g() > 0.f );
                                fprintf( ioPUN, "\t%.6e", AnuUnit( (*st).energy().Ryd() ) );
                        }
                }
                fprintf( ioPUN, "\n");
                return;
        }
 
        if( strcmp( chKey, "POPU" )==0 )
        {
                if( save.lgPunHeader[ipPun] )
                {
                        fprintf( ioPUN, "#depth [cm] species populations [cm-3]");
 
                        for( size_t i = 0; i < maxLevels; ++i )
                        {
                                fprintf( ioPUN, "\t%lu", (unsigned long)i );
                        }
                        fprintf( ioPUN, "\n");
                        save.lgPunHeader[ipPun] = false;
                }
 
                fprintf( ioPUN, "%.5e", radius.depth_mid_zone );
                fprintf( ioPUN, "\t%s", spg->label.c_str() );
 
                if( sp->levels == NULL || sp->levels->size() == 0 )
                {
                        PrintShortZero( ioPUN, sp->den );
                }
                else
                {
                        bool lgZeroHit = false;
                        // loop over levels for this species
                        for( qList::const_iterator st = sp->levels->begin(); st != sp->levels->end(); ++st )
                        {
                                // don't print high levels which can have 0 abundance at low temperature
                                if( !lgZeroHit )
                                        PrintShortZero( ioPUN, (*st).Pop() );
                                if( (*st).Pop() == 0.)
                                        lgZeroHit = true;
                        }
                }
                fprintf( ioPUN, "\n");
        }
        else if( strcmp( chKey, "COLU" )==0 )
        {
 
                if( save.lgPunHeader[ipPun] )
                {
                        fprintf( ioPUN, "#species column density [cm-2]");
 
                        for( size_t i = 0; i < maxLevels; ++i )
                        {
                                fprintf( ioPUN, "\t%lu", (unsigned long)i );
                        }
                        fprintf( ioPUN, "\n");
                        save.lgPunHeader[ipPun] = false;
                }
 
                fprintf( ioPUN, "%s", spg->label.c_str() );
 
                if( sp->levels == NULL || sp->levels->size() == 0 )
                {
                        PrintShortZero( ioPUN, sp->column );
                }
                else
                {
                        // loop over levels
                        bool lgZeroHit = false;
                        for( qList::const_iterator st = sp->levels->begin(); st != sp->levels->end(); ++st )
                        {
                                // don't print high levels which can have 0 abundance at low temperature
                                if( !lgZeroHit )
                                        PrintShortZero( ioPUN, (*st).ColDen() );
                                if( (*st).ColDen() == 0.)
                                        lgZeroHit = true;
                        }
                }
                fprintf( ioPUN, "\n");
        }
        
        return;
}