/home66/gary/public_html/cloudy/c08_branch/source/parse_interp.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 */
/*ParseInterp parse parameters on interpolate command */
#include "cddefines.h"
#include "called.h"
#include "rfield.h"
#include "trace.h"
#include "input.h"
#include "parse.h"

void ParseInterp(char *chCard, 
  bool *lgEOF)
{
        char chLab4[5];
        bool lgDONE, 
          lgEOL,
          lgHit;
        long int i, 
          k,
          npairs;
        double cmax, 
          cmin, 
          fac;

        DEBUG_ENTRY( "ParseInterp()" );

        /* 
         * this sub reads in the "interpolate" command, and has
         * logic for the "continue" lines as well.
         * OUTPUT:
         * TNU is vector of energies where the grid is defined
         * TSLOP initially is vector of log fnu at each freq
         * converted into slopes here too
         */

        if( rfield.nspec >= LIMSPC-1 )
        {
                fprintf( ioQQQ, " Too many spectra entered.  Increase LIMSPC\n" );
                cdEXIT(EXIT_FAILURE);
        }
        /* >>chng 06 jul 16, fix memory leak when optimizing, PvH */
        if( !rfield.lgContMalloc[rfield.nspec] )
        {
                rfield.tNuRyd[rfield.nspec] = (realnum*)MALLOC((size_t)(NCELL*sizeof(realnum)) );
                rfield.tslop[rfield.nspec] = (realnum*)MALLOC((size_t)(NCELL*sizeof(realnum)) );
                rfield.tFluxLog[rfield.nspec] = (realnum*)MALLOC((size_t)(NCELL*sizeof(realnum)) );
                rfield.lgContMalloc[rfield.nspec] = true;
        }

        strcpy( rfield.chSpType[rfield.nspec], "INTER" );

        /* read all of the numbers on a line */
        npairs = 0;

        /* this is flag saying that all numbers are in */
        lgDONE = false;

        /* this flag says we hit end of command stream */
        *lgEOF = false;
        while( !lgDONE && !*lgEOF )
        {
                i = 5;
                lgEOL = false;

                /* keep scanning numbers until we hit eol for current line image */
                /* >>chng 01 aug 10, add check on npairs not exceeding NC_ELL as per 
                 * Ilse van Bemmel bug report */
                /*while( !lgEOL && npairs<rfield.nupper )*/
                while( !lgEOL && npairs<NCELL )
                {
                        rfield.tNuRyd[rfield.nspec][npairs] = 
                                (realnum)FFmtRead(chCard ,&i,INPUT_LINE_LENGTH,&lgEOL);
                        rfield.tFluxLog[rfield.nspec][npairs] = 
                                (realnum)FFmtRead(chCard ,&i,INPUT_LINE_LENGTH,&lgEOL);
                        ++npairs;
                }
                /* check if we fell through due to too many points, did not hit EOL */
                if( !lgEOL )
                {
                        fprintf( ioQQQ, "Too many continuum points were entered.\n" );
                        fprintf( ioQQQ, 
                                "The current logic limits the number of possible points to the value of NCELL, which is %i.\n",NCELL );
                        fprintf( ioQQQ, 
                                "Increase the value of NCELL in rfield.h.\nSorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                /* added too many above, so back off */
                --npairs;

                /* read a new line, checking for EOF */
                input_readarray(chCard,lgEOF);

                /* option to ignore all lines starting with #, *, or %. 
                 * >>chng 06 sep 04 use routine to check for comments */
                while( !*lgEOF && lgInputComment(chCard)/*((chCard[0] == '#' || chCard[0] == '*') || chCard[0] == '%')*/ )
                {
                        input_readarray(chCard,lgEOF);
                }

                /* get capital first four char of chCard */
                cap4( chLab4 , chCard );

                /* print the line, but only if it is a continue line */
                if( called.lgTalk && (strncmp(chLab4,"CONT",4)==0) )
                {
                        fprintf( ioQQQ, "                       * ");
                        /* next logic will print command plus spaces to right justified * */
                        k=0;
                        while( chCard[k]!='\0' )
                        {
                                fprintf(ioQQQ,"%c",chCard[k]);
                                ++k;
                        }
                        while( k<80 )
                        {
                                fprintf(ioQQQ,"%c",' ');
                                ++k;
                        }
                        fprintf( ioQQQ,"*\n"); 
                }

                /* now convert entire line to caps */
                caps(chCard);

                /* is this a continue line? */
                if( strncmp(chCard,"CONT",4) != 0 )
                {
                        /* we have a line but next command, not continue */
                        lgDONE = true;
                }

                /* this is another way to hit end of input stream - blank lines */
                if( chCard[0] == ' ' )
                        *lgEOF = true;
        }

        /* if valid next line, backup one line */
        if( lgDONE )
        {
                input.nRead -= input.iReadWay;
        }

        /* done reading all of the possible lines */
        --npairs;
        if( npairs < 1 )
        {
                fprintf( ioQQQ, "There must be at least 2 pairs to interpolate,\nSorry\n" );
                cdEXIT(EXIT_FAILURE);
        }

        /* >>chng 02 apr 19, from > to >=, test did not trigger due to overrun protection
         * in main loop */
        else if( npairs >= NCELL - 2 )
        {
                fprintf( ioQQQ, " Too many continuum points entered.\n" );
                fprintf( ioQQQ, 
                        "The current logic limits the number of possible points to the value of NCELL, which is %i.\nSorry.\n",NCELL );
                fprintf( ioQQQ, " Increase NCELL (in cddefines.h) to more than the number of continuum points.\n" );
                cdEXIT(EXIT_FAILURE);
        }

        if( rfield.tNuRyd[rfield.nspec][0] == 0. )
        {
                /* special case - if first energy is zero then it is low energy */
                if( rfield.tNuRyd[rfield.nspec][1] > 0. )
                {
                        /* second energy positive, assume linear Ryd */
                        rfield.tNuRyd[rfield.nspec][0] = (realnum)rfield.emm;
                }
                else if( rfield.tNuRyd[rfield.nspec][1] < 0. )
                {
                        /* second energy negative, assume log of Ryd */
                        rfield.tNuRyd[rfield.nspec][0] = (realnum)log10(rfield.emm);
                }
                else
                {
                        /* second energy zero, not allowed */
                        fprintf( ioQQQ, 
                                "An energy of zero was entered for element%3ld in INTERPOLATE and is not allowed.\nSorry\n", 
                          rfield.nspec );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        /* convert from log(Hz) to Ryd if first nu>5 */
        if( rfield.tNuRyd[rfield.nspec][0] >= 5. )
        {
                for( i=0; i < (npairs + 1); i++ )
                {
                        rfield.tNuRyd[rfield.nspec][i] = 
                                (realnum)pow(10.,rfield.tNuRyd[rfield.nspec][i] - 15.517);
                }
        }

        else if( rfield.tNuRyd[rfield.nspec][0] < 0. )
        {
                /* energies entered as logs */
                for( i=0; i < (npairs + 1); i++ )
                {
                        rfield.tNuRyd[rfield.nspec][i] = (realnum)pow(10.,(double)rfield.tNuRyd[rfield.nspec][i]);
                }
        }

        else
        {
                /* numbers are linear Rydbergs */
                for( i=0; i < (npairs + 1); i++ )
                {
                        if( rfield.tNuRyd[rfield.nspec][i] == 0. )
                        {
                                fprintf( ioQQQ, "An energy of zero was entered for element%3ld in INTERPOLATE and is not allowed.\nSorry\n", 
                                  i );
                                cdEXIT(EXIT_FAILURE);
                        }
                }
        }

        /* option to print debugging file then exit */
        {
                /* following should be set true to print information */
                enum {DEBUG_LOC=false};
                if( DEBUG_LOC  )
                {
                        for( i=0; i < npairs; i++ )
                        {
                                fprintf(ioQQQ,"%.4e\t%.3e\n",
                                        rfield.tNuRyd[rfield.nspec][i],
                                        pow(10.,(double)rfield.tFluxLog[rfield.nspec][i])  * rfield.tNuRyd[rfield.nspec][i]);
                        }
                        cdEXIT(EXIT_SUCCESS);
                }
        }

        for( i=0; i < npairs; i++ )
        {
                /* check that frequencies are monotonically increasing */
                if( rfield.tNuRyd[rfield.nspec][i+1] <= rfield.tNuRyd[rfield.nspec][i] )
                {
                        fprintf( ioQQQ, "The energies MUST be in increasing order.  Energy #%3ld=%10.2e Ryd was greater than or equal to the next one.\nSorry.\n", 
                          i, rfield.tNuRyd[rfield.nspec][i] );
                        cdEXIT(EXIT_FAILURE);
                }

                /* TFAC is energy, and TSLOP is slope in f_nu; not photons */
                rfield.tslop[rfield.nspec][i] = (realnum)((rfield.tFluxLog[rfield.nspec][i+1] - 
                  rfield.tFluxLog[rfield.nspec][i])/log10(rfield.tNuRyd[rfield.nspec][i+1]/
                  rfield.tNuRyd[rfield.nspec][i]));
        }

        /*>>chng 06 may 17, must insure that rNuRyd is sane through NCELL values, not just
         * nupper values.  This bit when stellar continuum had more cells than cloudy grid,
         * so npairs is much greater than nupper */
        /*if( npairs + 2 < rfield.nupper )*/
        if( npairs + 2 < NCELL )
        {
                /* zero out remainder of array */
                /*>>chng 06 may 17, same as above */
                /*for( i=npairs + 1; i < rfield.nupper; i++ )*/
                for( i=npairs + 1; i < NCELL; i++ )
                {
                        rfield.tNuRyd[rfield.nspec][i] = 0.;
                }
        }

        /* now check that array is defined over all energies */
        if( rfield.tNuRyd[rfield.nspec][0] > rfield.emm )
        {
                /* not defined over low energy part of array */
                fprintf( ioQQQ, 
                        "\n NOTE The incident continuum was not defined over the entire energy range. Some energies are set to zero.\n" );
                fprintf( ioQQQ, 
                        " NOTE You may be making a BIG mistake.\n\n" );
        }

        /* check on IR */
        if( rfield.tNuRyd[rfield.nspec][0] > rfield.emm )
                rfield.lgMMok = false;

        if( rfield.tNuRyd[rfield.nspec][0] > 1/36. )
                rfield.lgHPhtOK = false;

        /* gamma ray, EGAMRY is roughly 100MeV */
        if( rfield.tNuRyd[rfield.nspec][npairs] < rfield.egamry )
                rfield.lgGamrOK = false;

        /* EnerGammaRay is roughly 100keV; high is gamma ray */
        if( rfield.tNuRyd[rfield.nspec][npairs] < rfield.EnerGammaRay )
                rfield.lgXRayOK = false;

        /* find min and max of continuum */
        cmax = -38.;
        cmin = 28;

        /* rfield.tFluxLog can be very large or small */
        /* >>chng 01 jul 11, from <npairs to <=npairs, [npairs] is highest point in table */
        for( i=0; i <= npairs; i++ )
        {
                cmax = MAX2(cmax,rfield.tFluxLog[rfield.nspec][i]);
                cmin = MIN2(cmin,rfield.tFluxLog[rfield.nspec][i]);
        }

        /* check on dynamic range of input continuum */
        if( cmax - cmin > 74. )
        {
                fprintf( ioQQQ, "The dynamic range of the specified continuum is too large.\nSorry.\n" );
                cdEXIT(EXIT_FAILURE);
        }

        if( cmin < -37. )
        {
                fac = -cmin - 37.;
                /* >>chng 01 jul 11, from <npairs to <=npairs, [npairs] is highest point in table */
                for( i=0; i <= npairs; i++ )
                {
                        rfield.tFluxLog[rfield.nspec][i] += (realnum)fac;
                }
        }

        else if( cmax > 37. )
        {
                fac = cmax - 37.;
                /* >>chng 01 jul 11, from <npairs to <=npairs, [npairs] is highest point in table */
                for( i=0; i <= npairs; i++ )
                {
                        rfield.tFluxLog[rfield.nspec][i] -= (realnum)fac;
                }
        }

        /* option to print out results at this stage - "trace continuum" */
        if( trace.lgConBug && trace.lgTrace )
        {
                fprintf( ioQQQ, " Table for this continuum;\ni\tTNU\tTFAC\tTSLOP, npairs=%li\n",
                        npairs );
                for( i=0; i < npairs; i++ )
                {
                        fprintf( ioQQQ, "%li\t%.4e\t%.4e\t%.4e\n", 
                                i , rfield.tNuRyd[rfield.nspec][i], 
                          rfield.tFluxLog[rfield.nspec][i], rfield.tslop[rfield.nspec][i] );
                }
                i = npairs;
                fprintf( ioQQQ, "%li\t%.4e\t%.4e\n", 
                        i , rfield.tNuRyd[rfield.nspec][i], 
                        rfield.tFluxLog[rfield.nspec][i]);
        }

        /* renormalize continuum so that flux we will interpolate upon passes through unity
         * at near 1 Ryd.  but first we must find 1 Ryd in the array.
         * find 1 Ryd, npairs is one less than number of continuum pairs */
        i = 0;
        while( rfield.tNuRyd[rfield.nspec][i] <= 1. && i < npairs-1 )
        {
                i += 1;
        }
        /* i is now the table point where rfield.tNuRyd[i-1] is <= 1 ryd,
         * and rfield.tNuRyd[i] > 1 ryd */

        /* following is impossible but sanity check */
        if( i < 1 )
        {
                fprintf( ioQQQ, "ParseInput finds insane i after rfield.tNuRyd loop\n");
                ShowMe();
                cdEXIT(EXIT_FAILURE);
        }

        /* do the renormalization so 1 ryd is about unity */
        fac = rfield.tFluxLog[rfield.nspec][i-1];
        for( i=0; i <= npairs; i++ )
        {
                rfield.tFluxLog[rfield.nspec][i] -= (realnum)fac;
                /*fprintf(ioQQQ,"DEBUG parse %li %e \n", i , rfield.tFluxLog[rfield.nspec][i] );*/
        }

        /* finally check that we are within dynamic range of this cpu */
        cmin = log10( FLT_MIN );
        cmax = log10( FLT_MAX );
        lgHit = false;
        for( i=0; i <= npairs; i++ )
        {
                if( rfield.tFluxLog[rfield.nspec][i] <= cmin || 
                        rfield.tFluxLog[rfield.nspec][i] >= cmax )
                {
                        fprintf(ioQQQ,
                                " The log of the flux specified in interpolate pair %li is not within dynamic range of this CPU - please rescale.\n",i);
                        fprintf(ioQQQ,
                                " The frequency is %f and the log of the flux is %f.\n\n",
                                rfield.tNuRyd[rfield.nspec][i] , 
                                rfield.tFluxLog[rfield.nspec][i]);
                        lgHit = true;
                }
        }
        if( lgHit )
        {
                fprintf(ioQQQ,"\n NOTE The log of the flux given in an interpolate command is outside the range of this cpu.\n");
                fprintf(ioQQQ," NOTE I will try to renormalize it to be within the range of this cpu, but if I crash, this is a likely reason.\n");
                fprintf(ioQQQ," NOTE Note that the interpolate command only is used for the shape of the continuum.\n");
                fprintf(ioQQQ," NOTE The order of magnitude of the flux is not used in any way.\n");
                fprintf(ioQQQ," NOTE For safety this could be of order unity.\n\n");
        }

        /* zero out remainder of array */
        for( i=npairs+1; i < rfield.nupper; i++ )
        {
                /* >>chng 05 jul 27, next three indices had been ipspec, changed to nspec
                 * bug caught by Ralf Quast */
                rfield.tFluxLog[rfield.nspec][i] = 0.;
                rfield.tslop[rfield.nspec][i] = 0.;
                rfield.tNuRyd[rfield.nspec][i] = 0.;
        }

        /* now increment number of continua */
        ++rfield.nspec;
        if( rfield.nspec >= LIMSPC )
        {
                /* too many continua were entered */
                fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
                cdEXIT(EXIT_FAILURE);
        }

        return;
}

---