cont_createmesh.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 */
/*ContCreateMesh calls fill to set up continuum energy mesh if first call, 
 * otherwise reset to original mesh */
/*fill define the continuum energy grid over a specified range */
/*ChckFill perform sanity check confirming that the energy array has been properly filled */
/*rfield_opac_malloc MALLOC space for opacity arrays */
/*read_continuum_mesh read the continuum definition from the file continuum_mesh.ini */
#include "cddefines.h"
#include "rfield.h"
#include "iterations.h"
#include "physconst.h"
#include "dense.h"
#include "trace.h"
#include "opacity.h"
#include "ipoint.h"
#include "geometry.h"
#include "continuum.h"

/* read the continuum definition from the file continuum_mesh.ini */
STATIC void read_continuum_mesh( void );

/*fill define the continuum energy grid over a specified range */
STATIC void fill(double fenlo, 
  double fenhi, 
  double resolv, 
  long int *n0, 
  long int *ipnt,
  /* this says only count, do not fill */
  bool lgCount );

/*rfield_opac_malloc MALLOC space for opacity arrays */
STATIC void rfield_opac_malloc(void);

/*ChckFill perform sanity check confirming that the energy array has been properly filled */
STATIC void ChckFill(void);

void ContCreateMesh(void)
{
        long int 
          i, 
          ipnt, 
          n0;

        /* flag to say whether pointers have ever been evaluated */
        static bool lgPntEval = false;

        DEBUG_ENTRY( "ContCreateMesh()" );

        /* lgPntEval is local static variable defined false when defined. 
         * it is set true below, so that pointers only created one time in the
         * history of this coreload. */
        if( lgPntEval )
        {
                if( trace.lgTrace )
                {
                        fprintf( ioQQQ, " ContCreateMesh called, not evaluating.\n" );
                }
                /* now save current form of energy array */
                for( i=0; i < rfield.nupper; i++ )
                {
                        rfield.anu[i] = rfield.AnuOrg[i];
                        rfield.anu2[i] = rfield.anu[i]*rfield.anu[i];
                }
                return;
        }
        else
        {
                if( trace.lgTrace )
                {
                        fprintf( ioQQQ, " ContCreateMesh called first time.\n" );
                }
                lgPntEval = true;
        }

        /* set size of arrays that continuum iteration information */

        /* read in the continuum mesh resolution definition */
        /* >>chng 01 sep 29, add external file "continuum_mesh.ini" with fill parameters */
        read_continuum_mesh();

        /* fill in continuum with freq points
         * arg are range pointer, 2 energy limits, resolution
         * first argument is lower energy of the continuum range to be defined
         * second argument is upper range; both are in Rydbergs
         * third number is the relative energy resolution, dnu/nu,
         * for that range of the continuum
         * last two numbers are internal book keeping
         * N0 is number of energy cells used so far
         * IPNT is a counter for the number of fills used
         *
         * if this is changed, then also change warning in GetTable about using
         * transmitted continuum - it says version number where continuum changed
         * */
        n0 = 1;
        ipnt = 0;
        /* this is number of ranges that will be introduced below*/
        continuum.nrange = 0;

        /* ================================================================ */
        /* NB - this block must agree exactly with the one that follows     */
        n0 = 1;
        ipnt = 0;
        /* this is number of ranges that will be introduced below*/
        continuum.nrange = 0;
        continuum.StoredEnergy[continuum.nStoredBands-1] = rfield.egamry;

        fill(rfield.emm, continuum.StoredEnergy[0] , continuum.StoredResolution[0],&n0,&ipnt,true );
        for(i=1; i<continuum.nStoredBands; ++i )
        {
                fill(continuum.StoredEnergy[i-1]      , 
                        continuum.StoredEnergy[i],  
                        continuum.StoredResolution[i],
                        &n0,&ipnt,true);
        }
        /* ================================================================ */

        /* at this point debugger shows that anu and widflx are defined 
         * up through n0-2 - due to c offset -1 from Fortran! */
        rfield.nupper = n0 - 1;
        /* there must be a cell above nflux for us to pass unity through the vol integrator */
        if( rfield.nupper >= NCELL )
        {
                fprintf(ioQQQ," Currently the arrays that hold interpolated tables can only hold %i points.\n",NCELL);
                fprintf(ioQQQ," This continuum mesh really needs to have %li points.\n",rfield.nupper);
                fprintf(ioQQQ," Please increase the value of NCELL in rfield.h and recompile.\n Sorry.");
                cdEXIT(EXIT_FAILURE);
        }

        /*>>chng 04 oct 10, from nflux = nupper to nupper-1 since vectors must malloc to nupper, but
         * will address [nflux] for unit continuum test */
        rfield.nflux = rfield.nupper-1;

        /* allocate space for continuum arrays within rfield.h and opacity arrays in opacity.h
         * sets lgRfieldMalloced true */
        rfield_opac_malloc();

        /* geometry.nend_max is largest number of zones needed on any iteration,
         * will use it to malloc arrays that save source function as function of zone */
        /* now change all limits, for all iterations, to this value */
        geometry.nend_max = geometry.nend[0];
        for( i=1; i < iterations.iter_malloc; i++ )
        {
                geometry.nend_max = MAX2( geometry.nend_max , geometry.nend[i] );
        }
        /* nend_max+1 because search phase is zone 0, first zone at illumin face is 1 */
        rfield.ConEmitLocal = (realnum**)MALLOC( (size_t)(geometry.nend_max+1)*sizeof(realnum *)  );
        rfield.ConSourceFcnLocal = (realnum**)MALLOC( (size_t)(geometry.nend_max+1)*sizeof(realnum *)  );
        for( i=0;i<(geometry.nend_max+1); ++i )
        {
                rfield.ConEmitLocal[i] = (realnum*)MALLOC( (size_t)rfield.nupper*sizeof(realnum)  );
                rfield.ConSourceFcnLocal[i] = (realnum*)MALLOC( (size_t)rfield.nupper*sizeof(realnum)  );
        }
        for( i=0;i<(geometry.nend_max+1); ++i )
        {
                for( long j=0; j<rfield.nupper; ++j)
                {
                        rfield.ConSourceFcnLocal[i][j] = 1.;
                }
        }

        /* ================================================================ */
        n0 = 1;
        ipnt = 0;

        /* this is number of ranges that will be introduced below*/
        continuum.nrange = 0;

        /* the default array values are set in continuum_mesh.ini */
        fill(rfield.emm, continuum.StoredEnergy[0] , continuum.StoredResolution[0] , 
                &n0,&ipnt,false);
        for(i=1; i<continuum.nStoredBands; ++i )
        {
                fill(continuum.StoredEnergy[i-1]      , 
                        continuum.StoredEnergy[i],  
                        continuum.StoredResolution[i],
                        &n0,&ipnt,false);
        }

        /* ================================================================ */

        /* fill in the false highest cell used for unit verification */
        rfield.widflx[rfield.nupper] = rfield.widflx[rfield.nupper-1];
        rfield.anu[rfield.nupper] = rfield.anu[rfield.nupper-1] + 
                rfield.widflx[rfield.nupper];

        /* there must be a cell above nflux for us to pass unity through the vol integrator
         * as a sanity check.  assert that this is true so we will crash if ever changed 
        ASSERT( rfield.nupper +1 <= rfield.nupper );*/

        /* this is done here when the space is first allocated, 
         * then done on every subsequent initialization in zero.c */
        rfield_opac_zero( 0 , rfield.nupper );

        /* this is a sanity check for results produced above by fill */
        ChckFill();

        /* now fix widflx array so that it is correct */
        for( i=1; i<rfield.nupper-1; ++i )
        {
                rfield.widflx[i] = ((rfield.anu[i+1] - rfield.anu[i]) + (rfield.anu[i] - 
                        rfield.anu[i-1]))/2.f;
        }

        ipnt = 0;
        /* now save current form of array, and define some quantities related to it */
        for( i=0; i < rfield.nupper; i++ )
        {
                double alf , bet;

                rfield.AnuOrg[i] = rfield.anu[i];
                rfield.anusqr[i] = (realnum)sqrt(rfield.AnuOrg[i]);
                /* following are Compton exchange factors from Tarter */
                /* this code also appears in highen, but coef needed before that routine called. */
                alf = 1./(1. + rfield.anu[i]*(1.1792e-4 + 7.084e-10*rfield.anu[i]));
                bet = 1. - alf*rfield.anu[i]*(1.1792e-4 + 2.*7.084e-10*rfield.anu[i])/4.;
                rfield.csigh[i] = (realnum)(alf*rfield.anu[i]*rfield.anu[i]*3.858e-25);
                rfield.csigc[i] = (realnum)(alf*bet*rfield.anu[i]*3.858e-25);
                rfield.anu2[i] = rfield.anu[i]*rfield.anu[i];

                /* >>chng 05 feb 28, add transmission and mapping coef */
                /* map these coarse continua into fine continuum grid */
                if( rfield.anu[i] < rfield.fine_ener_lo || rfield.anu[i]>rfield.fine_ener_hi )
                {
                        /* 0 (false) says not defined */
                        rfield.ipnt_coarse_2_fine[i] = 0;
                }
                else
                {
                        if( ipnt==0 )
                        {
                                /* this is the first one that maps onto the fine array */
                                rfield.ipnt_coarse_2_fine[i] = 0;
                                ipnt = 1;
                        }
                        else
                        {
                                /* find first fine frequency that is greater than this coarse value */
                                while (ipnt < rfield.nfine_malloc && rfield.fine_anu[ipnt]<rfield.anu[i] )
                                {
                                        ++ipnt;
                                }
                                rfield.ipnt_coarse_2_fine[i] = ipnt;
                        }
                }
                /*fprintf(ioQQQ," coarse %li nu= %.3e points to fine %li nu=%.3e\n",
                        i, rfield.anu[i] , rfield.ipnt_coarse_2_fine[i] , rfield.fine_anu[rfield.ipnt_coarse_2_fine[i]] );*/
        }
        rfield.resetCoarseTransCoef();
        return;
}

/*fill define the continuum energy grid over a specified range, called by ContCreateMesh */
STATIC void fill(
  /* lower bounds to this energy range */
  double fenlo, 
  /* upper bounds to this continuum range */
  double fenhi, 
  /* relative energy resolution */
  double resolv, 
  /* starting index within frequency grid */
  long int *n0, 
  /* which energy band this is */
  long int *ipnt,  
  /* this says only count, do not fill */
  bool lgCount )
{
        long int i, 
          nbin;
        realnum widtot;
        double aaa , bbb;

        DEBUG_ENTRY( "fill()" );

        ASSERT( fenlo>0. && fenhi>0. && resolv>0. );

        /* this is the number of cells needed to fill the array with numbers at the requested resolution */
        nbin = (long int)(log(10.)*log10(fenhi/fenlo)/resolv + 1);

        if( lgCount )
        {
                /* true means only count number of cells, don't do anything */
                *n0 += nbin;
                return;
        }

        if( *ipnt > 0 && fabs(1.-fenlo/continuum.filbnd[*ipnt]) > 1e-4 )
        {
                fprintf( ioQQQ, " FILL improper bounds.\n" );
                fprintf( ioQQQ, " ipnt=%3ld fenlo=%11.4e filbnd(ipnt)=%11.4e\n", 
                  *ipnt, fenlo, continuum.filbnd[*ipnt] );
                cdEXIT(EXIT_FAILURE);
        }

        ASSERT( *ipnt < continuum.nStoredBands );

        continuum.ifill0[*ipnt] = *n0 - 1;
        continuum.filbnd[*ipnt] = (realnum)fenlo;
        continuum.filbnd[*ipnt+1] = (realnum)fenhi;

        /* this is the number of cells needed to fill the array with numbers 
        nbin = (long int)(log(10.)*log10(fenhi/fenlo)/resolv + 1);*/
        continuum.fildel[*ipnt] = (realnum)(log10(fenhi/fenlo)/nbin);

        if( continuum.fildel[*ipnt] < 0.01 )
        {
                continuum.filres[*ipnt] = (realnum)(log(10.)*continuum.fildel[*ipnt]);
        }
        else
        {
                continuum.filres[*ipnt] = (realnum)((pow(10.,2.*continuum.fildel[*ipnt]) - 1.)/2./
                        pow((realnum)10.f,continuum.fildel[*ipnt]));
        }

        if( (*n0 + nbin-2) > rfield.nupper )
        {
                fprintf( ioQQQ, " Fill would need %ld cells to get to an energy of %.3e\n", 
                  *n0 + nbin, fenhi );
                fprintf( ioQQQ, " This is a major logical error in fill.\n");
                ShowMe();
                cdEXIT(EXIT_FAILURE);
        }

        widtot = 0.;
        for( i=0; i < nbin; i++ )
        {
                bbb = continuum.fildel[*ipnt]*((realnum)(i) + 0.5);
                aaa = pow( 10. , bbb );

                rfield.anu[i+continuum.ifill0[*ipnt]] = (realnum)(fenlo*aaa);

                rfield.widflx[i+continuum.ifill0[*ipnt]] = rfield.anu[i+continuum.ifill0[*ipnt]]*
                  continuum.filres[*ipnt];

                widtot += rfield.widflx[i+continuum.ifill0[*ipnt]];
        }

        *n0 += nbin;
        if( trace.lgTrace && (trace.lgConBug || trace.lgPtrace) )
        {
                fprintf( ioQQQ, 
                        " FILL range%2ld from%10.3e to%10.3eR in%4ld cell; ener res=%10.3e WIDTOT=%10.3e\n", 
                  *ipnt, 
                  rfield.anu[continuum.ifill0[*ipnt]] - rfield.widflx[continuum.ifill0[*ipnt]]/2., 
                  rfield.anu[continuum.ifill0[*ipnt]+nbin-1] + rfield.widflx[continuum.ifill0[*ipnt]+nbin-1]/2., 
                  nbin, 
                  continuum.filres[*ipnt], 
                  widtot );

                fprintf( ioQQQ, " The requested range was%10.3e%10.3e The requested resolution was%10.3e\n", 
                  fenlo, fenhi, resolv );
        }

        /* nrange is number of ranges  */
        *ipnt += 1;
        continuum.nrange = MAX2(continuum.nrange,*ipnt);
        return;
}

/*ChckFill perform sanity check confirming that the energy array has been properly filled */
STATIC void ChckFill(void)
{
        bool lgFail;
        long int i, 
          ipnt;
        double energy;

        DEBUG_ENTRY( "ChckFill()" );

        ASSERT( rfield.anu[0] >= rfield.emm*0.99 );
        ASSERT( rfield.anu[rfield.nupper-1] <= rfield.egamry*1.01 );

        lgFail = false;
        for( i=0; i < continuum.nrange; i++ )
        {
                /* test middle of energy bound */
                energy = (continuum.filbnd[i] + continuum.filbnd[i+1])/2.;
                ipnt = ipoint(energy);
                if( energy < rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5 )
                {
                        fprintf( ioQQQ, " ChckFill middle test low fail\n" );
                        lgFail = true;
                }

                /* >>chng 02 jul 16, add second test - when "set resol 10" used,
                 * very large values of cell width, combined with fact that cells
                 * are log increasing, causes problem.  */
                else if( (energy > rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]*0.5) &&
                        ( energy > rfield.anu[ipnt] - rfield.widflx[ipnt]*0.5 ) )
                {
                        fprintf( ioQQQ, " ChckFill middle test high fail\n" );
                        lgFail = true;
                }

                /* test near low bound */
                energy = continuum.filbnd[i]*0.99 + continuum.filbnd[i+1]*0.01;
                ipnt = ipoint(energy);
                if( energy < rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5 )
                {
                        fprintf( ioQQQ, " ChckFill low test low fail\n" );
                        lgFail = true;
                }

                else if( energy > rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]* 0.5 )
                {
                        fprintf( ioQQQ, " ChckFill low test high fail\n" );
                        lgFail = true;
                }

                /* test near high bound */
                energy = continuum.filbnd[i]*0.01 + continuum.filbnd[i+1]*0.99;
                ipnt = ipoint(energy);

                if( energy < rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5 )
                {
                        fprintf( ioQQQ, " ChckFill high test low fail\n" );
                        lgFail = true;
                }
                /* >>chng 02 jul 16, add second test - when "set resol 10" used,
                 * very large values of cell width, combined with fact that cells
                 * are log increasing, causes problem.  */
                else if( (energy > rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]*0.5) &&
                        ( energy > rfield.anu[ipnt] - rfield.widflx[ipnt]*0.5 ) )
                {
                        fprintf( ioQQQ, " ChckFill high test high fail\n" );
                        lgFail = true;
                }
        }

        if( lgFail )
        {
                cdEXIT(EXIT_FAILURE);
        }
        return;
}

/* MALLOC arrays within rfield */
STATIC void rfield_opac_malloc(void)
{
        long i;

        DEBUG_ENTRY( "rfield_opac_malloc()" );

        /* allocate one more than we use for the unit integration,
         * will back up at end of routine */
        ++rfield.nupper;

        /* >>chng 03 feb 12, add fine mesh fine grid fine opacity array to keep track of line overlap */
        /* frequency range in Rydberg needed for all resonance lines */
        rfield.fine_ener_lo = rfield.emm;
        rfield.fine_ener_hi = 1500.f;

        /* set resolution of fine continuum mesh. 
         * rfield.fine_opac_velocity_width is width per cell, cm/s 
         * choose width so that most massive species (usually Fe) is well resolved
         * 
         * rfield.fine_opac_nelem is the most massive (hence sharpest line)
         * we will worry about.  By default this is iron but can be changed
         * with SET FINE CONTINUUM command 
         * 
         * TeLowestFineOpacity of 1e4 K is temperature were line width is 
         * evaluated.  Tests were done using the stop temperature in its place
         * Te below 1e4 K made fine opacity grid huge 
         * do not let temp get higher than 1e4 either - code run with stop temp 10 set
         * stop temp of 1e10K and assert thrown at line 204 of cont_createpointers.c 
         * simply use 1e4 K as a characteristic temperature */
        double TeLowestFineOpacity = 1e4;
        rfield.fine_opac_velocity_width = 
                (realnum)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*TeLowestFineOpacity/
                dense.AtomicWeight[rfield.fine_opac_nelem] ) / 
                /* we want fine_opac_nresolv continuum elements across this line
                 * default is 1, changed with SET FINE CONTINUUM command */
                rfield.fine_opac_nresolv;

        /* we are at first zone so velocity shift is zero */
        rfield.ipFineConVelShift = 0;

        /* dimensionless resolution, dE/E, this is used in ipoint to get offset in find mesh */
        rfield.fine_resol = rfield.fine_opac_velocity_width / SPEEDLIGHT;

        /* the number of cells needed */
        rfield.nfine_malloc = (long)(log10( rfield.fine_ener_hi / rfield.fine_ener_lo ) / log10( 1. + rfield.fine_resol ) );
        if( rfield.nfine_malloc <= 0 )
                TotalInsanity();
        rfield.nfine = rfield.nfine_malloc;

        /* this is the fine opacity array to ghost the main low-resolution array */
        rfield.fine_opac_zone = (realnum *)MALLOC(sizeof(realnum)*(unsigned)rfield.nfine_malloc );
        memset(rfield.fine_opac_zone , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) );

        /* this is the fine total optical array to ghost the main low-resolution array */
        rfield.fine_opt_depth = (realnum *)MALLOC(sizeof(realnum)*(unsigned)rfield.nfine_malloc );
        memset(rfield.fine_opt_depth , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) );

        rfield.fine_anu = (realnum *)MALLOC(sizeof(realnum)*(unsigned)rfield.nfine_malloc );

        /* now fill in energy array */
        ASSERT( rfield.fine_ener_lo > 0. && rfield.fine_resol > 0 );
        for( i=0;i<rfield.nfine_malloc; ++i )
        {
                rfield.fine_anu[i] = rfield.fine_ener_lo * (realnum)pow( (1.+rfield.fine_resol), (i+1.) );
        }
        /* done with fine array */

        /* used to count number of lines per cell */
        rfield.line_count = (long *)MALLOC(sizeof(long)*(unsigned)NCELL );
        for( i=0; i<rfield.nupper; ++i)
        {
                rfield.line_count[i] = 0;
        }
        rfield.anu = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
        rfield.AnuOrg = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
        rfield.widflx = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.anulog = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.anusqr = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.anu2 = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.anu3 = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.flux_beam_time = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.flux_isotropic = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.flux_beam_const = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.flux_accum = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.ExtinguishFactor = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.convoc = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.OccNumbBremsCont = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.OccNumbIncidCont = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.OccNumbDiffCont = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.OccNumbContEmitOut = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.ConInterOut = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.SummedCon = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
        rfield.SummedDif = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.SummedDifSave = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.SummedOcc = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.ConOTS_local_photons = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.DiffuseEscape = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.TotDiff2Pht = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.ConOTS_local_OTS_rate = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.otslin = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.otscon = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.outlin_noplot = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.flux_beam_const_save = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.flux_time_beam_save = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.flux_isotropic_save = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.setCoarseTransCoefPtr((realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) ));
        rfield.DiffuseLineEmission = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.ipnt_coarse_2_fine = (long int*)MALLOC((size_t)(rfield.nupper*sizeof(long int)) );

        /* possibly save cumulative flux */
        rfield.flux = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
        rfield.ConEmitReflec = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
        rfield.ConEmitOut = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
        rfield.ConRefIncid = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
        rfield.flux_total_incident = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
        rfield.reflin = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );
        rfield.outlin = (realnum**)MALLOC((size_t)(2*sizeof(realnum*)) );

        for( i=0; i<2; ++i )
        {
                rfield.flux[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
                rfield.ConEmitReflec[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
                rfield.ConEmitOut[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
                rfield.ConRefIncid[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
                rfield.flux_total_incident[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
                rfield.reflin[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
                rfield.outlin[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        }
        // the cumulative (time integral) emission
        memset(rfield.flux[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
        memset(rfield.ConEmitReflec[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
        memset(rfield.ConEmitOut[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
        memset(rfield.ConRefIncid[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
        memset(rfield.flux_total_incident[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
        memset(rfield.reflin[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );
        memset(rfield.outlin[1] , 0 , (unsigned long)rfield.nupper*sizeof(realnum) );

        /* chng 02 may 16, by Ryan...added array for gaunt factors for ALL charges, malloc here.        */
        /* First index is EFFECTIVE CHARGE MINUS ONE!   */
        rfield.gff = (realnum**)MALLOC((size_t)((LIMELM+1)*sizeof(realnum*)) );
        for( i = 1; i <= LIMELM; i++ )
        {
                rfield.gff[i] = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        }

        rfield.csigh = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        rfield.csigc = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );

        rfield.comdn = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
        rfield.comup = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
        rfield.ContBoltz = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );

        /*realnum rfield.otssav[NC_ELL][2];*/
        rfield.otssav = (realnum**)MALLOC((size_t)(rfield.nupper*sizeof(realnum*)));
        for( i=0; i<rfield.nupper; ++i)
        {
                rfield.otssav[i] = (realnum*)MALLOC(2*sizeof(realnum));
        }


        /* char rfield.chLineLabel[NLINES][5];*/
        rfield.chLineLabel = (char**)MALLOC((size_t)(rfield.nupper*sizeof(char*)));
        rfield.chContLabel = (char**)MALLOC((size_t)(rfield.nupper*sizeof(char*)));

        /* now allocate all the labels for each of the above lines */
        for( i=0; i<rfield.nupper; ++i)
        {
                rfield.chLineLabel[i] = (char*)MALLOC(5*sizeof(char));
                rfield.chContLabel[i] = (char*)MALLOC(5*sizeof(char));
        }

        opac.TauAbsFace = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        memset( opac.TauAbsFace , 0 , rfield.nupper*sizeof(realnum) );

        opac.TauScatFace = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        opac.E2TauAbsFace = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        opac.E2TauAbsTotal = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        opac.TauAbsTotal = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        opac.E2TauAbsOut = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        opac.ExpmTau = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );
        opac.tmn = (realnum*)MALLOC((size_t)(rfield.nupper*sizeof(realnum)) );

        opac.opacity_abs = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
        opac.opacity_abs_savzon1 = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
        opac.OldOpacSave = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
        opac.opacity_sct = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
        opac.albedo = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
        opac.opacity_sct_savzon1 = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
        opac.OpacStatic = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
        opac.FreeFreeOpacity = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );
        opac.ExpZone = (double*)MALLOC((size_t)(rfield.nupper*sizeof(double)) );

        opac.TauAbsGeo = (realnum**)MALLOC((size_t)(2*sizeof(realnum *)) );
        opac.TauScatGeo = (realnum**)MALLOC((size_t)(2*sizeof(realnum *)) );
        opac.TauTotalGeo = (realnum**)MALLOC((size_t)(2*sizeof(realnum *)) );

        for( i=0; i<2; ++i)
        {
                opac.TauAbsGeo[i] = (realnum*)MALLOC(rfield.nupper*sizeof(realnum));
                opac.TauScatGeo[i] = (realnum*)MALLOC(rfield.nupper*sizeof(realnum));
                opac.TauTotalGeo[i] = (realnum*)MALLOC(rfield.nupper*sizeof(realnum));
        }

        /* fix allocate trick for one more than we use for the unit integration */
        --rfield.nupper;

        /* say that space exists */
        lgRfieldMalloced = true;
        return;
}


/* read the continuum definition from the file continuum_mesh.ini */
STATIC void read_continuum_mesh( void )
{
        FILE *ioDATA;
        char chLine[INPUT_LINE_LENGTH];
        long i;
        bool lgEOL;
        long i1 , i2 , i3;

        DEBUG_ENTRY( "read_continuum_mesh()" );

        if( trace.lgTrace )
                fprintf( ioQQQ," read_continuum_mesh opening continuum_mesh.ini:");

        ioDATA = open_data( "continuum_mesh.ini", "r" );

        /* first line is a version number and does not count */
        if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
        {
                fprintf( ioQQQ, " read_continuum_mesh could not read first line of continuum_mesh.ini.\n");
                cdEXIT(EXIT_FAILURE);
        }
        /* count how many lines are in the file, ignoring all lines
         * starting with '#' */
        continuum.nStoredBands = 0;
        while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
        {
                /* we want to count the lines that do not start with #
                 * since these contain data */
                if( chLine[0] != '#')
                        ++continuum.nStoredBands;
        }

        /* we now have number of lines containing pairs of bounds,
         * allocate space for the arrays we will need */
        continuum.filbnd = 
                ((realnum *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(realnum )));
        continuum.fildel = 
                ((realnum *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(realnum )));
        continuum.filres = 
                ((realnum *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(realnum )));
        continuum.ifill0 = 
                ((long *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(long )));
        continuum.StoredEnergy = 
                ((double *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(double )));
        continuum.StoredResolution = 
                ((double *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(double )));

        /* now rewind the file so we can read it a second time*/
        if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
        {
                fprintf( ioQQQ, " read_continuum_mesh could not rewind continuum_mesh.ini.\n");
                cdEXIT(EXIT_FAILURE);
        }

        /* check that magic number is ok */
        if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
        {
                fprintf( ioQQQ, " read_continuum_mesh could not read first line of continuum_mesh.ini.\n");
                cdEXIT(EXIT_FAILURE);
        }

        i = 1;
        /* continuum mesh magic number */
        i1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
        i2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
        i3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);

        bool lgResPower;

        /* the following is the set of numbers that appear at the start of continuum_mesh.ini */
        if( i1 == 1 && i2 == 9 && i3 == 29 )
                // old version of the file (c08 and older), this has pairs: upper limit freq range, resolution
                // this format is still supported to accomodate users with existing continuum_mesh.ini files.
                lgResPower = false;
        else if( i1 == 10 && i2 == 8 && i3 == 8 )
                // new version of the file (c10 and newer), this has pairs: upper limit freq range, resolving power
                // resolving power = 1./resolution
                lgResPower = true;
        else
        {
                fprintf( ioQQQ, 
                        " read_continuum_mesh: the version of continuum_mesh.ini is not supported.\n" );
                fprintf( ioQQQ, 
                        " I found version number %li %li %li.\n" ,
                         i1 , i2 , i3 );
                fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
                cdEXIT(EXIT_FAILURE);
        }

        /* this starts at 1 not 0 since zero is reserved for the
         * dummy line */
        continuum.nStoredBands = 0;
        while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
        {
                /* only look at lines without '#' in first col */
                if( chLine[0] != '#')
                {
                        i = 1;
                        continuum.StoredEnergy[continuum.nStoredBands] = FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
                        continuum.StoredResolution[continuum.nStoredBands] = FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);

                        // continuum energy could be 0 to indicate low or high energy bounds of code
                        // but none can be negative
                        if( continuum.StoredEnergy[continuum.nStoredBands]<0. ||
                                        continuum.StoredResolution[continuum.nStoredBands]<=0. )
                        {
                                fprintf(ioQQQ, "DISASTER PROBLEM continuum_mesh.ini has a non-positive number.\n");
                                cdEXIT(EXIT_FAILURE);
                        }

                        // convert resolving power (entered quantity) into resolution
                        if( lgResPower )
                                continuum.StoredResolution[continuum.nStoredBands] = 1./
                                        continuum.StoredResolution[continuum.nStoredBands];

                        /* this is option to rescale resolution with set resolution command */
                        continuum.StoredResolution[continuum.nStoredBands] *= continuum.ResolutionScaleFactor;

                        ++continuum.nStoredBands;
                }
        }

        fclose( ioDATA );

        /* now verify continuum grid is ok - first are all values but the last positive? */
        for( i=1; i<continuum.nStoredBands-1; ++i )
        {
                if( continuum.StoredEnergy[i-1]>=continuum.StoredEnergy[i] )
                {
                        fprintf( ioQQQ, 
                                " read_continuum_mesh: The continuum definition array energies must be in increasing order.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }
        if( continuum.StoredEnergy[continuum.nStoredBands-1]!=0 )
        {
                fprintf( ioQQQ, 
                        " read_continuum_mesh: The last continuum array energies must be zero.\n" );
                cdEXIT(EXIT_FAILURE);
        }
        return;
}

/*rfield_opac_zero zero out rfield arrays between certain limits */
void rfield_opac_zero( 
                                          /* index for first element in arrays to be set to zero */
                                          long lo , 
                                          /* array index for highest element to be set */
                                          long ihi )
{
        long int i;

        /* >>chng 01 aug 19, space not allocated yet,
        * following code must also be present in contcreatemesh where
        * space allocated for the first time */
        if( lgRfieldMalloced )
        {
                unsigned long n=(unsigned long)(ihi-lo+1);
                memset(&rfield.OccNumbDiffCont[lo]      , 0 , n*sizeof(realnum) );
                memset(&rfield.OccNumbContEmitOut[lo]   , 0 , n*sizeof(realnum) );
                memset(&rfield.ContBoltz[lo]            , 0 , n*sizeof(double) );
                /*>>chng 06 aug 15, this is now 2D array, saving diffuse continuum
                * over all zones for use in exact RT */
                /*memset(&rfield.ConEmitLocal[lo]         , 0 , n*sizeof(realnum) );*/
                memset(&rfield.ConEmitReflec[0][lo]        , 0 , n*sizeof(realnum) );
                memset(&rfield.ConEmitOut[0][lo]           , 0 , n*sizeof(realnum) );
                memset(&rfield.reflin[0][lo]               , 0 , n*sizeof(realnum) );
                memset(&rfield.ConRefIncid[0][lo]          , 0 , n*sizeof(realnum) );
                memset(&rfield.SummedCon[lo]            , 0 , n*sizeof(double) );
                memset(&rfield.OccNumbBremsCont[lo]     , 0 , n*sizeof(realnum) );
                memset(&rfield.convoc[lo]               , 0 , n*sizeof(realnum) );
                memset(&rfield.flux[0][lo]                 , 0 , n*sizeof(realnum) );
                memset(&rfield.flux_total_incident[0][lo]  , 0 , n*sizeof(realnum) );
                memset(&rfield.flux_beam_const_save[lo] , 0 , n*sizeof(realnum) );
                memset(&rfield.flux_time_beam_save[lo]  , 0 , n*sizeof(realnum) );
                memset(&rfield.flux_isotropic_save[lo]  , 0 , n*sizeof(realnum) );
                memset(&rfield.SummedOcc[lo]            , 0 , n*sizeof(realnum) );
                memset(&rfield.SummedDif[lo]            , 0 , n*sizeof(realnum) );
                memset(&rfield.flux_accum[lo]           , 0 , n*sizeof(realnum) );
                memset(&rfield.otslin[lo]               , 0 , n*sizeof(realnum) );
                memset(&rfield.otscon[lo]               , 0 , n*sizeof(realnum) );
                memset(&rfield.ConInterOut[lo]          , 0 , n*sizeof(realnum) );
                memset(&rfield.outlin[0][lo]               , 0 , n*sizeof(realnum) );
                memset(&rfield.outlin_noplot[lo]        , 0 , n*sizeof(realnum) );
                memset(&rfield.ConOTS_local_OTS_rate[lo], 0 , n*sizeof(realnum) );
                memset(&rfield.ConOTS_local_photons[lo] , 0 , n*sizeof(realnum) );
                memset(&opac.OldOpacSave[lo]            , 0 , n*sizeof(double) );
                memset(&opac.opacity_abs[lo]            , 0 , n*sizeof(double) );
                memset(&opac.opacity_sct[lo]            , 0 , n*sizeof(double) );
                memset(&opac.albedo[lo]                 , 0 , n*sizeof(double) );
                memset(&opac.FreeFreeOpacity[lo]        , 0 , n*sizeof(double) );

                /* these are not defined on first iteration */
                memset( &opac.E2TauAbsTotal[lo]        , 0 , n*sizeof(realnum) );
                memset( &opac.E2TauAbsOut[lo]          , 0 , n*sizeof(realnum) );
                memset( &opac.TauAbsTotal[lo]          , 0 , n*sizeof(realnum) );

                for( i=lo; i <= ihi; i++ )
                {
                        opac.TauTotalGeo[0][i] = opac.taumin;
                        opac.TauAbsGeo[0][i] = opac.taumin;
                        opac.TauScatGeo[0][i] = opac.taumin;
                        opac.tmn[i] = 1.;
                        opac.ExpZone[i] = 1.;
                        opac.E2TauAbsFace[i] = 1.;
                        opac.ExpmTau[i] = 1.;
                        opac.OpacStatic[i] = 1.;
                }
                /* also zero out fine opacity fine grid fine mesh array */
                memset(rfield.fine_opac_zone , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) );
                /* also zero out fine opacity array */
                memset(rfield.fine_opt_depth , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) );
        }
        return;
}