/home66/gary/public_html/cloudy/c08_branch/source/iso_error.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 */
/*HeLikeError fills uncertainty arrays */
#include "cddefines.h" 
#include "iso.h"

/* This routine handles errors when that option is turned on (via the command
 * "atom he-like error generation" */
void iso_put_error(long int ipISO,
                          long int nelem,
                          long int ipHi,
                          long int ipLo,
                          long int whichData,
                          realnum errorToPut)
{

        DEBUG_ENTRY( "iso_put_error()" );

        if( iso.lgRandErrGen[ipISO] )
        {
                /* whichData is either IPRAD, IPCOLLIS, or IPENERGY */
                ASSERT( whichData <= 2 );
                ASSERT( ipISO <= NISO );
                ASSERT( nelem < LIMELM );
                ASSERT( ipHi <= iso.numLevels_max[ipISO][nelem] );
                ASSERT( ipLo <= iso.numLevels_max[ipISO][nelem] );
                ASSERT( errorToPut >= 0. );

                iso.Error[ipISO][nelem][ipHi][ipLo][whichData] = errorToPut;
        }
        return;
}

void iso_error_generation( long ipISO, long nelem )
{

        long ipHi, ipLo, typeOfRate;
        realnum ErrorToPut;

        DEBUG_ENTRY( "iso_error_generation()" );

        /* put recombination errors into iso.Error[ipISO] array */
        for( ipHi=0; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
        {
                if( N_(ipHi)<=5 )
                {
                        ErrorToPut = 0.02f;
                }
                else if( N_(ipHi)>iso.n_HighestResolved_max[ipISO][nelem] )
                {
                        ErrorToPut = 0.02f;
                }
                else if( L_(ipHi)>=3 )
                {
                        ErrorToPut = 0.005f;
                }
                else if( L_(ipHi)==2 )
                {
                        ErrorToPut = 0.01f;
                }
                else if( L_(ipHi)==1 )
                {
                        ErrorToPut = 0.025f;
                }
                else
                {
                        ErrorToPut = 0.05f;
                }
                iso_put_error(ipISO,nelem,iso.numLevels_max[ipISO][nelem],ipHi,IPRAD,ErrorToPut);
        }
        /* Error for total recombination */
        iso_put_error(ipISO,nelem,iso.numLevels_max[ipISO][nelem],iso.numLevels_max[ipISO][nelem],IPRAD,0.02f);

        for( ipHi=1; ipHi<= iso.numLevels_max[ipISO][nelem]; ipHi++ )
        {
                /* >>chng 06 mar 15, the upper limit incorrectly went to numLevels_max */
                for( ipLo=0; ipLo<= ipHi; ipLo++ )
                {
                        for( typeOfRate=0; typeOfRate<=1; typeOfRate++ )
                        {
                                if( iso.Error[ipISO][nelem][ipHi][ipLo][typeOfRate] >= 0. )
                                {
                                        iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][typeOfRate] =  
                                                (realnum)MyGaussRand( iso.Error[ipISO][nelem][ipHi][ipLo][typeOfRate] );
                                        ASSERT( iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][typeOfRate] > 0. );
                                }
                                else
                                {
                                        iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][typeOfRate] = 1.0f;
                                }
                        }
                }
        }

        /* set flag saying that error generation has been done.  */
        iso.lgErrGenDone[ipISO][nelem] = true;
        return;
}

---