/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * 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.dat */
#include "cddefines.h"
#include "rfield.h"
#include "physconst.h"
#include "path.h"
#include "doppvel.h"
#include "dense.h"
#include "trace.h"
#include "opacity.h"
#include "ipoint.h"
#include "continuum.h"

/* read the continuum definition from the file continuum_mesh.dat */
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 */
  int 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 int lgPntEval = FALSE;

#	ifdef DEBUG_FUN
	fputs( "<+>ContCreateMesh()\n", debug_fp );
#	endif

	/* 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];
		}
		
#		ifdef DEBUG_FUN
		fputs( " <->ContCreateMesh()\n", debug_fp );
#		endif
		return;
	}
	else
	{
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, " ContCreateMesh called first time.\n" );
		}
		lgPntEval = TRUE;
	}

	/* read in the continuum mesh resolution definition */
	/* >>chng 01 sep 29, add external file "continuum_mesh.dat" 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
	++rfield.nupper;*/
	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.");
		puts( "[Stop in ContCreateMesh]" );
		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();

	/* ================================================================ */
	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.dat */
	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] = (float)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] = (float)(alf*rfield.anu[i]*rfield.anu[i]*3.858e-25);
		rfield.csigc[i] = (float)(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.trans_coef[i] = 1.f;
	}

#	ifdef DEBUG_FUN
	fputs( " <->ContCreateMesh()\n", debug_fp );
#	endif
	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 */
  int lgCount )
{
	long int i, 
	  nbin;
	float widtot;
	double aaa , bbb;

#	ifdef DEBUG_FUN
	fputs( "<+>fill()\n", debug_fp );
#	endif

	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;
#		ifdef DEBUG_FUN
		fputs( " <->fill()\n", debug_fp );
#		endif
		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] );
		puts( "[Stop in fill]" );
		cdEXIT(EXIT_FAILURE);
	}

	ASSERT( *ipnt < continuum.nStoredBands );

	continuum.ifill0[*ipnt] = *n0 - 1;
	continuum.filbnd[*ipnt] = (float)fenlo;
	continuum.filbnd[*ipnt+1] = (float)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] = (float)(log10(fenhi/fenlo)/nbin);

	if( continuum.fildel[*ipnt] < 0.01 )
	{
		continuum.filres[*ipnt] = (float)(log(10.)*continuum.fildel[*ipnt]);
	}
	else
	{
		continuum.filres[*ipnt] = (float)((pow(10.,2.*continuum.fildel[*ipnt]) - 1.)/2./
			pow(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();
		puts( "[Stop in fill]" );
		cdEXIT(EXIT_FAILURE);
	}

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

		rfield.anu[i+continuum.ifill0[*ipnt]] = (float)(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);

#	ifdef DEBUG_FUN
	fputs( " <->fill()\n", debug_fp );
#	endif
	return;
}

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

#	ifdef DEBUG_FUN
	fputs( "<+>ChckFill()\n", debug_fp );
#	endif

	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 )
	{
		puts( "[Stop in ChckFill]" );
		cdEXIT(EXIT_FAILURE);
	}

}

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

#	ifdef DEBUG_FUN
	fputs( "<+>rfield_opac_malloc()\n", debug_fp );
#	endif

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

	/*lint -e737 loss of sign in promotion from long to unsigned long - LOTS of these */

	/* >>chng 03 feb 12, add fine mesh fine grid fine opacity array to keep track of line overlap */
	/*TODO	3	consider making the fine opacity array a double.  with a float, the opacity 
	 * itself often becomes a denormalized number, it then becomes significant when multiplied
	 * by dr - can cause numerical noise.  this is why the coarse opacity array is a double */

	/* frequency range in rydberg needed for all resonance lines */
	rfield.fine_ener_lo = 0.1f;
	rfield.fine_ener_hi = 1500.f;

	/* width will be exactly 1 km/s - FINE_OPAC_VELOCITY_WIDTH is 1e5 cm/s
	rfield.fine_resol = FINE_OPAC_VELOCITY_WIDTH / SPEEDLIGHT; */

	/* >>chng 05 feb 27, have resolution of fine mesh depend on lowest temperature 
	 * 1 km/s resolution was arbitrary 
	 * first find width (cm/s) of Fe line at lowest temperature to be considered */
	/* rfield.fine_opac_velocity_width is width per cell, cm/s */
	/* do not let temperature fall below 1e4 K for choosing stop temperature,
	 * to prevent fine opacity grid from becoming huge */
	/*>>chng 06 feb 08, 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 */
	rfield.fine_opac_velocity_width = 
		/*(float)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*MAX2(1e4,StopCalc.tend)/dense.AtomicWeight[rfield.fine_opac_nelem]+*/
		(float)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*  1e4  /dense.AtomicWeight[rfield.fine_opac_nelem]+
		POW2(DoppVel.TurbVel)) / 
		/* we want fine_opac_nresolv continuum elements across this line */
		rfield.fine_opac_nresolv;

	/* 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 oopacity array to ghost the main low-resolution array */
	if( (rfield.fine_opac = (float *)MALLOC(sizeof(float)*(unsigned)rfield.nfine_malloc )  )==NULL )
		BadMalloc();
	memset(rfield.fine_opac , 0 , (unsigned long)rfield.nfine_malloc*sizeof(float) );

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

	if( (rfield.fine_anu = (float *)MALLOC(sizeof(float)*(unsigned)rfield.nfine_malloc )  )==NULL )
		BadMalloc();

	/* 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 * (float)pow( (1.+rfield.fine_resol), (i+1.) );
	}
	/* done with fine array */

	/* used to count number of lines per cell */
	if( (rfield.line_count = (long *)MALLOC(sizeof(long)*(unsigned)NCELL ) )==NULL )
		BadMalloc();
	for( i=0; i<rfield.nupper; ++i)
	{
		rfield.line_count[i] = 0;
	}
	if( (rfield.anu = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.AnuOrg = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.widflx = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.anulog = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL ) 
		BadMalloc();
	if( (rfield.anusqr = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL ) 
		BadMalloc();
	if( (rfield.anu2 = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL ) 
		BadMalloc();
	if( (rfield.anu3 = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.flux = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.flux_accum = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.convoc = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.OccNumbBremsCont = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.OccNumbIncidCont = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.OccNumbDiffCont = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.OccNumbContEmitOut = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.ConEmitLocal = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.ConEmitReflec = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.ConEmitOut = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.ConInterOut = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.ConRefIncid = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.SummedCon = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.SummedDif = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.SummedDifSave = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.SummedOcc = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.ConOTS_local_photons = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.TotDiff2Pht = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.ConOTS_local_OTS_rate = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.otslin = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL ) 
		BadMalloc();
	if( (rfield.otscon = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL ) 
		BadMalloc();
	if( (rfield.otslinNew = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.otsconNew = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.outlin = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL ) 
		BadMalloc();
	if( (rfield.outlin_noplot = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL ) 
		BadMalloc();
	if( (rfield.reflin = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.FluxSave = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.trans_coef = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.ipnt_coarse_2_fine = (long int*)MALLOC((size_t)(rfield.nupper*sizeof(long int)) ) ) == NULL )
		BadMalloc();
	
	/* chng 02 may 16, by Ryan...added array for gaunt factors for ALL charges, malloc here.	*/
	/* First index is EFFECTIVE CHARGE MINUS ONE!	*/
	if( (rfield.gff = (float**)MALLOC((size_t)((LIMELM+1)*sizeof(float*)) ) ) == NULL )
		BadMalloc();
	for( i = 1; i <= LIMELM; i++ )
	{
		if( (rfield.gff[i] = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
			BadMalloc();
	}

	if( (rfield.csigh = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL ) 
		BadMalloc();
	if( (rfield.csigc = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.ConTabRead = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();

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

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


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

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

	if( (opac.TauAbsFace = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	memset( opac.TauAbsFace , 0 , rfield.nupper*sizeof(float) );

	if( (opac.TauScatFace = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (opac.e2TauAbs = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (opac.ExpmTau = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (opac.tmn = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();

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

	if( (opac.TauAbsGeo = (float**)MALLOC((size_t)(2*sizeof(float *)) ) ) == NULL )
		BadMalloc();
	if( (opac.TauScatGeo = (float**)MALLOC((size_t)(2*sizeof(float *)) ) ) == NULL )
		BadMalloc();
	if( (opac.TauTotalGeo = (float**)MALLOC((size_t)(2*sizeof(float *)) ) ) == NULL )
		BadMalloc();

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

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

	/* say that space exists */
	lgRfieldMalloced = TRUE;

#	ifdef DEBUG_FUN
	fputs( " <->rfield_opac_malloc()\n", debug_fp );
#	endif
}
/*lint +e737 loss of sign */


/* read the continuum definition from the file continuum_mesh.dat */
static void read_continuum_mesh( void )
{
	FILE *ioDATA;
	char chFilename[FILENAME_PATH_LENGTH_2],
		chLine[INPUT_LINE_LENGTH];
	long i;
	int lgEOL;
	long i1 , i2 , i3;

	/* check on path is file not here and path set */
	/* path was parsed in getset */
	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , "continuum_mesh.dat" );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , "continuum_mesh.dat" );
	}

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

	if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " read_continuum_mesh could not open %s\n",chFilename );
		if( lgDataPathSet == TRUE )
		{
			fprintf( ioQQQ, " even tried path\n" );
			fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
			fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
		}

		puts( "[Stop in read_continuum_mesh]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* first line is a version number and does not count */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
	{
		fprintf( ioQQQ, " read_continuum_mesh could not read first line of continuum_mesh.dat.\n");
		puts( "[Stop in read_continuum_mesh]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* count how many lines are in the file, ignoring all lines
	 * starting with '#' */
	continuum.nStoredBands = 0;
	while( fgets( 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 */
	if( (continuum.filbnd = 
		((float *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(float )))) == NULL )
		BadMalloc();
	if( (continuum.fildel = 
		((float *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(float )))) == NULL )
		BadMalloc();
	if( (continuum.filres = 
		((float *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(float )))) == NULL )
		BadMalloc();
	if( (continuum.ifill0 = 
		((long *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(long )))) == NULL )
		BadMalloc();
	if( (continuum.StoredEnergy = 
		((double *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(double )))) == NULL )
		BadMalloc();
	if( (continuum.StoredResolution = 
		((double *)MALLOC( (size_t)(continuum.nStoredBands+1)*sizeof(double )))) == NULL )
		BadMalloc();

	/* 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.dat.\n");
		puts( "[Stop in read_continuum_mesh]" );
		cdEXIT(EXIT_FAILURE);
	}

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

	i = 1;
	/* level 1 magic number */
	i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
	i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
	i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);

	/* the following is the set of numbers that appear at the start of continuum_mesh.dat 01 08 10 */
	if( ( i1 != 1 ) || ( i2 != 9 ) || ( i3 != 29 ) )
	{
		fprintf( ioQQQ, 
			" read_continuum_mesh: the version of continuum_mesh.dat is not the current version.\n" );
		fprintf( ioQQQ, 
			" I expected to find the number 01 09 29 and got %li %li %li instead.\n" ,
			i1 , i2 , i3 );
		fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
		puts( "[Stop in read_continuum_mesh]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* this starts at 1 not 0 since zero is reserved for the
	 * dummy line */
	continuum.nStoredBands = 0;
	while( fgets( 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,INPUT_LINE_LENGTH,&lgEOL);
			continuum.StoredResolution[continuum.nStoredBands] = FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);

			/* option to enter numbers as logs if less than zero */
			if( continuum.StoredEnergy[continuum.nStoredBands]<0. )
				continuum.StoredEnergy[continuum.nStoredBands] = 
				pow(10.,continuum.StoredEnergy[continuum.nStoredBands]);

			if( continuum.StoredResolution[continuum.nStoredBands]<0. )
				continuum.StoredResolution[continuum.nStoredBands] = 
				pow(10.,continuum.StoredResolution[continuum.nStoredBands]);

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

			if( continuum.StoredResolution[continuum.nStoredBands] == 0. )
			{
				fprintf( ioQQQ, 
					" read_continuum_mesh: A continuum resolution was zero - this is not allowed.\n" );
				puts( "[Stop in read_continuum_mesh]" );
				cdEXIT(EXIT_FAILURE);
			}

			++continuum.nStoredBands;
		}
	}

	/* 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" );
			puts( "[Stop in read_continuum_mesh]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
	if( continuum.StoredEnergy[continuum.nStoredBands-1]!=0 )
	{
		fprintf( ioQQQ, 
			" read_continuum_mesh: The last continuum array energies must be zero.\n" );
		puts( "[Stop in read_continuum_mesh]" );
		cdEXIT(EXIT_FAILURE);
	}
}
