/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*ContSetIntensity derive intensity of incident continuum */
/*extin do extinction of incident continuum as set by extinguish command */
/*sumcon sums L and Q for net incident continuum */
/*ptrcer show continuum pointers in real time following drive pointers command */
/*conorm normalize continuum to proper intensity */
/*qintr integrates Q for any continuum between two limits, used for normalization */
/*pintr integrates L for any continuum between two limits, used for normalization */
#include "cddefines.h"
#include "physconst.h"
#include "iso.h"
#include "extinc.h"
#include "noexec.h"
#include "ionbal.h"
#include "hextra.h"
#include "trace.h"
#include "dense.h"
#include "oxy.h"
#include "prt.h"
#include "heavy.h"
#include "rfield.h"
#include "phycon.h"
#include "called.h"
#include "hydrogenic.h"
#include "timesc.h"
#include "neutrn.h"
#include "bit32.h"
#include "secondaries.h"
#include "opacity.h"
#include "thermal.h"
#include "ipoint.h"
#include "atmdat.h"
#include "rt.h"
#include "pressure.h"
#include "radius.h"
#include "geometry.h"
#include "grainvar.h"
#include "continuum.h"

/* these are weights used for continuum integration */
static double aweigh[4]={-0.4305682,-0.1699905, 0.1699905, 0.4305682}; 
static double fweigh[4]={ 0.1739274, 0.3260726, 0.3260726, 0.1739274};

/*conorm normalize continuum to proper intensity */
static void conorm(void);

/*pintr integrates L for any continuum between two limits, used for normalization */
static double pintr(double penlo, 
  double penhi);

/*qintr integrates Q for any continuum between two limits, used for normalization */
static double qintr(double *qenlo, 
  double *qenhi);


/*sumcon sums L and Q for net incident continuum */
static void sumcon(long int il, 
  long int ih, 
  float *q, 
  float *p, 
  float *panu);

/*extin do extinction of incident continuum as set by extinguish command */
static void extin(float *ex1ryd);

/*ptrcer show continuum pointers in real time following drive pointers command */
static void ptrcer(void);

/* return value is zero if all is ok, 1 if something bad happened */
int ContSetIntensity(void)
{
	int lgCheckOK;

	long int i, 
	  ip, 
	  j, 
	  n;

	float EdenHeav, 
	  ex1ryd, 
	  factor, 
	  occ1, 
	  p, 
	  p1, 
	  p2, 
	  p3, 
	  p4, 
	  p5, 
	  p6, 
	  p7, 
	  pgn, 
	  phe, 
	  pheii, 
	  qgn, 
	  temp;

	float xIoniz;

	double rec,
	  wanu[4],
	  alf, 
	  bet, 
	  fntest, 
	  fsum, 
	  ecrit, 
	  tcompr, 
	  tcomp, 
	  r2ov1, 
	  r3ov2;

	double amean, 
	  amean2, 
	  amean3, 
	  peak, 
	  wfun[4];

	long int nelem , ion;

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

	/* set continuum */
	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " ContSetIntensity called.\n" );
	}

	/* find normalization factors for the continua - this decides whether continuum is
	 * per unit area of luminoisity, and which is desired final product */
	conorm();

	/* define factors to convert rfeld.flux array into photon occupation array OCCNUM
	 * by multiplication */
	factor = (float)(EN1RYD/PI4/FR1RYD/HNU3C2);

	/*------------------------------------------------------------- */
	lgCheckOK = TRUE;
	for( i=0; i < rfield.nupper; i++ )
	{
		/* this was original anu array with no continuum averaging */
		rfield.anu[i] = rfield.AnuOrg[i];
		rfield.ContBoltz[i] = 0.;
		fsum = 0.;
		amean = 0.;
		amean2 = 0.;
		amean3 = 0.;

		for( j=0; j < 4; j++ )
		{
			wanu[j] = rfield.anu[i] + rfield.widflx[i]*aweigh[j];
			/* >>chng 02 jul 16, add test on continuum limits -
			 * this was exceeded when resolution set very large */
			wanu[j] = MAX2( wanu[j] , rfield.emm );
			wanu[j] = MIN2( wanu[j] , rfield.egamry );

			wfun[j] = fweigh[j]*ffun(wanu[j]);
			fsum += wfun[j];
			amean += wanu[j]*wfun[j];
			amean2 += wanu[j]*wanu[j]*wfun[j];
			amean3 += wanu[j]*wanu[j]*wanu[j]*wfun[j];
		}

		rfield.flux[i] = (float)(fsum*rfield.widflx[i]);
		if( rfield.flux[i] > 0. )
		{
			rfield.anu[i] = (float)(amean2/amean);
			rfield.anu2[i] = (float)(amean3/amean);
			/* fix conversion factor for occupation number */
		}

		else if( rfield.flux[i] == 0. )
		{
			rfield.anu2[i] = rfield.anu[i]*rfield.anu[i];
		}

		else
		{
			rfield.anu2[i] = rfield.anu[i]*rfield.anu[i];
			fprintf( ioQQQ, " negative continuum returned at%6ld%10.2e%10.2e\n", 
			  i, rfield.anu[i], rfield.flux[i] );
			lgCheckOK = FALSE;
		}
		rfield.anu3[i] = rfield.anu2[i]*rfield.anu[i];

		rfield.ConEmitReflec[i] = 0.;
		rfield.ConEmitOut[i] = 0.;
		rfield.convoc[i] = factor/rfield.widflx[i]/rfield.anu2[i];

		/* following are Compton exchange factors from Tarter */
		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.anu2[i]*3.858e-25);
		rfield.csigc[i] = (float)(alf*bet*rfield.anu[i]*3.858e-25);
	}

#if 0
	/* commented out since we must conserve energy, and continuum was set with old widflx */
	/* 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.widflx[i] = ((rfield.anu[i+1] - rfield.anu[i]) + (rfield.anu[i] - rfield.anu[i-1]))/2.f;
	}
#endif

	if( !lgCheckOK )
	{
		ShowMe();
		puts( "[Stop in ContSetIntensity]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( trace.lgTrace && trace.lgComBug )
	{
		fprintf( ioQQQ, "\n\n Compton heating, cooling coefficients \n" );
		for( i=0; i < rfield.nupper; i += 2 )
		{
			fprintf( ioQQQ, "%6ld%10.2e%10.2e%10.2e", i, rfield.anu[i], 
			  rfield.csigh[i], rfield.csigc[i] );
		}
		fprintf( ioQQQ, "\n" );
	}

	/* option to check frequencies in real time, drive pointers command,
	 * routine is below, is file static */
	if( trace.lgPtrace )
		ptrcer();

	for( i=0; i < rfield.nupper; i++ )
	{
		/* define array of LOG10( nu(ryd)) */
		rfield.anulog[i] = (float)log10(rfield.anu[i]);
	}

	/* extinguish continuum if set on */
	extin(&ex1ryd);

	/* now find peak of hydrogen ionizing continuum - for PDR calculations
	 * this will remain equal to 1 since the loop will not execute */
	prt.ipeak = 1;
	peak = 0.;

	for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ )
	{
		if( rfield.flux[i]*rfield.anu[i]/rfield.widflx[i] > (float)peak )
		{
			/* prt.ipeak points to largest f_nu at H-ionizing energies
			 * and is passed to other parts of code */
			/* i+1 to keep ipeak on fortran version of energy array */
			prt.ipeak = i+1;
			peak = rfield.flux[i]*rfield.anu[i]/rfield.widflx[i];
		}
	}

	/* say what type of cpu this is, if desired */
	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " ContSetIntensity: The peak of the H-ion continuum is at%10.3e\n", 
		  rfield.anu[prt.ipeak-1] );
		if( bit32.lgBit32 )
		{
			fprintf( ioQQQ, " ContSetIntensity: this is a 32-bit cpu.\n" );
		}
		else
		{
			fprintf( ioQQQ, " ContSetIntensity: this is a long-word cpu.\n" );
		}
	}

	/* find highest energy to consider in continuum flux array
	 * peak is the peak product nu*flux */
	peak = rfield.flux[prt.ipeak-1]/rfield.widflx[prt.ipeak-1]*
	  rfield.anu2[prt.ipeak-1];

	if( peak > 1e38 && bit32.lgBit32 )
	{
		fprintf( ioQQQ, " ContSetIntensity: this appears to be a 32-bit cpu.\n" );
		fprintf( ioQQQ, " The continuum is too intense to compute with this cpu.  Use a long-word cpu or a fainter continuum. (This is the nu*f_nu test)\n" );
		fprintf( ioQQQ, " Sorry.\n" );
		puts( "[Stop in ContSetIntensity]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* FluxFaint set in zero.c; normally 1e-10 */
	/* this will be faiintest level of continuum we want to consider.
	 * peak was set above, is peak of hydrogen ionizing radiation field, 
	 * and is zero if no H-ionizing radiation */
	fntest = peak*rfield.FluxFaint;
	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		/* print flux array then quit */
		if( DEBUG_LOC )
		{
			for( i=0; i<rfield.nupper; ++i )
			{
				fprintf(ioQQQ," consetintensityBUGGG\t%.2e\t%.2e\n" , 
					rfield.anu[i] , rfield.flux[i]/rfield.widflx[i] );
			}
			cdEXIT(EXIT_SUCCESS);
		}
	}

	if( fntest > 0. )
	{
		/* this test is not done in pdr conditions where NO H-ionizing radiation,
		 * since fntest is zero*/
		i = rfield.nupper;
		while( i > prt.ipeak && 
			rfield.flux[i-1]*rfield.anu2[i-1]/rfield.widflx[i-1] < (float)fntest )
		{
			--i;
		}
	}
	else
	{
		/* when no H-ionizing radiation set to Lyman edge */
		i = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s];
	}

	/* 
	 * this line of code dates from 1979 and IOA Cambridge.  removed july 19 95
	 * I think it was the last line of the original Cambridge source
	   nflux = MAX( ineon(1)+4 , i )
	 */

	/* >>chng 99 apr 28, reinstate the rfield.FluxFaint limit with nflux */
	rfield.nflux = i;

	/* trim down nflux, was set to rfield.nupper, the dimension of all vectors, in zero.c,
	 * in ContCreatePointers was set to nupper, the number of cells needed to get up to the
	 * high energy limit of the code */
	while( rfield.flux[rfield.nflux-1] < SMALLFLOAT && rfield.nflux > 1 )
	{
		--rfield.nflux;
	}

	if( rfield.nflux == 1 )
	{
		fprintf( ioQQQ, " This incident continuum appears to have no radiation.\n" );
		fprintf( ioQQQ, " Sorry.\n" );
		puts( "[Stop in ContSetIntensity]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* >>chng 04 oct 10, add this limit - arrays will malloc to nupper, but will add unit
	 * continuum to [nflux] - this must be within array */
	rfield.nflux = MIN2( rfield.nupper-1 , rfield.nflux );

	/* check that continuum defined everywhere - look for zero's and comment if present */
	continuum.lgCon0 = FALSE;
	ip = 0;
	for( i=0; i < rfield.nflux; i++ )
	{
		if( rfield.flux[i] == 0. )
		{
			if( called.lgTalk && !continuum.lgCon0 )
			{
				fprintf( ioQQQ, " Setcon: continuum has zero intensity starting at %11.4e Ryd.\n", 
				  rfield.anu[i] );
				continuum.lgCon0 = TRUE;
			}
			++ip;
		}
	}

	if( continuum.lgCon0 && called.lgTalk )
	{
		fprintf( ioQQQ, 
			"%6ld cells in the incident continuum have zero intensity.  Problems???\n", 
		  ip );
	}

	/*begin sanity check */
	lgCheckOK = TRUE;
	for( i=1; i < rfield.nflux; i++ )
	{
		if( rfield.flux[i] < 0. )
		{
			fprintf( ioQQQ, 
				" Continuum has negative intensity at%11.4e Ryd=%10.2e %4.4s %4.4s\n", 
			  rfield.anu[i], rfield.flux[i], rfield.chLineLabel[i]
			  , rfield.chContLabel[i] );
			lgCheckOK = FALSE;
		}
		else if( rfield.anu[i] <= rfield.anu[i-1] )
		{
			fprintf( ioQQQ, 
				" Continuum energies not in increasing order: energies follow\n" );
			fprintf( ioQQQ, 
				"%3ld %10.2e%3ld %10.2e%3ld %10.2e\n", 
			  i -1 , rfield.anu[i-1], i, rfield.anu[i], i +1, rfield.anu[i+1] );
			lgCheckOK = FALSE;
		}
	}

	if( !lgCheckOK )
	{
		ShowMe();
		puts( "[Stop in ContSetIntensity]" );
		cdEXIT(EXIT_FAILURE);
	}
	/*end sanity check */

	/* turn off recoil ionization if high energy < 190R */
	if( rfield.anu[rfield.nflux-1] <= 190 )
	{
		ionbal.lgCompRecoil = FALSE;
	}

	/* sum photons and energy, save mean */

	/* sum from low energy to Balmer edge */
	sumcon(1,iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p]-1,&rfield.qrad,&prt.pradio,&p1);

	/* sum over Balmer continuum */
	sumcon(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2],iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1,&rfield.qbal,&prt.pbal,&p1);

	/* sum from Lyman edge to HeI edge */
	sumcon(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s],
		iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]-1,&prt.q,&p,&p2);

	/* sum from HeI to HeII edges */
	sumcon(iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0],
		iso.ipIsoLevNIonCon[ipH_LIKE][1][ipH1s]-1,&rfield.qhe,&phe,&p3);

	/* sum from Lyman edge to carbon k-shell */
	sumcon(iso.ipIsoLevNIonCon[ipH_LIKE][1][ipH1s],opac.ipCKshell-1,&rfield.qheii,&pheii,&p4);

	/* sum from c k-shell to gamma ray - where pairs start */
	sumcon( opac.ipCKshell , rfield.ipEnerGammaRay-1 , &prt.qx,
		&prt.xpow ,  &p5);

	/* complete sum up to high energy limit */
	sumcon(rfield.ipEnerGammaRay,rfield.nflux,&prt.qgam,&prt.GammaLumin, &p6);

	/* find to estimate photoerosion timescale */
	n = ipoint(7.35e5);
	sumcon(n,rfield.nflux,&qgn,&pgn,&p7);
	timesc.TimeErode = qgn;

	/* find Compton temp */
	tcompr = (p1 + p2 + p3 + p4 + p5 + p6)/(prt.pradio + prt.pbal + 
	  p + phe + pheii + prt.xpow + prt.GammaLumin);

	tcomp = tcompr/(4.*6.272e-6);

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, 
			" mean photon energy=%10.3eR =%10.3eK low, high nu=%12.4e%12.4e\n", 
		  tcompr, tcomp, rfield.anu[0] - rfield.widflx[0]/2., rfield.anu[rfield.nflux-1] + 
		  rfield.widflx[rfield.nflux-1]/2. );
	}

	/* this is total power in ionizing radiation */
	prt.powion = p + phe + pheii + prt.xpow + prt.GammaLumin;

	/* this is the total photon luminosity */
	continuum.TotalLumin = prt.pradio + prt.powion + prt.pbal;

	/* this is placed into the line stack on the first zone, then
	 * reset to zero, to end up with luminosity in the emission lines array.
	 * at end of iteration it is reset to TotalLumin */
	continuum.totlsv = continuum.TotalLumin;

	/* total H-ionizing photon number, */
	rfield.qhtot = prt.q + rfield.qhe + rfield.qheii + prt.qx + prt.qgam;

	/* ftotal photon number, all energies  */
	rfield.qtot = rfield.qhtot + rfield.qbal + rfield.qrad;

	if( prt.powion <= 0. && called.lgTalk )
	{
		rfield.lgHionRad = TRUE;
		fprintf( ioQQQ, " >>>>There is no hydrogen-ionizing radiation.<<<<\n" );
		fprintf( ioQQQ, " >>>>Was this intended?<<<<\n" );
		fprintf( ioQQQ, "  \n" );
		/* >>chng 97 mar 18, add following sanity check - stop if no Paschen
		 * ionizing radiaion since even metals will be totally neutral */
		sumcon(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][3],iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1,&factor,&temp,&p1);
		if( factor <= 0. )
		{
			fprintf( ioQQQ, " >>>>There is no sodium-ionizing radiation.<<<<\n" );
			fprintf( ioQQQ, " >>>>This code is not approriate for these conditions.<<<<\n" );
			
#			ifdef DEBUG_FUN
			fputs( " <->ContSetIntensity()\n", debug_fp );
#			endif
			return(1);
		}
	}

	else
	{
		rfield.lgHionRad = FALSE;
	}

	/* option to add energy deposition due to fast neutrons, as frc of tot lum
	 * efficiency default is unity */
	if( neutrn.lgNeutrnHeatOn )
	{
		neutrn.totneu = (float)(neutrn.effneu*continuum.TotalLumin*pow(10.f,neutrn.frcneu));
	}
	else
	{
		neutrn.totneu = (float)0.;
	}

	/* temp correspond to energy density, printed in STARTR */
	phycon.TEnerDen = (float)pow(continuum.TotalLumin/SPEEDLIGHT/7.56464e-15,0.25);

	/* sanity check for single blackbody, that energy density temperature
	 * is smaller than black body temperature */
	if( rfield.nspec==1 && 
		strcmp( rfield.chSpType[rfield.nspec-1], "BLACK" )==0 )
	{
		/* single black body, now confirm that TEnerDen is less than this temperature,
		 * >>>chng 99 may 02,
		 * in lte these are very very close, factor of 1.00001 allows for numerical
		 * errors, and apparently slightly different atomic coef in different parts
		 * of code.  eventaully all mustuse physonst.h and agree exactly */
		if( phycon.TEnerDen > 1.0001f*rfield.slope[rfield.nspec-1] )
		{
			fprintf( ioQQQ,
				"\n WARNING:  The energy density temperature (%g) is greater than the"
				" black body temperature (%g).  This is unphysical.\n\n",
				phycon.TEnerDen , rfield.slope[rfield.nspec-1]);
		}
	}

	/* incident continuum nu*f_nu at Hbeta and Ly-alpha */
	continuum.cn4861 = (float)(ffun(0.1875)*HPLANCK*FR1RYD*0.1875*0.1875);
	continuum.cn1216 = (float)(ffun(0.75)*HPLANCK*FR1RYD*0.75*0.75);
	continuum.sv4861 = continuum.cn4861;
	continuum.sv1216 = continuum.cn1216;

	/* flux density nu*Fnu = erg / s / cm2
	 * EX1RYD is optional extinction factor at 1 ryd */
	prt.fx1ryd = (float)(ffun(1.000)*HPLANCK*ex1ryd*FR1RYD);

	/* check for plasma frequency - then zero out incident continuum
	 * for energies below this
	 * this is critical electron density, shielding of incident continuum
	 * if electron density is greater than this */
	ecrit = POW2(rfield.anu[0]/2.729e-12);

	if( dense.gas_phase[ipHYDROGEN]*1.2 > ecrit )
	{
		rfield.lgPlasNu = TRUE;
		rfield.plsfrq = (float)(2.729e-12*sqrt(dense.gas_phase[ipHYDROGEN]*1.2));
		rfield.plsfrqmax = rfield.plsfrq;
		rfield.ipPlasma = ipoint(rfield.plsfrq);

		/* save max pointer too */
		rfield.ipPlasmax = rfield.ipPlasma;

		/* now loop over all affected energies, setting incident continuum
		 * to zero there, and counting all as reflected */
		/* >>chng 01 jul 14, from i < ipPlasma to ipPlasma-1 - 
		 * when ipPlasma is 1 plasma freq is not on energy scale */
		for( i=0; i < rfield.ipPlasma-1; i++ )
		{
			/* count as reflected incident continuum */
			rfield.ConRefIncid[i] = rfield.flux[i];
			/* set continuum to zero there */
			rfield.flux[i] = 0.;
		}
	}
	else
	{
		rfield.lgPlasNu = FALSE;
		/* >>chng 01 jul 14, from 0 to 1 - 1 is the first array element on the F scale,
		 * ipoint would return this, so rest of code assumes ipPlasma is 1 plus correct index */
		rfield.ipPlasma = 1;
		rfield.plsfrqmax = 0.;
		rfield.plsfrq = 0.;
	}

	if( rfield.ipPlasma > 1 && called.lgTalk )
	{
		fprintf( ioQQQ, 
			"           !The plasma frequency is %.2e Ryd.  The incident continuum is set to 0 below this.\n", 
		  rfield.plsfrq );
	}

	rfield.occmax = 0.;
	rfield.tbrmax = 0.;
	for( i=0; i < rfield.nupper; i++ )
	{
		/* set up occupation number array */
		rfield.OccNumbIncidCont[i] = rfield.flux[i]*rfield.convoc[i];
		if( rfield.OccNumbIncidCont[i] > rfield.occmax )
		{
			rfield.occmax = rfield.OccNumbIncidCont[i];
			rfield.occmnu = rfield.anu[i];
		}
		/* following product is continuum brightness temperature */
		if( rfield.OccNumbIncidCont[i]*TE1RYD*rfield.anu[i] > rfield.tbrmax )
		{
			rfield.tbrmax = (float)(rfield.OccNumbIncidCont[i]*TE1RYD*rfield.anu[i]);
			rfield.tbrmnu = rfield.anu[i];
		}
		/* save continuum for next iteration */
		rfield.FluxSave[i] = rfield.flux[i];
	}

	/* if continuum brightness temp is large, where does it fall below
	 * 1e4k??? */
	if( rfield.tbrmax > 1e4 )
	{
		i = ipoint(rfield.tbrmnu);
		while( i < rfield.nupper && (rfield.OccNumbIncidCont[i-1]*TE1RYD*
		  rfield.anu[i-1] > 1e4) )
		{
			i += 1;
		}
		rfield.tbr4nu = rfield.anu[i-1];
	}
	else
	{
		rfield.tbr4nu = 0.;
	}

	/* if continuum occ num is large, where does it fall below 1? */
	if( rfield.occmax > 1. )
	{
		i = ipoint(rfield.occmnu)-1;
		while( i < rfield.nupper && (rfield.OccNumbIncidCont[i] > 1.) )
		{
			++i;
		}
		rfield.occ1nu = rfield.anu[i];
	}
	else
	{
		rfield.occ1nu = 0.;
	}

	/* remember if incident radiation field is less than 10*Habing ISM */
	if( continuum.TotalLumin < 1.8e-2 )
	{
		rfield.lgHabing = TRUE;
	}

	/* fix ionization parameter (per hydrogen) at inner edge */
	rfield.uh = (float)(rfield.qhtot/dense.gas_phase[ipHYDROGEN]/SPEEDLIGHT);
	rfield.uheii = (float)((rfield.qheii + prt.qx)/dense.gas_phase[ipHYDROGEN]/SPEEDLIGHT);

	/* guess first temperature and neutral h density */
	if( thermal.ConstTemp > 0. )
	{
		phycon.te = thermal.ConstTemp;
	}
	else
	{
		if( rfield.uh > 0. )
		{
			phycon.te = (float)(20000.+log10(rfield.uh)*5000.);
			phycon.te = (float)MAX2(8000. , phycon.te );
		}
		else
		{
			phycon.te = (float)1000.;
		}
	}
	/* must not call tfidle at this stage since not all vectors have been allocated */
	/*tfidle(TRUE);*/

	/* this is an option to stop after printing header only */
	if( noexec.lgNoExec )
		return(0);

	/* following needed for tfidle */
	dense.eden = dense.gas_phase[ipHYDROGEN];
	dense.xIonDense[ipHYDROGEN][0] = 0.;
	/* this must be zero sinze PresTotCurrent will do radiation pressure due to H */
	iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s] = 0.;

	/* next two just to make sure some values are set */
	/* this sets values of pressure.PresTotlCurr */
	PresTotCurrent();

	/* estimate secondary ionization rate - probably 0, but possible extra
	 * SetCsupra set with "set csupra" */
	/* >>>chng 99 apr 29, added cosmic ray ionization since this is used to get
	 * helium ionization fraction, and was zero in pdr, so He turned off at start,
	 * and never turned back on */
	/* coef on cryden is from highen.c */
	for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
	{
		for( ion=0; ion<nelem+1; ++ion )
		{
			secondaries.csupra[nelem][ion] = 
				secondaries.SetCsupra + hextra.cryden*2e-9f;
		}
	}

	/*********************************************************************
	 *                                                                   *
	 * esimate hydrogen's level of ionization                            *
	 *                                                                   *
	 *********************************************************************/

	/* first estimate level of hydrogen ionization */
	rec = (-9.9765209 + 0.158607055*phycon.telogn[0] + 0.30112749*
	  phycon.telogn[1] - 0.063969007*phycon.telogn[2] + 0.0012691546*
	  phycon.telogn[3])/(1. + 0.035055422*phycon.telogn[0] - 
	  0.037621619*phycon.telogn[1] + 0.0076175048*phycon.telogn[2] - 
	  0.00023432613*phycon.telogn[3]);

	rec = pow(10.,rec)/phycon.te*dense.eden;
	xIoniz = (float)(rfield.qhtot*2e-18 + atmdat_coll_ion(1,1,phycon.te)*dense.eden);
	if( rec < SMALLFLOAT )
		fprintf(ioQQQ," consetintensity insanity, rec is %.3e ionize %.3e\n", 
		rec , xIoniz );
	/* these prints stopped optimized compile in dotnet from blowing the code */
	/*fprintf(ioQQQ," rec isss %.3e\n",rec);*/
	/*fprintf(ioQQQ," xIoniz isss %.3e\n",xIoniz);*/
	r2ov1 = xIoniz/rec;
	/*fprintf(ioQQQ," r2ov1 isss %.3e\n",r2ov1);*/
	dense.xIonDense[ipHYDROGEN][0] = (float)(dense.gas_phase[ipHYDROGEN]/(1. + r2ov1));
	/* >>chng 00 aug 26 add following logic, had used difference between
	 * hden and hi, as in second branch.  this was zero for very low
	 * ionization conditions */
	if( r2ov1 > SMALLFLOAT )
	{
		dense.xIonDense[ipHYDROGEN][1] = (float)(dense.gas_phase[ipHYDROGEN]/( 1. + 1./r2ov1 ) );
	}
	else
	{
		dense.xIonDense[ipHYDROGEN][1] = (float)dense.gas_phase[ipHYDROGEN] - dense.xIonDense[ipHYDROGEN][0];
		/* >>chng 02 nov 25, add this for optimized compile in dotnet and bug mentioned above */
		dense.xIonDense[ipHYDROGEN][1] = MAX2(0.f,dense.xIonDense[ipHYDROGEN][1]);
	}
	ASSERT( dense.xIonDense[ipHYDROGEN][0] >0 && dense.xIonDense[ipHYDROGEN][1]>= 0.);

	if( dense.xIonDense[ipHYDROGEN][1] > 1e-30 )
	{
		iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s] = dense.xIonDense[ipHYDROGEN][0]/dense.xIonDense[ipHYDROGEN][1];
	}
	else
	{
		iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s] = 0.;
	}

	/* now save estimates of whether induced recombination is going
	 * to be important -this is a code pacesetter since GammaBn is slower
	 * than GammaK */
	hydro.lgHInducImp = FALSE;
	for( i=ipH1s; i < iso.numLevels[ipH_LIKE][ipHYDROGEN]; i++ )
	{
		if( rfield.OccNumbIncidCont[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][i]-1] > 0.01 )
			hydro.lgHInducImp = TRUE;
	}

	/*********************************************************************
	 *                                                                   *
	 * esimate helium's level of ionization                            *
	 *                                                                   *
	 *********************************************************************/

	/* only if helium is turned on */
	if( dense.lgElmtOn[ipHELIUM] )
	{
		/* next estimate level of helium singly ionized */
		xIoniz = (float)atmdat_coll_ion(2,2,phycon.te);
		xIoniz = (float)(xIoniz*dense.eden + rfield.qhe*1e-18);
		r2ov1 = xIoniz/rec;

		/* now estimate level of helium doubly ionized */
		xIoniz = (float)atmdat_coll_ion(2,1,phycon.te);
		xIoniz = (float)(xIoniz*dense.eden + rfield.qheii*1e-18);

		/* rough charge dependence */
		rec *= 4.;
		r3ov2 = xIoniz/rec;

		/* now set level of helium ionization */
		if( r2ov1 > 0. )
		{
			dense.xIonDense[ipHELIUM][1] = (float)(dense.gas_phase[ipHELIUM]/(1./r2ov1 + 1. + r3ov2));
			dense.xIonDense[ipHELIUM][0] = (float)(dense.xIonDense[ipHELIUM][1]/r2ov1);
		}
		else
		{
			/* no He ionizing radiation */
			dense.xIonDense[ipHELIUM][1] = 0.;
			dense.xIonDense[ipHELIUM][0] = dense.gas_phase[ipHELIUM];
		}

		dense.xIonDense[ipHELIUM][2] = (float)(dense.xIonDense[ipHELIUM][1]*r3ov2);

		if( dense.xIonDense[ipHELIUM][2] > 1e-30 )
		{
			iso.Pop2Ion[ipH_LIKE][1][ipH1s] = dense.xIonDense[ipHELIUM][1]/dense.xIonDense[ipHELIUM][2];
		}
		else
		{
			iso.Pop2Ion[ipH_LIKE][1][ipH1s] = 0.;
		}
	}
	else
	{
		/* case where helium is turned off */
		dense.xIonDense[ipHELIUM][1] = 0.;
		dense.xIonDense[ipHELIUM][0] = 0.;
		dense.xIonDense[ipHELIUM][2] = 0.;
		iso.Pop2Ion[ipH_LIKE][1][ipH1s] = 0.;
	}

	/* fix number of stages of ionization */
	for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
	{
		if( dense.lgElmtOn[nelem] )
		{
			dense.IonLow[nelem] = 0;
			/* 
			 * IonHigh[n] is the highest stage of ionization present
			 * the IonHigh array index is on the C scake, so [0] is hydrogen
			 * the value is also on the C scale, so element [nelem] can range
			 * from 0 to nelem+1 
			 */
			dense.IonHigh[nelem] = nelem + 1;
			/* >>chng 04 jan 13, add this test, caught by Orly Gnat */
			/* check on actual zero abundances of lower stages - this will only 
			 * happen when ionization is set with element ionization command */
			/* >>chng 04 jun 03, move this test here from conv ionopac do */
			if( dense.lgSetIoniz[nelem] )
			{
				while( dense.SetIoniz[nelem][dense.IonLow[nelem]] == 0. )
					++dense.IonLow[nelem];
				while( dense.SetIoniz[nelem][dense.IonHigh[nelem]] == 0. )
					--dense.IonHigh[nelem];
			}
		}
		else
		{
			/* this element is turned off, make stages impossible */
			dense.IonLow[nelem] = -1;
			dense.IonHigh[nelem] = -1;
		}
	}

	/* estimate electrons from heavies, assuming each at least
	 * 1 times ionized */
	EdenHeav = 0.;
	for( i=2; i < LIMELM; i++ )
	{
		if( dense.lgElmtOn[i] )
		{
			EdenHeav += dense.gas_phase[i];
		}
	}

	/* estimate of electron density */
	dense.eden = 
		dense.xIonDense[ipHYDROGEN][1] + dense.xIonDense[ipHELIUM][1] + 
	  2.*dense.xIonDense[ipHELIUM][2] + EdenHeav + dense.EdenExtra;

	if( dense.EdenSet > 0. )
	{
		dense.eden = dense.EdenSet;
	}

	dense.EdenHCorr = dense.eden;

	if( dense.eden < 0. )
	{
		fprintf( ioQQQ, " Negative electron density results in ContSetIntensity.\n" );
		fprintf( ioQQQ, "%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e\n", 
		  dense.eden, dense.xIonDense[ipHYDROGEN][1], dense.xIonDense[ipHELIUM][1], 
		  dense.xIonDense[ipHELIUM][2], dense.gas_phase[ipCARBON], dense.EdenExtra );
		ShowMe();
		puts( "[Stop in ContSetIntensity]" );
		cdEXIT(EXIT_FAILURE);
	}

	dense.EdenTrue = dense.eden;

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, 
			" ContSetIntensity sets initial EDEN to %.4e, contributors H+=%.2e He+, ++= %.2e %.2e Heav %.2e extra %.2e\n", 
		  dense.eden ,
		  dense.xIonDense[ipHYDROGEN][1],
		  dense.xIonDense[ipHELIUM][1],
          2.*dense.xIonDense[ipHELIUM][2],
		  EdenHeav,
		  dense.EdenExtra);
	}

	occ1 = (float)(prt.fx1ryd/HNU3C2/PI4/FR1RYD);

	/* what is occupation number at 1 Ryd? */
	if( occ1 > 1. )
	{
		rfield.lgOcc1Hi = TRUE;
	}
	else
	{
		rfield.lgOcc1Hi = FALSE;
	}

	if( trace.lgTrace && trace.lgConBug )
	{
		/*  print some useful pointers to ionization edges */
		fprintf( ioQQQ, " H2,1=%5ld%5ld NX=%5ld IRC=%5ld\n", 
		  iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2], 
		  iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s],
		  opac.ipCKshell, 
		  ionbal.ipCompRecoil[ipHYDROGEN][0] );

		fprintf( ioQQQ, " CARBON" );
		for( i=0; i < 6; i++ )
		{
			fprintf( ioQQQ, "%5ld", Heavy.ipHeavy[ipCARBON][i] );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, " OXY" );
		for( i=0; i < 8; i++ )
		{
			fprintf( ioQQQ, "%5ld", Heavy.ipHeavy[ipOXYGEN][i] );
		}
		fprintf( ioQQQ, "%5ld%5ld%5ld\n", opac.ipo3exc[0], 
		  oxy.i2d, oxy.i2p );

		fprintf( ioQQQ, 
			"\n\n                                                  PHOTONS PER CELL (NOT RYD)\n" );
		fprintf( ioQQQ, 
			"\n\n                                        nu, flux, wid, occ \n" );
		fprintf( ioQQQ, 
			" " );

		for( i=0; i < rfield.nflux; i++ )
		{
			fprintf( ioQQQ, "%4ld%10.2e%10.2e%10.2e%10.2e", i, 
			  rfield.anu[i], rfield.flux[i], rfield.widflx[i], 
			  rfield.OccNumbIncidCont[i] );
		}
		fprintf( ioQQQ, " \n" );
	}

	/* zero out some continua related to the ots rates,
	 * prototype and routine in RT_OTS_Update.  This is done here since summed cont will
	 * be set to rfield */
	RT_OTS_Zero();

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " ContSetIntensity returns, nflux=%5ld anu(nflux)=%11.4e eden=%10.2e\n", 
		  rfield.nflux, rfield.anu[rfield.nflux-1], dense.eden );
	}

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

	return(0);
}

/*sumcon sums L and Q for net incident continuum */
static void sumcon(long int il, 
  long int ih, 
  float *q, 
  float *p, 
  float *panu)
{
	long int i, 
	  iupper; /* used as upper limit to the sum */

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

	*q = 0.;
	*p = 0.;
	*panu = 0.;

	/* soft continua may not go as high as the requested bin */
	iupper = MIN2(rfield.nflux,ih);

	/* n.b. - in F77 loop IS NOT executed when IUPPER < IL */
	for( i=il-1; i < iupper; i++ )
	{
		/* sum photon number */
		*q += rfield.flux[i];
		/* en1ryd is needed to stop overflow */
		/* sum flux */
		*p += (float)(rfield.flux[i]*(rfield.anu[i]*EN1RYD));
		/* this sum needed for means */
		*panu += (float)(rfield.flux[i]*(rfield.anu2[i]*EN1RYD));
	}

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

	return;
}

/*ptrcer show continuum pointers in real time following drive pointers command */
static void ptrcer(void)
{
	char chCard[INPUT_LINE_LENGTH];
	/* in case of checking everything, will write errors to this file */
	FILE * ioERRORS=NULL;
	int lgEOL;
	char chKey;
	long int i, 
	  ipnt, 
	  j;
	double pnt, 
	  t1, 
	  t2;

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

	fprintf( ioQQQ, " There are two ways to do this:\n");
	fprintf( ioQQQ, " do you want me to test all the pointers (enter y)\n");
	fprintf( ioQQQ, " or do you want to enter energies yourself? (enter n)\n" );

	if( fgets( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
	{
		fprintf( ioQQQ, " error getting input \n" );
		puts( "[Stop in ptrcer]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* this must be either y or n */
	chKey = chCard[0];

	if( chKey == 'n' )
	{
		/* this branck, enter energies by hand, and see what happens */
		fprintf( ioQQQ, " Enter energy (Ryd); 0 to stop; negative is log.\n" );
		pnt = 1.;
		while( pnt!=0. )
		{
			if( fgets( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
			{
				fprintf( ioQQQ, " error getting input2 \n" );
				puts( "[Stop in ptrcer]" );
				cdEXIT(EXIT_FAILURE);
			}
			/* now get the number off the line */
			i = 1;
			pnt = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

			/* bail if no number at all, or it is zero*/
			if( lgEOL || pnt==0. )
			{
				break;
			}

			/* if number negative then interpret as log */
			if( pnt < 0. )
			{
				pnt = pow(10.,pnt);
			}

			/* get pointer to call */
			ipnt = ipoint(pnt);
			fprintf( ioQQQ, " Cell num%4ld center:%10.2e width:%10.2e low:%10.2e hi:%10.2e convoc:%10.2e\n", 
			  ipnt, rfield.anu[ipnt-1], rfield.widflx[ipnt-1], 
			  rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]/2., 
			  rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]/2., 
			  rfield.convoc[ipnt-1] );
		}
	}

	else if( chKey == 'y' )
	{
		/* first check that ipoint will not crash due to out of range call*/
		if( rfield.anu[0] - rfield.widflx[0]/2.*0.9 < continuum.filbnd[0] )
		{
			fprintf( ioQQQ," ipoint would crash since lowest desired energy of %e ryd is below limit of %e\n",
				rfield.anu[0] - rfield.widflx[0]/2.*0.9 , continuum.filbnd[0] );
			fprintf( ioQQQ," width of cell is %e\n",rfield.widflx[0]);
			puts( "[Stop in ptrcer]" );
			cdEXIT(EXIT_FAILURE);
		}

		else if( rfield.anu[rfield.nflux-1] + rfield.widflx[rfield.nflux-1]/2.*0.9 > 
			continuum.filbnd[continuum.nrange] )
		{
			fprintf( ioQQQ," ipoint would crash since highest desired energy of %e ryd is above limit of %e\n",
				rfield.anu[rfield.nflux-1] + rfield.widflx[rfield.nflux-1]/2.*0.9 , 
				continuum.filbnd[continuum.nrange-1] );
			fprintf( ioQQQ," width of cell is %e\n",rfield.widflx[rfield.nflux]);
			fprintf( ioQQQ," this, previous cells are %e %e\n",
				rfield.anu[rfield.nflux-1],rfield.anu[rfield.nflux-2]);
			puts( "[Stop in ptrcer]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* this branch check everything, write errors to error file */
		fprintf( ioQQQ, " errors output on errors.txt\n");
		fprintf( ioQQQ, " IP(cor),IP(fount),nu lower, upper of found, desired cell.\n" );

		/* error file not open, set to null so we can check later */
		ioERRORS = NULL;
		for( i=0; i < rfield.nflux-1; i++ )
		{
			t1 = rfield.anu[i] - rfield.widflx[i]/2.*0.9;
			t2 = rfield.anu[i] + rfield.widflx[i]/2.*0.9;

			j = ipoint(t1);
			if( j != i+1 )
			{
				/* open file for errors if not already open */
				if( ioERRORS == NULL )
				{
					ioERRORS = fopen("errors.txt" , "w" );
					if( ioERRORS==NULL ) 
					{
						fprintf( ioQQQ," could not create1 errors.txt file\n");
						puts( "[Stop in ptrcer]" );
						cdEXIT(EXIT_FAILURE);
					}
					else
					{
						fprintf( ioQQQ," created errors.txt file with error summary\n");
					}
				}

				fprintf( ioQQQ, " Pointers do not agree for lower bound of cell%4ld, %e\n",
					i, rfield.anu[i]);
				fprintf( ioERRORS, " Pointers do not agree for lower bound of cell%4ld, %e\n",
					i, rfield.anu[i] );
			}

			j = ipoint(t2);
			if( j != i+1 )
			{
				/* open file for errors if not already open */
				if( ioERRORS == NULL )
				{
					ioERRORS = fopen("errors.txt" , "w" );
					if( ioERRORS==NULL ) 
					{
						fprintf( ioQQQ," could not create2 errors.txt file\n");
						puts( "[Stop in ptrcer]" );
						cdEXIT(EXIT_FAILURE);
					}
					else
					{
						fprintf( ioQQQ," created errors.txt file with error summary\n");
					}
				}
				fprintf( ioQQQ, " Pointers do not agree for upper bound of cell%4ld, %e\n", 
					i , rfield.anu[i]);
				fprintf( ioERRORS, " Pointers do not agree for upper bound of cell%4ld, %e\n", 
					i , rfield.anu[i]);
			}

		}
	}

	else
	{
		fprintf( ioQQQ, "I do not understand this key, sorry.  %c\n", chKey );
		puts( "[Stop in ptrcer]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( ioERRORS!=NULL )
		fclose( ioERRORS );
	puts( "[Stop in ptrcer]" );
	cdEXIT(EXIT_FAILURE);
}

/*extin do extinction of incident continuum as set by extinguish command */
static void extin(float *ex1ryd)
{
	long int i, 
	  low;
	double absorb, 
	  factor;

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


	/* modify input continuum by leaky absorber
	 * power law fit to 
	 * >>refer	XUV	extinction	Cruddace et al. 1974 ApJ 187, 497. */
	if( extinc.excolm == 0. )
	{
		*ex1ryd = 1.;
	}
	else
	{
		absorb = 1. - extinc.exleak;
		/*factor = extinc.excolm*6.22e-18;*/
		/* >>chng 01 dec 19, use variable for the constant that gives extinction */
		factor = extinc.excolm*extinc.cnst_col2optdepth;
		/* extinction at 1 and 4 Ryd */
		*ex1ryd = (float)(extinc.exleak + absorb*sexp(factor));

		low = ipoint(extinc.exlow);
		for( i=low-1; i < rfield.nflux; i++ )
		{
			rfield.flux[i] *= (float)(extinc.exleak + absorb*sexp(factor*
			  /* >>chng 01 dec 19, use var with this constant */
			  /*(pow(rfield.anu[i],-2.43))));*/
			  (pow(rfield.anu[i],extinc.cnst_power))));
		}
	}

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

	return;
}

/*conorm normalize continuum to proper intensity */
static void conorm(void)
{
	long int i , nd;
	double xLog_radius_inner, 
	  diff, 
	  f, 
	  flx1, 
	  flx2, 
	  pentrd, 
	  qentrd;

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

	xLog_radius_inner = log10(radius.Radius);

	/* this is a sanity check, it can't happen */
	for( i=0; i < rfield.nspec; i++ )
	{
		if( strcmp(rfield.chRSpec[i],"UNKN") == 0 )
		{
			fprintf( ioQQQ, " UNKN spectral normalization cannot happen.\n" );
			fprintf( ioQQQ, " conorm punts.\n" );
			puts( "[Stop in conorm]" );
			cdEXIT(EXIT_FAILURE);
		}

		else if( strcmp(rfield.chRSpec[i],"SQCM") != 0 && 
			      strcmp(rfield.chRSpec[i],"4 PI") != 0 )
		{
			fprintf( ioQQQ, " chRSpec must be SQCM or 4 PI, and it was %4.4s.  This cannot happen.\n", 
			  rfield.chRSpec[i] );
			fprintf( ioQQQ, " conorm punts.\n" );
			puts( "[Stop in conorm]" );
			cdEXIT(EXIT_FAILURE);
		}


		/* this sanity check makes sure that atlas.mod or werner.mod grids
		 * are for the current version of the code */
		if( strcmp(rfield.chSpType[i],"VOLK ") == 0 )
		{
			/* check that wavelength scale is actually defined outside here */
			ASSERT( rfield.AnuOrg[rfield.nupper-1]>0. );

			diff = fabs(rfield.tNuRyd[rfield.nupper-1][i]-rfield.AnuOrg[rfield.nupper-1])/
			  rfield.AnuOrg[rfield.nupper-1];

			/* this was read from a binary file, so match should be precise */
			if( diff > 10.*FLT_EPSILON )
			{
				if( continuum.ResolutionScaleFactor == 1. )
				{
					fprintf( ioQQQ, "%10.2e%10.2e\n", rfield.AnuOrg[rfield.nupper-1], 
					rfield.tNuRyd[rfield.nupper-1][i] );

					fprintf( ioQQQ,"conorm: The energy grid of the stellar atmosphere file does not agree with the grid in this version of the code.\n" );
					fprintf( ioQQQ,"A stellar atmosphere grid from an old version of the code is probably in place.\n" );
					fprintf( ioQQQ,"A grid for the current version of Cloudy must be generated and used.\n" );
					fprintf( ioQQQ,"This is done with the COMPILE STARS command.\n" );
					fprintf( ioQQQ,"Sorry.\n" );

					puts( "[Stop in conorm]" );
					cdEXIT(EXIT_FAILURE);
				}
				else
				{
					fprintf( ioQQQ,"\n\nThe continuum resolution has been chnaged with the SET CONTINUUM RESOLUTION command.\n" );
					fprintf( ioQQQ,"The compiled stellar continua are not consistent with this changed continuum.\n" );
					fprintf( ioQQQ,"Sorry.\n" );
					puts( "[Stop in conorm]" );
					cdEXIT(EXIT_FAILURE);
				}
			}
		}
	}

	/* this sanity check is that the grains we have read in from opacity files agree
	 * with the energy grid in this version of cloudy */
	for( nd=0; nd < gv.nBin; nd++ )
	{
		int lgErrorFound = FALSE;

		/* these better agree */
		if( gv.bin[nd]->NFPCheck != rfield.nupper )
		{
			fprintf( ioQQQ, 
				" conorm:  expected:%ld found: %ld: number of frequency points in grains data do not match\n",
				 rfield.nupper, gv.bin[nd]->NFPCheck );
			lgErrorFound = TRUE;
		}

		/* check that wavelength scale is actually defined outside here */
		ASSERT( rfield.AnuOrg[rfield.nupper-1]>0. );

		diff = fabs( gv.bin[nd]->EnergyCheck - rfield.AnuOrg[rfield.nupper-1] )/
		  rfield.AnuOrg[rfield.nupper-1];

		/* this was read from an ascii file, so we have to be more lenient
		 * the last constant is determined by the number of decimal places in the .opc file */
		if( diff > MAX2(10.*FLT_EPSILON,3.e-6f) )
		{
			fprintf( ioQQQ, "%14.6e%14.6e: frequencies of last grid point do not match\n",
				 rfield.AnuOrg[rfield.nupper-1], gv.bin[nd]->EnergyCheck );
			lgErrorFound = TRUE;
		}

		if( lgErrorFound )
		{
			if( continuum.ResolutionScaleFactor == 1. )
			{
				fprintf( ioQQQ,"conorm: The energy grid of the grain opacity file does not agree with the grid in this version of the code.\n" );
				fprintf( ioQQQ,"A compiled grain grid from an old version of the code is probably in place.\n" );
				fprintf( ioQQQ,"A grid for the current version of Cloudy must be generated and used.\n" );
				fprintf( ioQQQ,"This is done with the COMPILE ALL GRAINS command.\n" );
				fprintf( ioQQQ,"Sorry.\n" );

				puts( "[Stop in conorm]" );
				cdEXIT(EXIT_FAILURE);
			}
			else
			{
				fprintf( ioQQQ,"\n\nThe continuum resolution has been chnaged with the SET CONTINUUM RESOLUTION command.\n" );
				fprintf( ioQQQ,"The compiled grain opacities are not consistent with this changed continuum, so cannot be used.\n" );
				fprintf( ioQQQ,"Sorry.\n" );

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

	/* default is is to predict line intensities, 
	 * but if any continuum specified as luminosity then would override this -
	 * following two values are correct for intensities */
	radius.pirsq = 0.;
	radius.lgPredLumin = FALSE;

	/* check whether ANY luminosities are present */
	for( i=0; i < rfield.nspec; i++ )
	{
		if( strcmp(rfield.chRSpec[i],"4 PI") == 0 )
		{
			radius.pirsq = (float)(1.0992099 + 2.*xLog_radius_inner);
			radius.lgPredLumin = TRUE;
			/* convert down to intensity */
			rfield.totpow[i] -= radius.pirsq;

			if( trace.lgTrace )
			{
				fprintf( ioQQQ, 
					" conorm converts continuum %ld from luminosity to intensity.\n", 
					i );
			}
		}
	}

	/* if total luminosities are present, must have specified a starting radius */
	if( radius.lgPredLumin && !radius.lgRadiusKnown )
	{
		fprintf(ioQQQ,"conorm: - A continuum source was specified as a luminosity, but the inner radius of the cloud was not set.\n");
		fprintf(ioQQQ,"Please set an inner radius.\nSorry.\n");
		puts( "[Stop in conorm]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* convert ionization parameters to number of photons, called "q(h)"
	 * at this stage q(h) and "PHI " are the same */
	for( i=0; i < rfield.nspec; i++ )
	{
		if( strcmp(rfield.chSpNorm[i],"IONI") == 0 )
		{
			/* the log of the ionization parameter was stored here, this converts
			 * it to the log of the number of photons per sq cm */
			rfield.totpow[i] += log10(dense.gas_phase[ipHYDROGEN]) + log10(SPEEDLIGHT);
			strcpy( rfield.chSpNorm[i], "Q(H)" );
			if( trace.lgTrace )
			{
				fprintf( ioQQQ, 
					" conorm converts continuum %ld from ionizat par to q(h).\n", 
					i );
			}
		}
	}

	/* convert x-ray ionization parameter xi to intensity */
	for( i=0; i < rfield.nspec; i++ )
	{
		if( strcmp(rfield.chSpNorm[i],"IONX") == 0 )
		{
			/* this converts it to an intensity */
			rfield.totpow[i] += log10(dense.gas_phase[ipHYDROGEN]) - log10(PI4);
			strcpy( rfield.chSpNorm[i], "LUMI" );
			if( trace.lgTrace )
			{
				fprintf( ioQQQ, " conorm converts continuum%3ld from x-ray ionizat par to I.\n", 
				  i );
			}
		}
	}

	/* indicate whether we ended up with luminosity or intensity */
	if( trace.lgTrace )
	{
		if( radius.lgPredLumin )
		{
			fprintf( ioQQQ, " Cloudy will predict lumin into 4pi\n" );
		}
		else
		{
			fprintf( ioQQQ, " Cloudy will do surface flux for lumin\n" );
		}
	}

	/* if intensity per unit area is predicted then geometric
	 * covering factor must be unity
	 * variable can also be set elsewhere */
	if( !radius.lgPredLumin )
	{
		geometry.covgeo = 1.;
	}

	/* main loop over all continuum shapes to find continuum normalization 
	 * for each one */
	for( i=0; i < rfield.nspec; i++ )
	{
		rfield.ipspec = i;

		/* check that, if laser, bounds include laser energy */
		if( strcmp(rfield.chSpType[rfield.ipspec],"LASER") == 0 )
		{
			if( !( rfield.range[0][rfield.ipspec] < rfield.slope[rfield.ipspec] &&
				rfield.range[1][rfield.ipspec] > rfield.slope[rfield.ipspec]) )
			{
				fprintf(ioQQQ,"DISASTER, a continuum source is a laser at %f Ryd, but the intensity was specified over a range from %f to %f Ryd.\n",
					rfield.slope[rfield.ipspec],
					rfield.range[0][rfield.ipspec],
					rfield.range[1][rfield.ipspec]);
				fprintf(ioQQQ,"Please specify the continuum flux where the laser is active.\n");
				puts( "[Stop in conorm]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		if( trace.lgTrace )
		{
			fprintf( ioQQQ, " conorm continuum number%3ld is shape %s range is %.2e %.2e\n", 
			  i, 
			  rfield.chSpType[i], 
			  rfield.range[0][i], 
			  rfield.range[1][i] );
		}

		if( strcmp(rfield.chSpNorm[i],"RATI") == 0 )
		{
			/* option to scale relative to previous continua
			 * this must come first since otherwise may be too late
			 * BUT ratio cannot be the first continuum source */
			if( trace.lgTrace )
			{
				fprintf( ioQQQ, " conorm this is ratio to 1st con\n" );
			}

			/* check that this is not the first continuum source, we must ratio */
			if( i == 0 )
			{
				fprintf( ioQQQ, " I cant form a ratio if continuum is first source.\n" );
				puts( "[Stop in conorm]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* first find photon flux and Q of prevous continuum */
			rfield.ipspec -= 1;
			flx1 = ffun1(rfield.range[0][i])*rfield.spfac[rfield.ipspec]*
			  rfield.range[0][i];

			/* check that previous continua were not zero where ratio is formed */
			if( flx1 <= 0. )
			{
				fprintf( ioQQQ, " Previous continua were zero where ratio is desired.\n" );
				puts( "[Stop in conorm]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* return pointer to previous (correct) value, find F, Q */
			rfield.ipspec += 1;

			/* we want a continuum totpow as powerful, flx is now desired flx */
			flx1 *= rfield.totpow[i];

			/*.        find flux of this new continuum at that point */
			flx2 = ffun1(rfield.range[1][i])*rfield.range[1][i];

			/* this is ratio of desired to actual */
			rfield.spfac[i] = flx1/flx2;
			if( trace.lgTrace )
			{
				fprintf( ioQQQ, " conorm ratio will set scale fac to%10.3e at%10.2e Ryd.\n", 
				  rfield.totpow[i], rfield.range[0][i] );
			}
		}

		else if( strcmp(rfield.chSpNorm[i],"FLUX") == 0 )
		{
			/* specify flux density
			 * option to use arbitrary frequency or range */
			f = ffun1(rfield.range[0][i]); 
			f = MAX2(1e-37,f); 
			f = log10(f) + log10(rfield.range[0][i]*EN1RYD/FR1RYD);

			f = rfield.totpow[i] - f;
			/* >>chng 96 dec 31, added following test */
			if( bit32.lgBit32 )
			{
				if( f > 35. )
				{
					fprintf( ioQQQ, " Continuum source %ld is too intense for this cpu - is it normalized correctly?\n", 
					  i );
					puts( "[Stop in conorm]" );
					cdEXIT(EXIT_FAILURE);
				}
			}

			rfield.spfac[i] = pow(10.,f);
			if( trace.lgTrace )
			{
				fprintf( ioQQQ, " conorm will set log fnu to%10.3e at%10.2e Ryd.  Factor=%11.4e\n", 
				  rfield.totpow[i], rfield.range[0][i], rfield.spfac[i] );
			}
		}

		else if( strcmp(rfield.chSpNorm[i],"Q(H)") == 0 || 
			strcmp(rfield.chSpNorm[i],"PHI ") == 0 )
		{
			/* some type of photon density entered */
			if( trace.lgTrace )
			{
				fprintf( ioQQQ, " conorm calling qintr range=%11.3e %11.3e desired val is %11.3e\n", 
				  rfield.range[0][i], 
				  rfield.range[1][i] ,
				  rfield.totpow[i]);
			}

			/* the total number of photons over the specified range in
			 * the arbitrary system of units that the code save the continuum shape */
			qentrd = qintr(&rfield.range[0][i],&rfield.range[1][i]);
			/* this is the log of the scale factor that must multiply the 
			 * continuum shape to get the final set of numbers */
			diff = rfield.totpow[i] - qentrd;

			/* >>chng 03 mar 13, from diff < -25 to <-35,
			 * tripped for very low U models used for H2 simulations */
			/*if( bit32.lgBit32 && (diff < -25. || diff > 35.) )*/
			if( bit32.lgBit32 && (diff < -35. || diff > 35.) )
			{
				fprintf( ioQQQ, " Continuum source specified is too extreme for a 32 bit cpu.\n" );
				fprintf( ioQQQ, 
					" The integral over the continuum shape gave (log) %.3e photons, and the command requested (log) %.3e.\n" ,
					qentrd , rfield.totpow[i]);
				fprintf( ioQQQ, 
					" The difference in the log is %.3e.\n" ,
					diff );
				if( diff>0. )
				{
					fprintf( ioQQQ, " The continuum source is too bright.\n" );
				}
				else
				{
					fprintf( ioQQQ, " The continuum source is too faint.\n" );
				}
				/* explain what happened */
				fprintf( ioQQQ, " The usual cause for this problem is an incorrect continuum intensity/luminosity command.\n" );
				fprintf( ioQQQ, " There were a total of %li continuum shape commands entered - the problem is with number %li.\n",
					rfield.nspec , i+1 );
				fprintf( ioQQQ, " If this command is correct then please show this input stream to Gary Ferland\n");
				puts( "[Stop in conorm]" );
				cdEXIT(EXIT_FAILURE);
			}

			else
			{
				rfield.spfac[i] = pow(10.,diff);
			}

			if( trace.lgTrace )
			{
				fprintf( ioQQQ, " conorm finds Q over range from%11.4e-%11.4e Ryd, integral= %10.4e Factor=%11.4e\n", 
				  rfield.range[0][i], 
				  rfield.range[1][i], 
				  qentrd ,
				  rfield.spfac[i] );
			}
		}

		else if( strcmp(rfield.chSpNorm[i],"LUMI") == 0 )
		{
			/* luminosity entered, special since default is TOTAL lumin */
			/*pintr integrates L for any continuum between two limits, used for normalization,
			 * return units are log of ryd cm-2 s-1, last log conv to ergs */
			pentrd = pintr(rfield.range[0][i],rfield.range[1][i]) + 
			  log10(EN1RYD);
			f = rfield.totpow[i] - pentrd;

			/* >>chng 96 dec 31, added following test */
			if( bit32.lgBit32 )
			{
				if( f > 35. )
				{
					fprintf( ioQQQ, " Continuum source%3ld is too intense for this cpu - is it normalized correctly?\n", 
					  i );
					puts( "[Stop in conorm]" );
					cdEXIT(EXIT_FAILURE);
				}
			}

			rfield.spfac[i] = pow(10.,f);
			if( trace.lgTrace )
			{
				fprintf( ioQQQ, " conorm finds luminosity range is%10.3e to %9.3e Ryd, factor is%11.4e\n", 
				  rfield.range[0][i], rfield.range[1][i], 
				  rfield.spfac[i] );
			}
		}

		else
		{
			fprintf( ioQQQ, "What chSpNorm label is this? =%s=\n", rfield.chSpNorm[i]);
			fprintf( ioQQQ, "Insanity has been detected, I cant go on.\n");
			puts( "[Stop in conorm]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* sec part after .or. added June 93 because sometimes spfac=0
		 * got past first test */
		if( 1./rfield.spfac[i] == 0. || rfield.spfac[i] == 0. )
		{
			fprintf( ioQQQ, "conorm finds infinite continuum scale factor.\n" );
			fprintf( ioQQQ, "The continuum is too intense to compute with this cpu.\n" );
			fprintf( ioQQQ, "Were the intensity and luminosity commands switched?\n" );
			fprintf( ioQQQ, "Sorry, but I cannot go on.\n" );
			puts( "[Stop in conorm]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* this is conversion factor for final units of line intensities or luminosities in printout,
	 * will be intensities (==0) unless luminosity is to be printed, or flux at Earth 
	 * pirsq is the log of 4 pi r_in^2 */
	radius.Conv2PrtInten = radius.pirsq;

	/* >>chng 02 apr 25, add option for slit on aperture command */
	if( geometry.iEmissPower == 1  )
	{
		if( radius.lgPredLumin )
		{
			/* factor should be divided by 2 r_in */
			radius.Conv2PrtInten -= (log10(2.) + xLog_radius_inner);
		}
		else if( !radius.lgPredLumin )
		{
			/* this is an error - slit requested but radius is not known */
			fprintf( ioQQQ, "conorm: Aperture slit specified, but not predicting luminosity.\n" );
			fprintf( ioQQQ, "conorm: Please specify an inner radius to determine L.\nSorry\n" );
			puts( "[Stop in conorm]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
	if( geometry.iEmissPower == 0 && radius.lgPredLumin)
	{
		/* leave Conv2PrtInten at zero if not predicting luminosity */
		radius.Conv2PrtInten = log10(2.);
	}

	/* this is option go give final absolute results as flux observed at Earth */
	if( radius.distance > 0. && radius.lgRadiusKnown && prt.lgPrintFluxEarth )
	{
		/*radius.Conv2PrtInten -= 2.*xLog_radius_inner - 2.*log10(radius.distance) ;*/
		/* div (in log) by 4 pi dist^2 */
		radius.Conv2PrtInten -= log10( 4.*PI*POW2(radius.distance) ) ;
	}

	/* normally lines are into 4pi, this is option to do per sr or arcsec^2 */
	if( prt.lgSurfaceBrightness )
	{
		if( radius.pirsq!= 0. )
		{
			/* make sure we are predicting line intensities, not luminosity */
			fprintf( ioQQQ, " Sorry, but both luminosity and surface brightness have been requested for lines.\n" );
			fprintf( ioQQQ, " the PRINT LINE SURFACE BRIGHTNESS command can only be used when lines are predicted per unit cloud area.\n" );
			puts( "[Stop in conorm]" );
			cdEXIT(EXIT_FAILURE);
		}
		if( prt.lgSurfaceBrightness_SR )
		{
			/* we want final units to be per sr */
			radius.Conv2PrtInten -= log10( PI4 ) ;
		}
		else
		{
			/* we want final units to be per square arcsec,
			 * first 4 pi converts to per sr,
			 * there are 4pi sr in 360 deg, so term in square
			 * goes from sr to square sec of arc */
			radius.Conv2PrtInten -= log10( PI4 *POW2( (360./(PI2)*3600.) ) ) ;
		}
	}

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

/*qintr integrates Q for any continuum between two limits, used for normalization */
static double qintr(double *qenlo, 
  double *qenhi)
{
	long int i, 
	  ipHi, 
	  ipLo, 
	  j;
	double qintr_v, 
	  sum, 
	  wanu;

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

	/* returns LOG of number of photons over energy interval */
	sum = 0.;
	/* >>chng 02 oct 27, do not use qg32, always use same method as
	 * routine that does final set of continuum */

	/* this is copy of logic that occurs three times across code */
	ipLo = ipoint(*qenlo);
	ipHi = ipoint(*qenhi);
	/* this is actual sum of photons within band */
	for( i=ipLo-1; i < (ipHi - 1); i++ )
	{
		/*sum += ffun1(rfield.anu[i])*rfield.widflx[i];*/
		for( j=0; j < 4; j++ )
		{
			wanu = rfield.anu[i] + rfield.widflx[i]*aweigh[j];
			/* >>chng 02 jul 16, add test on continuum limits -
			* this was exceeded when resolution set very large */
			wanu = MAX2( wanu , rfield.emm );
			wanu = MIN2( wanu , rfield.egamry );
			sum += fweigh[j]*ffun1(wanu)*rfield.widflx[i];
		}
	}

	if( sum <= 0. )
	{
		fprintf( ioQQQ, " Photon number sum in QINTR is %.3e\n", 
		  sum );
		fprintf( ioQQQ, " This source has no ionizing radiation, and the number of ionizing photons was specified.\n" );
		fprintf( ioQQQ, " This was continuum source number%3ld\n", 
		  rfield.ipspec );
		fprintf( ioQQQ, " Sorry, but I cannot go on.  ANU and FLUX arrays follow.  Enjoy.\n" );
		for( i=0; i < rfield.nupper; i++ )
		{
			fprintf( ioQQQ, "%.2e\t%.2e\n", 
				rfield.anu[i], 
				rfield.flux[i] );
		}
		puts( "[Stop in qintr]" );
		cdEXIT(EXIT_FAILURE);
	}
	else
	{
		qintr_v = log10(sum);
	}

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

/*pintr integrates L for any continuum between two limits, used for normalization,
 * return units are log of ryd cm-2 s-1 */
static double pintr(double penlo, 
  double penhi)
{
	long int i, 
	  j;
	double fsum, 
	  pintr_v, 
	  sum, 
	  wanu, 
	  wfun;
	long int ip1 , ip2;

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

	/* computes log of luminosity in radiation over some intergal
	 * answer is in Ryd per sec */

	sum = 0.;
	/* >>chng 02 oct 27, do not call qg32, do same type sum as
		* final integration */
	/* laser is special since delta function, this is center of laser */
	/* >>chng 01 jul 01, was +-21 cells, change to call to ipoint */
	ip1 = ipoint( penlo );

	ip2 = ipoint( penhi );

	for( i=ip1-1; i < ip2-1; i++ )
	{
		fsum = 0.;
		for( j=0; j < 4; j++ )
		{
			wanu = rfield.anu[i] + rfield.widflx[i]*aweigh[j];
			wfun = fweigh[j]*ffun1(wanu)*wanu;
			fsum += wfun;
		}
		sum += fsum*rfield.widflx[i];
	}

	if( sum > 0. )
	{
		pintr_v = log10(sum);
	}
	else
	{
		pintr_v = -38.;
	}

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

	return( pintr_v );
}

