/home66/gary/public_html/cloudy/c08_branch/source/grid_xspec.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 */
/*gridXspec handles all grid calculations, called by griddo */
/*gridFunc */
/*GridGatherInCloudy - gathers spectra for each grid calculation to save for end */
#include "cddefines.h"
#include "punch.h"
#include "warnings.h"
#include "optimize.h"
#include "cddrive.h"
#include "rfield.h"
#include "grid.h"
#include "called.h"
#include "prt.h"

/*gridXspec handles all grid calculations, called by grid_do */
void gridXspec(realnum xc[], long int nInterpVars)
{
        long int i, j, k;
        double averageChi2;

        DEBUG_ENTRY( "gridXspec()" );

        if( nInterpVars > LIMPAR )
        {
                fprintf( ioQQQ, "grid_do: too many parameters are varied, increase LIMPAR\n" );
                cdEXIT(EXIT_FAILURE);
        }

        optimize.nOptimiz = 0;
        grid.nintparm = nInterpVars;

        /* if this is changed there must be some change made to actually
         * stuff the additional parameter information. */
        grid.naddparm = 0;

        ASSERT( grid.nintparm + grid.naddparm <= LIMPAR );

        grid.totNumModels = 1;
        /* >>chng 06 aug 21, allow the number of parameter values to be different for different parameters. */
        for( i=0; i<nInterpVars; i++ )
        {
                /* >>chng 06 sep 4, use grid variable instead of passing to routine. */
                grid.totNumModels *= grid.numParamValues[i];
        }
        /* grid.totNumModels = (long)pow((double)nParVals, (double)nInterpVars); */
        ASSERT( grid.totNumModels > 1 );

        grid.paramNames = (char**)MALLOC(sizeof(char*)*(unsigned)(nInterpVars+grid.naddparm) );
        grid.paramMethods = (long*)MALLOC(sizeof(long)*(unsigned)(nInterpVars+grid.naddparm) );
        grid.paramRange = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nInterpVars+grid.naddparm) );
        grid.paramData = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nInterpVars+grid.naddparm) );
        grid.interpParameters = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(grid.totNumModels) );
        /* save abort flag */
        grid.lgAbort = (bool*)MALLOC(sizeof(bool)*(unsigned)(grid.totNumModels) );
        /* number of warnings */
        grid.lgWarn = (bool*)MALLOC(sizeof(bool)*(unsigned)(grid.totNumModels) );

        for( i=0; i<nInterpVars+grid.naddparm; i++ )
        {
                grid.paramNames[i] = (char*)MALLOC(sizeof(char)*(unsigned)(12) );
                grid.paramRange[i] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(6) );
                grid.paramData[i] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(grid.numParamValues[i]) );

                sprintf( grid.paramNames[i],    "%s%ld", "PARAM", i+1 );
                /* Method is 0 for linear, 1 for logarithmic, for now all linear. */
                grid.paramMethods[i] = 0;
                /* Initial */
                grid.paramRange[i][0] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f)/2.f;
                /* Delta */
                grid.paramRange[i][1] = grid.paramIncrements[i]/10.f;
                /* Minimum */
                grid.paramRange[i][2] = xc[i]-grid.paramIncrements[i]/10.f;
                /* Bottom */
                grid.paramRange[i][3] = xc[i]-grid.paramIncrements[i]/10.f;
                /* Top */
                grid.paramRange[i][4] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f)+grid.paramIncrements[i]/10.f;
                /* Maximum */
                grid.paramRange[i][5] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f)+grid.paramIncrements[i]/10.f;

                for( j=0; j<grid.numParamValues[i]; j++ )
                {
                        grid.paramData[i][j] = xc[i]+grid.paramIncrements[i]*j;
                }
        }

        for( i=0; i<grid.totNumModels; i++ )
        {
                grid.interpParameters[i] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nInterpVars) );
        }

        /* >>chng 06 aug 23, this logic now allows non-square parameter spaces.  */
        for( i=0; i< grid.totNumModels; i++ )
        {
                realnum variableVector[LIMPAR];

                for( j=0; j<nInterpVars; j++ )
                {
                        int index;
                        long volumeOtherDimensions = 1;

                        /* by "volume", we simply mean the product of the parameter dimensions 
                         * AFTER the present index, i.e.:
                         * first "volume" is product of grid.numParamValues[1]*grid.numParamValues[2]*....grid.numParamValues[n]
                         * second "volume" is product of grid.numParamValues[2]*grid.numParamValues[3]*....grid.numParamValues[n]
                         * last "volume" is unity.  */
                        for( k=j+1; k<nInterpVars; k++ )
                        {
                                volumeOtherDimensions *= grid.numParamValues[k];
                        }

                        /* For each successive parameter, the "volume" is less than the previous one.
                         * So the left-hand side of this modulus operation increases for each parameter,
                         * which causes the index of each parameter to be incremented more often than the
                         * index of the previous parameter.  Thus, the last dimension is incremented
                         * every time, then the second to last dimension is incremented with each repeat
                         * of the last dimension.  This repeats until, finally, the first dimension is
                         * incremented.  For example, the indices of the parameter vectors for a 2x2x3
                         * box would be ordered as such:
                         * [0][0][0]
                         * [0][0][1]
                         * [0][0][2]
                         * [0][1][0]
                         * [0][1][1]
                         * [0][1][2]
                         * [1][0][0]
                         * [1][0][1]
                         * [1][0][2]
                         * [1][1][0]
                         * [1][1][1]
                         * [1][1][2]
                         */
                        index = (int)( (i/volumeOtherDimensions)%(grid.numParamValues[j]) );

                        /* this prevents parameter incrementation for debugging purposes.  */
                        if( grid.lgStrictRepeat )
                                variableVector[j] = xc[j];
                        else
                                variableVector[j] = xc[j] + grid.paramIncrements[j]*index;

                        grid.interpParameters[i][j] = variableVector[j];
                }

                for( j=nInterpVars; j<LIMPAR; j++ )
                {
                        variableVector[j] = xc[j];
                }

                if( i == grid.totNumModels - 1 )
                {
                        called.lgTalk = true;
                        called.lgTalkIsOK = true;
                        prt.lgFaintOn = true;
                        grid.lgGridDone = true;
                }

                averageChi2 = optimize_func(variableVector);
                /* silly test to use the var averageChi2 for something - keeps lint and some
                 * compilers happy */
                if( averageChi2 < 0. && averageChi2 == 0 )
                        fprintf( ioQQQ , "DEBUG interesting impossible print has occurred.\n");

        }
        return;
}

/* GridGatherAfterCloudy - gather information just before
 * punch output generated */
void GridGatherAfterCloudy(
   /* chTime is null terminated 4 char string, either "MIDL" or "LAST" */
   const char *chTime)
{
        DEBUG_ENTRY( "GridGatherAfterCloudy()" );

        if( !grid.lgGrid )
        {
                /* not grid - return */
                return;
        }

        /* this has to be either MIDL or LAST */
        if( chTime[0] == 'L' )
        {
                ASSERT( optimize.nOptimiz>= 0 && optimize.nOptimiz<grid.totNumModels );
                grid.lgAbort[optimize.nOptimiz] = lgAbort;
                grid.lgWarn[optimize.nOptimiz] = warnings.lgWarngs;
        }
        else if( chTime[0] != 'M' )
                TotalInsanity();
        return;
}

/*GridGatherInCloudy - gathers spectra for each grid calculation to save for end */
void GridGatherInCloudy( void )
{
        long i;
        static bool lgFirstRun = true;

        DEBUG_ENTRY( "GridGatherInCloudy()" );

        if( !grid.lgGrid )
        {
                TotalInsanity();
                /* not grid - return */
                return;
        }

        /* malloc some arrays if first call and save continuum energies. */
        if( lgFirstRun )
        {
                long i1, i2;
                grid.numEnergies = rfield.nupper-2;
                grid.Energies = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(grid.numEnergies) );
                grid.Spectra = (realnum***)MALLOC(sizeof(realnum**)*(unsigned)(NUM_OUTPUT_TYPES) );

                for( i1=0; i1< NUM_OUTPUT_TYPES; i1++ )
                {
                        if( grid.lgOutputTypeOn[i1] )
                        {
                                grid.Spectra[i1] = (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(grid.totNumModels) );
                                for( i2=0; i2<grid.totNumModels; i2++ )
                                {
                                        grid.Spectra[i1][i2] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(grid.numEnergies) );
                                }
                        }
                }

                for( i1=0; i1<grid.numEnergies; i1++ )
                {
                        grid.Energies[i1] = rfield.AnuOrg[i1];
                }
                lgFirstRun = false;
        }

        ASSERT( lgFirstRun == false );
        ASSERT( optimize.nOptimiz>=0 && optimize.nOptimiz<grid.totNumModels);

        for( i=1; i< NUM_OUTPUT_TYPES; i++ )
        {
                /* Grab spectrum for xspec printout 
                 * at this point nOptimiz has already been incremented for first model */
                if( grid.lgOutputTypeOn[i] )
                        cdSPEC2( i, grid.numEnergies, grid.Spectra[i][optimize.nOptimiz]);
        }
        return;
}

---