/* 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 */
/*GrainMakeDiffuse main routine for generating the grain diffuse emission, called by RT_diffuse */
#include "cddefines.h"
#include "physconst.h"
#include "radius.h"
#include "rfield.h"
#include "phycon.h"
#include "dense.h"
#include "hmi.h"
#include "thermal.h"
#include "trace.h"
#include "interpolate.h"
#include "iterations.h"
#include "grainvar.h"
#include "grains.h"
/*lint -e817 Conceivably negative subscript*/
/*lint -e661 possible access of out-of-bounds pointer */

#ifdef STATIC
#undef STATIC
#endif

#ifdef USE_GPROF
#define STATIC
#else
#define STATIC static
#endif


#define NO_ATOMS(ND) (gv.bin[ND]->AvVol*gv.bin[ND]->dustp[0]/ATOMIC_MASS_UNIT/gv.bin[ND]->atomWeight)

/* NB NB -- the sequence below has been carefully chosen and should NEVER be
 *          altered unless you really know what you are doing !! */
/* >>chng 03 jan 16, introduced QH_THIGH_FAIL and started using splint_safe and spldrv_safe
 *                   throughout the code; this solves the parispn_a031_sil.in bug, PvH */
/* >>chng 03 jan 16, rescheduled QH_STEP_FAIL as non-fatal; it can be triggered at very low temps, PvH */
typedef enum {
	/* the following are OK */
	/* 0        1               2              3    */
	QH_OK, QH_ANALYTIC, QH_ANALYTIC_RELAX, QH_DELTA, 
	/* the following are mild errors we already recovered from */
	/*     4              5              6             7        */
	QH_NEGRATE_FAIL, QH_LOOP_FAIL, QH_ARRAY_FAIL, QH_THIGH_FAIL,
	/* any of these errors will prompt qheat to do another iteration */
	/*  8          9             10             11       */
	QH_RETRY, QH_CONV_FAIL, QH_BOUND_FAIL, QH_DELTA_FAIL,
	/* any error larger than QH_NO_REBIN will cause GetProbDistr_LowLimit to return
	 * before even attempting to rebin the results; we may be able to recover though... */
	/*   12          13            14            15       */
	QH_NO_REBIN, QH_LOW_FAIL, QH_HIGH_FAIL, QH_STEP_FAIL,
	/* any case larger than QH_FATAL is truly pathological; there is no hope of recovery */
	/* 16          17           18             19      */
	QH_FATAL, QH_WIDE_FAIL, QH_NBIN_FAIL, QH_REBIN_FAIL
} QH_Code;

/*================================================================================*/
/* definitions relating to quantum heating routines */

/* this is the minimum number of bins for quantum heating calculation to be valid */
static const long NQMIN = 10L;

/* this is the lowest value for dPdlnT that should be included in the modeling */
static const double PROB_CUTOFF_LO = 1.e-15;
static const double PROB_CUTOFF_HI = 1.e-15;
static const double SAFETY = 1.e+8;

/* during the first NSTARTUP steps, use special algorithm to calculate stepsize */
static const long NSTARTUP = 5L;

/* if the average number of multiple events is above this number
 * don't try full quantum heating treatment. */
static const double MAX_EVENTS = 150.;

/* maximum number of tries for quantum heating routine */
/* >>chng 02 jan 30, changed LOOP_MAX from 10 -> 20, PvH */
static const long LOOP_MAX = 20L;

/* if all else fails, divide temp estimate by this factor */
static const double DEF_FAC = 3.;

/* total probability for all bins should not deviate more than this from 1. */
static const double PROB_TOL = 0.02;

/* after NQTEST steps, make an estimate if prob distr will converge in NQGRID steps */
/* >>chng 02 jan 30, change NQTEST from 1000 -> 500, PvH */
static const long NQTEST = 500L;

/* if the ratio fwhm/Umax is lower than this number
 * don't try full quantum heating treatment. */
static const double FWHM_RATIO = 0.1;
/* if the ratio fwhm/Umax is lower than this number
 * don't even try analytic approximation, use delta function instead */
static const double FWHM_RATIO2 = 0.007;

/* maximum number of steps for integrating quantum heating probability distribution */
static const long MAX_LOOP = 2*NQGRID;

/* this is the tolerance used while integrating the temperature distribution of the grains */
static const double QHEAT_TOL = 5.e-3;

/* maximum number of tries before we declare that probability distribution simply won't fit */
static const long WIDE_FAIL_MAX = 3;

/* multipliers for PROB_TOL used in GetProbDistr_HighLimit */
static const double STRICT = 1.;
static const double RELAX = 15.;

/* when rebinning quantum heating results, make ln(qtemp[i]/qtemp[i-1]) == QT_RATIO */
/* this constant determines the accuracy with which the Wien tail of the grain emission
 * is calculated; if x = h*nu/k*T_gr and d = QT_RATIO-1., then the relative accuracy of
 * the flux in the Wien tail is Acc = fabs(x*d^2/12 - x^2*d^2/24). A typical value of x
 * would be x = 15, so QT_RATIO = 1.03 should converge the spectrum to better than 1% */
static const double QT_RATIO = 1.03;


/*================================================================================*/
/* global variables */

/* these data define the enthalpy function for silicates
 * derived from:
 * >>refer	grain	physics	Guhathakurta & Draine, 1989, ApJ, 345, 230
 * coefficients converted to rydberg energy units, per unit atom
 * assuming a density of 3.3 g/cm^3 and pure MgSiFeO4.
 * this is not right, but the result is correct because number
 * of atoms will be calculated using the same assumption. */

/* this is the specific density of silicate in g/cm^3 */
#define DEN_SIL     3.30

/* these are the mean molecular weights per atom for MgSiFeO4 and SiO2 in amu */
#define MW_SIL      24.6051
/*#define MW_SIO2     20.0283*/

static const double tlim[5]={0.,50.,150.,500.,DBL_MAX};
static const double power[4]={2.00,1.30,0.68,0.00};
static const double cval[4]={
	1.40e3/DEN_SIL*ATOMIC_MASS_UNIT*MW_SIL/EN1RYD,
	2.20e4/DEN_SIL*ATOMIC_MASS_UNIT*MW_SIL/EN1RYD,
	4.80e5/DEN_SIL*ATOMIC_MASS_UNIT*MW_SIL/EN1RYD,
	3.41e7/DEN_SIL*ATOMIC_MASS_UNIT*MW_SIL/EN1RYD};


/* initialize phiTilde */
STATIC void qheat_init(long,/*@out@*/double[],/*@out@*/double*);
/* worker routine, this implements the algorithm of Guhathakurtha & Draine */
STATIC void GetProbDistr_LowLimit(long,double,double,double,/*@in@*/double[],/*@in@*/double[],
				  /*@out@*/double[],/*@out@*/double[],/*@out@*/double[],
				  /*@out@*/long*,/*@out@*/double*,long*,/*@out@*/QH_Code*);
/* try two consecutive integration steps using stepsize "step/2." (yielding p[k]),
 * and also one double integration step using stepsize "step" (yielding p2k). */
STATIC double TryDoubleStep(double[],double[],double[],double[],double[],double[],
			    double[],double,/*@out@*/double*,long,long,/*@out@*/int*);
/* calculate logarithmic integral from (x1,y1) to (x2,y2) */
STATIC double log_integral(double,double,double,double);
/* scan the probability distribution for valid range */
STATIC void ScanProbDistr(double[],double[],long,double,long,/*@out@*/long*,/*@out@*/long*,
			  /*@out@*/long*,long*,QH_Code*);
/* rebin the quantum heating results to speed up RT_diffuse */
STATIC long RebinQHeatResults(long,long,long,double[],double[],double[],double[],double[],
			      double[],double[],QH_Code*);
/* calculate approximate probability distribution in high intensity limit */
STATIC void GetProbDistr_HighLimit(long,double,double,double,/*@out@*/double[],/*@out@*/double[],
				   /*@out@*/double[],/*@out@*/double*,/*@out@*/long*,
				   /*@out@*/double*,/*@out@*/QH_Code*);
/* derivative of the enthalpy function dU/dT */
STATIC double uderiv(double,long);
/* enthalpy function */
STATIC double ufunct(double,long,/*@out@*/int*);
/* inverse enthalpy function */
STATIC double inv_ufunct(double,long,/*@out@*/int*);
/* helper function for calculating enthalpy, uses Debye theory */
STATIC double DebyeDeriv(double,long);

/* >>chng 01 oct 29, introduced gv.bin[nd]->cnv_H_pGR, cnv_GR_pH, etc. PvH */


/* main routine for generating the grain diffuse emission, called by RT_diffuse */
void GrainMakeDiffuse()
{
	long i,
	  j,
	  nd;
	double bfac,
	  f,
	  factor,
	  flux,
	  x;

	/* >>chng 01 sep 11, replace array allocation on stack with
	 * MALLOC to avoid bug in gcc 3.0 on Sparc platforms, PvH */
	double *qtemp/*[NQGRID]*/, *qprob/*[NQGRID]*/;

#	ifndef NDEBUG
	int lgNoTdustFailures;
	double BolFlux,
	  Comparison1,
	  Comparison2;
#	endif

	/* this assures 6-digit precision in the evaluation of the exponential below */
	const double LIM1 = pow(2.e-6,1./2.);
	const double LIM2 = pow(6.e-6,1./3.);

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

	factor = 2.*PI4*POW2(FR1RYD/SPEEDLIGHT)*FR1RYD;
	/* >>chng 96 apr 26 upper limt chng from 15 to 75 Peter van Hoof  */
	/* >>chng 00 apr 10 use constants appropriate for double precision, by PvH */
	x = log(0.999*DBL_MAX);

	/* add grain dust emission in IR */
	for( i=0; i < rfield.nflux; i++ )
	{
		/* this is added to one of the local outward beams below, for only one
		 * of the diffuse transfer methods */
		gv.GrainEmission[i] = 0.;
	}

	if( ( qtemp = (double*)MALLOC((size_t)(NQGRID*sizeof(double))) ) == NULL )
		BadMalloc();
	if( ( qprob = (double*)MALLOC((size_t)(NQGRID*sizeof(double))) ) == NULL )
		BadMalloc();
	
	for( nd=0; nd < gv.nBin; nd++ )
	{
		strg_type scase;
		int lgLocalQHeat;
		long qnbin=-200;
		float *grn, threshold;
		double xx;

		/* save emission from this species */
		scase = gv.which_strg[gv.bin[nd]->matType];
		switch( scase )
		{
		case STRG_SIL:
			grn = gv.SilicateEmission;
			break;
		case STRG_CAR:
			grn = gv.GraphiteEmission;
			break;
		default:
			fprintf( ioQQQ, " GrainMakeDiffuse called with unknown storage class: %d\n", scase );
			puts( "[Stop in GrainMakeDiffuse]" );
			cdEXIT(1);
		}

		/* this local copy is necessary to keep lint happy */
		lgLocalQHeat = gv.bin[nd]->lgQHeat;
		/* >>chng 04 nov 09, do not evaluate quantum heating if abundance is negligible, PvH
		 * this prevents PAH's deep inside molecular regions from failing if GrnVryDepth is used */
		/* >>chng 04 dec 31, introduced separate thresholds near I-front and in molecular region, PvH */
		threshold = ( dense.xIonDense[ipHYDROGEN][0]+dense.xIonDense[ipHYDROGEN][1] > hmi.H2_total ) ?
			gv.dstAbundThresholdNear : gv.dstAbundThresholdFar;

		if( lgLocalQHeat && gv.bin[nd]->dstAbund >= threshold )
		{
			qheat(qtemp,qprob,&qnbin,nd);

			if( gv.bin[nd]->lgUseQHeat )
			{
				ASSERT( qnbin > 0 );
			}
		}
		else
		{
			/* >> chng 04 dec 31, repaired bug lgUseQHeat not set when abundance below threshold, PvH */
			gv.bin[nd]->lgUseQHeat = FALSE;
		}

		flux = 1.;
		/* flux can only become zero in the Wien tail */
		/* >>chng 04 jan 25, replaced flux > 0. with (float)flux > 0.f, PvH */
		for( i=0; i < rfield.nflux && (float)flux > 0.f; i++ )
		{
			flux = 0.;
			if( lgLocalQHeat && gv.bin[nd]->lgUseQHeat )
			{
				xx = 1.;
				/* we start at high temperature end and work our way down
				 * until contribution becomes negligible */
				for( j=qnbin-1; j >= 0 && xx > flux*DBL_EPSILON; j-- )
				{
					f = TE1RYD/qtemp[j]*rfield.anu[i];
					if( f < x )
					{
						/* want the number exp(hnu/kT) - 1, two expansions */
						/* >>chng 04 jan 25, changed 1.e-5 -> LIM1, LIM2, PvH */
						if( f > LIM2 )
							bfac = exp(f) - 1.;
						else if( f > LIM1 )
							bfac = (1. + 0.5*f)*f;
						else
							bfac = f;
						xx = qprob[j]*gv.bin[nd]->dstab1[i]*rfield.anu2[i]*
							rfield.widflx[i]/bfac;
						flux += xx;
					}
					else
					{
						xx = 0.;
					}
				}
			}
			else
			{
				f = TE1RYD/gv.bin[nd]->tedust*rfield.anu[i];
				if( f < x )
				{
					/* >>chng 04 jan 25, changed 1.e-5 -> LIM1, LIM2, PvH */
					if( f > LIM2 )
						bfac = exp(f) - 1.;
					else if( f > LIM1 )
						bfac = (1. + 0.5*f)*f;
					else
						bfac = f;
					flux = gv.bin[nd]->dstab1[i]*rfield.anu2[i]*rfield.widflx[i]/bfac;
				}
			}

			/* >>chng 04 jan 25, moved these multiplications outside loop, PvH */
			flux *= factor*gv.bin[nd]->cnv_H_pCM3;

			/* >>chng 96 jul 13, save grain emission array
			 * >>chng 97 feb 28, GrainEmission was set to last evaluated grain,
			 * zeroed out above, and now increment grain instead of setting it */

			/* only local emission since zeroed out on each zone */
			gv.GrainEmission[i] += (float)flux;
			/* total emission */
			grn[i] += (float)(flux*radius.dVolOutwrd);
		}
	}

	free( qprob );
	free( qtemp );

#	ifndef NDEBUG
	/*********************************************************************************
	 *
	 * Following are three checks on energy and charge conservation by the grain code.
	 * Their primary function is as an internal consistency check, so that coding
	 * errors get caught as early as possible. This is why the program terminates as
	 * soon as any one of the checks fails.
	 *
	 * NB NB - despite appearances, these checks do NOT guarantee overall energy
	 *         conservation in the Cloudy model to the asserted tolerance, see note 1B !
	 *
	 * Note 1: there are two sources for energy imbalance in the grain code (see A & B).
	 *   A: Interpolation in dstems. The code calculates how much energy the grains
	 *      emit in thermal radiation (gv.bin[nd]->GrainHeat), and converts that into
	 *      an (average) grain temperature by reverse interpolation in dstems. If
	 *      quantum heating is not used, that temperature is used directly to generate
	 *      the local diffuse emission. Hence the finite resolution of the dstems grid
	 *      can lead to small errors in flux. This is tested in Check 1. The maximum
	 *      error of interpolating in dstems scales with NDEMS^-3. The same problem
	 *      can also occur when quantum heating is used, however, the fact that many
	 *      bins are used will probably lead to a cancellation effect of the errors.
	 *   B: RT_OTS_Update gets called AFTER grain() has completed, so the grain heating
	 *      was calculated using a different version of the OTS fields than the one
	 *      that gets fed into the RT routines (where the actual attenuation of the
	 *      radiation fields by the grain opacity is done). This can lead to an energy
	 *      imbalance, depending on how accurate the convergence of the OTS fields is.
	 *      This is outside the control of the grain code and is therefore NOT checked.
	 *      Rather, the grain code remembers the contribution from the old OTS fields
	 *      (through gv.bin[nd]->BolFlux) and uses that in Check 3. In most models the
	 *      difference will be well below 0.1%, but in AGN type models where OTS continua
	 *      are important, the energy imbalance can be of the order of 0.5% of the grain
	 *      heating (status nov 2001). On 04 jan 25 the initialization of phiTilde has
	 *      been moved to qheat, implying that phiTilde now uses the updated version of
	 *      the OTS fields. The total amount of radiated energy however is still based
	 *      on gv.bin[nd]->GrainHeat which uses the old version of the OTS fields.
	 *   C: Energy conservation for collisional processes is guaranteed by adding in
	 *      (very small) correction terms. These corrections are needed to cover up
	 *      small imperfection in the theory, and cannot be avoided without making the
	 *      already very complex theory even more complex.
	 *   D: Photo-electric heating and collisional cooling can have an important effect
	 *      on the total heating balance of the gas. Both processes depend strongly on
	 *      the grain charge, so assuring proper charge balance is important as well.
	 *      This is tested in Check 2.
	 *
	 * Note 2: for quantum heating it is important to resolve the Maxwell distribution
	 *   of the electrons and ions far enough into the tail that the total amount of
	 *   energy contained in the remaining tail is negligible. If this is not the case,
	 *   the assert at the beginning of the qheat() routine will fail. This is because
	 *   the code filling up the phiTilde array in GrainCollHeating() assumes a value for
	 *   the average particle energy based on a Maxwell distribution going to infinity.
	 *   If the maximum energy used is too low, the assumed average energy would be
	 *   incorrect.
	 *
	 *********************************************************************************/

	lgNoTdustFailures = TRUE;
	for( nd=0; nd < gv.nBin; nd++ )
	{
		if( !gv.bin[nd]->lgTdustConverged )
		{
			lgNoTdustFailures = FALSE;
			break;
		}
	}

	/* CHECK 1: does the grain thermal emission conserve energy ? */
	BolFlux = 0.;
	for( i=0; i < rfield.nflux; i++ )
	{
		BolFlux += gv.GrainEmission[i]*rfield.anu[i]*EN1RYD;
	}
	Comparison1 = 0.;
	for( nd=0; nd < gv.nBin; nd++ )
	{
		if( gv.bin[nd]->tedust < gv.bin[nd]->Tsublimat )
			Comparison1 += CONSERV_TOL*gv.bin[nd]->GrainHeat;
		else
			/* for high temperatures the interpolation in dstems
			 * is less accurate, so we have to be more lenient */
			Comparison1 += 10.*CONSERV_TOL*gv.bin[nd]->GrainHeat;
	}

	/* >>chng 04 mar 11, add constant grain temperature to pass assert */
	/* >>chng 04 jun 01, deleted test for constant grain temperature, PvH */
	ASSERT( fabs(BolFlux-gv.GrainHeatSum) < Comparison1 );

	/* CHECK 2: assert charging balance */
	for( nd=0; nd < gv.nBin; nd++ )
	{
		if( gv.bin[nd]->lgChrgConverged )
		{
			double ave = 0.5*(gv.bin[nd]->RateDn+gv.bin[nd]->RateUp);
			/* >>chng 04 dec 16, add lgAbort - do not throw assert if we are in 
			 * process of aborting */
			ASSERT( lgAbort || fabs(gv.bin[nd]->RateDn-gv.bin[nd]->RateUp) < CONSERV_TOL*ave );
		}
	}

	if( lgNoTdustFailures && gv.lgDHetOn && gv.lgDColOn && thermal.ConstGrainTemp == 0. )
	{
		/* CHECK 3: calculate the total energy donated to grains, must be balanced by
		 * the energy emitted in thermal radiation plus various forms of gas heating */
		Comparison1 = 0.;
		for( nd=0; nd < gv.nBin; nd++ )
		{
			Comparison1 += gv.bin[nd]->BolFlux;
		}
		/* add in collisional heating of grains by plasma (if positive) */
		Comparison1 += MAX2(gv.GasCoolColl,0.);
		/* add in net amount of chemical energy donated by recombining ions and molecule formation */
		Comparison1 += gv.GrainHeatChem;

		/*              thermal emis        PE effect          gas heating by coll    thermionic emis */
		Comparison2 = gv.GrainHeatSum+thermal.heating[0][13]+thermal.heating[0][14]+thermal.heating[0][25];

		ASSERT( lgAbort || fabs(Comparison1-Comparison2)/Comparison2 < CONSERV_TOL || gv.lgBakesPAH_heat );
	}
#	endif

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


/****************************************************************************
 *
 * qheat: driver routine for non-equilibrium heating of grains
 *
 * This routine calls GetProbDistr_LowLimit, GetProbDistr_HighLimit
 * (which do the actual non-equilibrium calculations), and does the
 * subsequent quality control.
 *
 * Written by Peter van Hoof (UK, CITA, QUB).
 *
 ****************************************************************************/

/* this is the new version of the quantum heating code, used starting Cloudy 96 beta 3 */

void qheat(/*@out@*/ double qtemp[],  /* qtemp[NQGRID] */
	   /*@out@*/ double qprob[],  /* qprob[NQGRID] */
	   /*@out@*/ long int *qnbin,
	   long int nd)
{
	int lgBoundErr,
	  lgDelta,
	  lgNegRate,
	  lgOK,
	  lgTried;
	long int i,
	  nWideFail;
	QH_Code ErrorCode,
	  ErrorCode2,
	  ErrorStart;
	double c0,
	  c1,
	  c2,
	  check,
	  DefFac,
	  deriv,
	  *dPdlnT/*[NQGRID]*/,
	  fwhm,
	  FwhmRatio,
	  integral,
	  minBracket,
	  maxBracket,
	  new_tmin,
	  NumEvents,
	  oldy,
	  *Phi/*[NC_ELL]*/,
	  *PhiDrv/*[NC_ELL]*/,
	  *phiTilde/*[NC_ELL]*/,
	  rel_tol,
	  Tmax,
	  tol,
	  Umax,
	  xx,
	  y;


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

	/* sanity checks */
	ASSERT( gv.bin[nd]->lgQHeat );
	ASSERT( nd >= 0 && nd < gv.nBin );

	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "\n >>>>qheat called for grain %s\n", gv.bin[nd]->chDstLab );
	}

	/* >>chng 01 aug 22, allocate space */
	/* phiTilde is continuum corrected for photo-electric effect, in events/H/s/cell, default depl */
	if( (phiTilde = (double*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(double)) ) ) == NULL )
		BadMalloc();
	if( (Phi = (double*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(double)) ) ) == NULL )
		BadMalloc();
	if( (PhiDrv = (double*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(double)) ) ) == NULL )
		BadMalloc();
	if( (dPdlnT = (double*)MALLOC((size_t)(NQGRID*sizeof(double)) ) ) == NULL )
		BadMalloc();

	qheat_init( nd, phiTilde, &check );

	check += gv.bin[nd]->GrainHeatColl-gv.bin[nd]->GrainCoolTherm;

	xx = integral = 0.;
	c0 = c1 = c2 = 0.;
	lgNegRate = FALSE;
	oldy = 0.;
	tol = 1.;

	/* >>chng 01 nov 29, rfield.nflux -> gv.qnflux, PvH */
	/* >>chng 03 jan 26, gv.qnflux -> gv.bin[nd]->qnflux, PvH */
	for( i=gv.bin[nd]->qnflux-1; i >= 0; i-- )
	{
		/* >>chng 97 jul 17, to summed continuum */
		/* >>chng 00 mar 30, to phiTilde, to keep track of photo-electric effect and collisions, by PvH */
		/* >>chng 01 oct 10, use trapezoidal rule for integrating Phi, reverse direction of integration.
		 * Phi[i] is now integral from exactly rfield.anu[i] to infinity to second-order precision, PvH */
		/* >>chng 01 oct 30, change normalization of Phi, PhiDrv from <unit>/cm^3 -> <unit>/grain, PvH */
		double fac = ( i < gv.bin[nd]->qnflux-1 ) ? 1./rfield.widflx[i] : 0.;
		/* phiTilde has units events/H/s, y has units events/grain/s/Ryd */
		y = phiTilde[i]*gv.bin[nd]->cnv_H_pGR*fac;
		/* PhiDrv[i] = (Phi[i+1]-Phi[i])/(anu[i+1]-anu[i]), units events/grain/s/Ryd */
		PhiDrv[i] = -0.5*(y + oldy);
		/* there is a minus sign here because we are integrating from infinity downwards */
		xx -= PhiDrv[i]*(rfield.anu[i+1]-rfield.anu[i]);
		/* Phi has units events/grain/s */
		Phi[i] = xx;

#		ifndef NDEBUG
		/* trapezoidal rule is not needed for integral, this is also second-order correct */
		integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu[i]*EN1RYD;
#		endif

		/* c<n> has units Ryd^(n+1)/grain/s */
		c0 += Phi[i]*rfield.widflx[i];
		c1 += Phi[i]*rfield.anu[i]*rfield.widflx[i];
		c2 += Phi[i]*POW2(rfield.anu[i])*rfield.widflx[i];

		lgNegRate = lgNegRate || ( phiTilde[i] < 0. );

		oldy = y;
	}

	/* sanity check */
	ASSERT( fabs(check-integral)/check <= CONSERV_TOL );

#	if 0
	{
		char fnam[50];
		FILE *file;

		sprintf(fnam,"Lambda_%2.2ld.asc",nd);
		file = fopen(fnam,"w");
		for( i=0; i < NDEMS; ++i )
			fprintf(file,"%e %e %e\n",
				exp(gv.dsttmp[i]),
				ufunct(exp(gv.dsttmp[i]),nd,&lgBoundErr),
				exp(gv.bin[nd]->dstems[i])*gv.bin[nd]->cnv_H_pGR/EN1RYD);
		fclose(file);

		sprintf(fnam,"Phi_%2.2ld.asc",nd);
		file = fopen(fnam,"w");
		for( i=0; i < gv.bin[nd]->qnflux; ++i )
			fprintf(file,"%e %e\n", rfield.anu[i],Phi[i]);
		fclose(file);
	}
#	endif

	/* Tmax is where p(U) will peak (at least in high intensity limit) */
	Tmax = gv.bin[nd]->tedust;
	/* grain enthalpy at peak of p(U), in Ryd */
	Umax = ufunct(Tmax,nd,&lgBoundErr);
	ASSERT( !lgBoundErr ); /* this should never happen */
	/* y is dln(Lambda)/dlnT at Tmax */
	spldrv_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,log(Tmax),&y,&lgBoundErr);
	ASSERT( !lgBoundErr ); /* this should never happen */
	/* deriv is dLambda/dU at Umax, in 1/grain/s */
	deriv = y*c0/(uderiv(Tmax,nd)*Tmax);
	/* estimated FWHM of probability distribution, in Ryd */
	fwhm = sqrt(8.*LN_TWO*c1/deriv);

	NumEvents = POW2(fwhm)*c0/(4.*LN_TWO*c2);
	FwhmRatio = fwhm/Umax;

	/* >>chng 01 nov 15, change ( NumEvents > MAX_EVENTS2 ) --> ( FwhmRatio < FWHM_RATIO2 ), PvH */
	lgDelta = ( FwhmRatio < FWHM_RATIO2 );
	/* high intensity case is always OK since we will use equilibrium treatment */
	lgOK = lgDelta;

	ErrorStart = QH_OK;

	if( lgDelta ) 
	{
		/* in this case we ignore negative rates, equilibrium treatment is good enough */
		lgNegRate = FALSE;
		ErrorStart = MAX2(ErrorStart,QH_DELTA);
	}

	if( lgNegRate )
		ErrorStart = MAX2(ErrorStart,QH_NEGRATE_FAIL);

	ErrorCode = ErrorStart;

	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "   grain heating: %.4e, integral %.4e, total rate %.4e\n",
			 gv.bin[nd]->GrainHeat,integral,Phi[0]);
		fprintf( ioQQQ, "   av grain temp %.4e av grain enthalpy (Ryd) %.4e\n",
			 gv.bin[nd]->tedust,Umax);
		fprintf( ioQQQ, "   fwhm^2/(4ln2*c2/c0): %.4e fwhm (Ryd) %.4e fwhm/Umax %.4e\n",
			 NumEvents,fwhm,FwhmRatio );
		fprintf( ioQQQ, "   lgQHTooWide %d\n", gv.bin[nd]->lgQHTooWide );
	}

	/* these two variables will bracket qtmin, they should only be needed during the initial search phase */
	minBracket = GRAIN_TMIN;
	maxBracket = gv.bin[nd]->tedust;

	/* >>chng 02 jan 30, introduced lgTried to avoid running GetProbDistr_HighLimit over and over..., PvH */
	lgTried = FALSE;
	/* >>chng 02 aug 06, introduced QH_WIDE_FAIL and nWideFail, PvH */
	nWideFail = 0;
	/* >>chng 03 jan 27, introduced DefFac to increase factor for repeated LOW_FAIL's, PvH */
	DefFac = DEF_FAC;
	/* >>chng 04 nov 10, introduced rel_tol to increase precision in case of repeated CONV_FAIL's, PvH */
	rel_tol = 1.;

	/* if average number of multiple photon events is too high, lgOK is set to TRUE */
	/* >>chng 02 aug 12, added gv.bin[nd]->lgQHTooWide to prevent unnecessary looping here.
	 * In general the number of integration steps needed to integrate the probability distribution
	 * will increase monotonically with depth into the cloud. Hence, once the distribution becomes
	 * too wide to fit into NQGRID steps (which only happens for extremely cold grains in deeply
	 * shielded conditions) there is no hope of ever converging GetProbDistr_LowLimit and the code
	 * will waste a lot of CPU time establishing this for every zone again. So once the distribution
	 * becomes too wide we immediately skip to the analytic approximation to save time, PvH */
	for( i=0; i < LOOP_MAX && !lgOK && !gv.bin[nd]->lgQHTooWide; i++ )
	{
		if( gv.bin[nd]->qtmin >= gv.bin[nd]->tedust )
		{
			/* >>chng 02 jul 31, was gv.bin[nd]->qtmin = 0.7*gv.bin[nd]->tedust */
			/* >>chng 03 nov 10, changed Umax/exp(+... to Umax*exp(-... to avoid overflow, PvH */
			double Ulo = Umax*exp(-sqrt(-log(PROB_CUTOFF_LO)/(4.*LN_TWO))*fwhm/Umax);
			double MinEnth = exp(gv.bin[nd]->DustEnth[0]);
			Ulo = MAX2(Ulo,MinEnth);
			gv.bin[nd]->qtmin = inv_ufunct(Ulo,nd,&lgBoundErr);
			ASSERT( !lgBoundErr ); /* this should never happen */
			/* >>chng 02 jul 30, added this test; prevents problems with ASSERT below, PvH */
			if( gv.bin[nd]->qtmin <= minBracket || gv.bin[nd]->qtmin >= maxBracket )
				gv.bin[nd]->qtmin = sqrt(minBracket*maxBracket);
		}
		gv.bin[nd]->qtmin = MAX2(gv.bin[nd]->qtmin,GRAIN_TMIN);

		ASSERT( minBracket <= gv.bin[nd]->qtmin && gv.bin[nd]->qtmin <= maxBracket );

		ErrorCode = ErrorStart;

		/* >>chng 01 nov 15, add ( FwhmRatio >= FWHM_RATIO ), PvH */
		if( FwhmRatio >= FWHM_RATIO && NumEvents <= MAX_EVENTS )
		{
			GetProbDistr_LowLimit(nd,rel_tol,Umax,fwhm,Phi,PhiDrv,qtemp,qprob,dPdlnT,qnbin,
					      &new_tmin,&nWideFail,&ErrorCode);

			/* >>chng 02 jan 07, various changes to improve convergence for very cold grains, PvH */
			if( ErrorCode == QH_DELTA_FAIL && fwhm < Umax && !lgTried )
			{
				double dummy;

				/* this situation can mean two things: either the photon rate is so high that
				 * the code needs too many steps to integrate the probability distribution,
				 * or alternatively, tmin is far too low and the code needs too many steps
				 * to overcome the rising side of the probability distribution.
				 * So we call GetProbDistr_HighLimit first to determine if the former is the
				 * case; if that fails then the latter must be true and we reset QH_DELTA_FAIL */
				ErrorCode = MAX2(ErrorStart,QH_ANALYTIC);
				/* use dummy to avoid losing estimate for new_tmin from GetProbDistr_LowLimit */
				/* >>chng 02 aug 06, introduced STRICT and RELAX, PvH */
				GetProbDistr_HighLimit(nd,STRICT,Umax,fwhm,qtemp,qprob,dPdlnT,&tol,qnbin,&dummy,
						       &ErrorCode);
						       
				if( ErrorCode >= QH_RETRY )
				{
					ErrorCode = QH_DELTA_FAIL;
					lgTried = TRUE;
				}
			}

			/* >>chng 02 aug 07 major rewrite of the logic below */
			if( ErrorCode < QH_NO_REBIN )
			{
				if( new_tmin < minBracket || new_tmin > maxBracket )
					++nWideFail;

				if( nWideFail < WIDE_FAIL_MAX )
				{
					if( new_tmin <= minBracket )
						new_tmin = sqrt(gv.bin[nd]->qtmin*minBracket);
					if( new_tmin >= maxBracket )
						new_tmin = sqrt(gv.bin[nd]->qtmin*maxBracket);
				}
				else
				{
					ErrorCode = MAX2(ErrorCode,QH_WIDE_FAIL);
				}

				if( ErrorCode == QH_CONV_FAIL )
				{
					rel_tol *= 0.9;
				}
			}
			else if( ErrorCode == QH_LOW_FAIL )
			{
				double help1 = gv.bin[nd]->qtmin*sqrt(DefFac);
				double help2 = sqrt(gv.bin[nd]->qtmin*maxBracket);
				minBracket = gv.bin[nd]->qtmin;
				new_tmin = MIN2(help1,help2);
				/* increase factor in case we get repeated LOW_FAIL's */
				DefFac += 1.5;
			}
			else if( ErrorCode == QH_HIGH_FAIL )
			{
				double help = sqrt(gv.bin[nd]->qtmin*minBracket);
				maxBracket = gv.bin[nd]->qtmin;
				new_tmin = MAX2(gv.bin[nd]->qtmin/DEF_FAC,help);
			}
			else
			{
				new_tmin = sqrt(minBracket*maxBracket);
			}
		}
		else
		{
			GetProbDistr_HighLimit(nd,STRICT,Umax,fwhm,qtemp,qprob,dPdlnT,&tol,qnbin,&new_tmin,
					       &ErrorCode);
		}

		gv.bin[nd]->qtmin = new_tmin;

		lgOK = ( ErrorCode < QH_RETRY );

		if( ErrorCode >= QH_FATAL )
			break;

		if( ErrorCode != QH_LOW_FAIL )
			DefFac = DEF_FAC;

		if( trace.lgTrace && trace.lgDustBug ) 
		{
			fprintf( ioQQQ, "   GetProbDistr returns code %d\n", ErrorCode );
			if( !lgOK )
			{
				fprintf( ioQQQ, " >>>>qheat trying another iteration, qtmin bracket %.4e %.4e",
					 minBracket,maxBracket );
				fprintf( ioQQQ, " nWideFail %ld\n", nWideFail );
			}
		}
	}

	if( ErrorCode == QH_WIDE_FAIL )
		gv.bin[nd]->lgQHTooWide = TRUE;

	/* >>chng 03 jan 17, added test for !lgDelta, PvH */
	/* if( gv.bin[nd]->lgQHTooWide ) */
	if( gv.bin[nd]->lgQHTooWide && !lgDelta )
		ErrorCode = MAX2(ErrorCode,QH_WIDE_FAIL);

/* 	if( ErrorCode >= QH_RETRY ) */
/* 		printf( "      *** PROBLEM  loop not converged, errorcode %d\n",ErrorCode ); */

	/* The quantum heating code tends to run into trouble when it goes deep into the neutral zone,
	 * especially if the original spectrum was very hard, as is the case in high excitation PNe or AGN.
	 * You then get a bipartition in the spectrum where most of the photons have low energy, while
	 * there still are hard X-ray photons left. The fact that the average energy per photon is low
	 * forces the quantum code to take tiny little steps when integrating the probability distribution,
	 * while the fact that X-ray photons are still present implies that big temperature spikes still
	 * occur and hence the temperature distribution is very wide. Therefore the code needs a zillion
	 * steps to integrate the probability distribution and simply runs out of room. As a last resort
	 * try the analytic approximation with relaxed constraints used below. */
	/* >>chng 02 oct 03, vary Umax and fwhm to force fit with fwhm/Umax remaining constant */
	/* >>chng 03 jan 17, changed test so that last resort always gets tried when lgOK = lgDelta = FALSE, PvH */
	/* if( !lgOK && FwhmRatio >= FWHM_RATIO && NumEvents <= MAX_EVENTS ) */
	if( !lgOK && !lgDelta )
	{
		double Umax2 = Umax*sqrt(tol);
		double fwhm2 = fwhm*sqrt(tol);

		for( i=0; i < LOOP_MAX; ++i )
		{
			double dummy;

			ErrorCode2 = MAX2(ErrorStart,QH_ANALYTIC);
			GetProbDistr_HighLimit(nd,RELAX,Umax2,fwhm2,qtemp,qprob,dPdlnT,&tol,qnbin,&dummy,
					       &ErrorCode2);

			lgOK = ( ErrorCode2 < QH_RETRY );
			if( lgOK )
			{
				gv.bin[nd]->qtmin = dummy;
				ErrorCode = ErrorCode2;
				break;
			}
			else
			{
				Umax2 *= sqrt(tol);
				fwhm2 *= sqrt(tol);
			}
		}
	}

	if( nzone == 1 )
		gv.bin[nd]->qtmin_zone1 = gv.bin[nd]->qtmin;

	gv.bin[nd]->lgUseQHeat = ( lgOK && !lgDelta );
	gv.bin[nd]->lgEverQHeat = ( gv.bin[nd]->lgEverQHeat || gv.bin[nd]->lgUseQHeat );

	if( lgOK )
	{
		if( trace.lgTrace && trace.lgDustBug )
			fprintf( ioQQQ, " >>>>qheat converged with code: %d\n", ErrorCode );
	}
	else
	{
		*qnbin = 0;
		++gv.bin[nd]->QHeatFailures;
		fprintf( ioQQQ, " PROBLEM  qheat did not converge grain %s in zone %ld, error code = %d\n",
			 gv.bin[nd]->chDstLab,nzone,ErrorCode );		
	}

	if( gv.QHPunchFile != NULL && ( iterations.lgLastIt || !gv.lgQHPunLast ) ) 
	{
		fprintf( gv.QHPunchFile, "\nDust Temperature Distribution: grain %s zone %ld\n",
			 gv.bin[nd]->chDstLab,nzone );

		fprintf( gv.QHPunchFile, "Equilibrium temperature: %.2f\n", gv.bin[nd]->tedust );

		if( gv.bin[nd]->lgUseQHeat ) 
		{
			/* >>chng 01 oct 09, remove qprob from output, it depends on step size, PvH */
			fprintf( gv.QHPunchFile, "Number of bins: %ld\n", *qnbin );
			fprintf( gv.QHPunchFile, "  Tgrain      dP/dlnT\n" );
			for( i=0; i < *qnbin; i++ ) 
			{
				fprintf( gv.QHPunchFile, "%.4e %.4e\n", qtemp[i],dPdlnT[i] );
			}
		}
		else 
		{
			fprintf( gv.QHPunchFile, "**** quantum heating was not used\n" );
		}
	}

	free ( phiTilde );
	free ( Phi );
	free ( PhiDrv );
	free ( dPdlnT );
 
#	ifdef DEBUG_FUN
	fputs( " <->qheat()\n", debug_fp );
#	endif
	return;
}


/* initialize phiTilde */
STATIC void qheat_init(long nd,
		       /*@out@*/ double phiTilde[],  /* phiTilde[rfield.nupper] */
		       /*@out@*/ double *check)
{
	long i,
	  nz;
	double sum = 0.;

	/*@-redef@*/
	enum {DEBUG_LOC=FALSE};
	/*@+redef@*/

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

	/* sanity checks */
	ASSERT( gv.bin[nd]->lgQHeat );
	ASSERT( nd >= 0 && nd < gv.nBin );

	*check = 0.;

	/* >>chng 01 nov 29, rfield.nflux -> gv.qnflux, PvH */
	/* >>chng 03 jan 26, gv.qnflux -> gv.bin[nd]->qnflux, PvH */
	for( i=0; i < gv.bin[nd]->qnflux; i++ )
	{
		phiTilde[i] = 0.;
	}

	/* fill in auxiliary array for quantum heating routine
	 * it reshuffles the input spectrum times the cross section to take
	 * the photo-electric effect into account. this prevents the quantum
	 * heating routine from having to calculate this effect over and over
	 * again; it can do a straight integration instead, making the code
	 * a lot simpler and faster. this initializes the array for non-ionizing
	 * energies, the reshuffling for higher energies is done in the next loop
	 * phiTilde has units events/H/s/cell at default depletion */

	for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
	{
		double check1 = 0.;
		ChargeBin *gptr = gv.bin[nd]->chrg[nz];

		/* integrate over incident continuum for non-ionizing energies */
		for( i=0; i < MIN2(gptr->ipThresInf,rfield.nflux); i++ )
		{
#			ifndef NDEBUG
			check1 += rfield.SummedCon[i]*gv.bin[nd]->dstab1[i]*rfield.anu[i];
#			endif

			phiTilde[i] += gptr->FracPop*rfield.SummedCon[i]*gv.bin[nd]->dstab1[i];
		}

		/* >>chng 01 mar 02, use new expresssions for grain cooling and absorption
		 * cross sections following the discussion with Joe Weingartner, PvH */
		for( i=gptr->ipThresInf; i < rfield.nflux; i++ )
		{
			double cs1,cs2,cs_tot,cool1,cool2,ehat1,ehat2;

			PE_init(nd,nz,i,&cs1,&cs2,&cs_tot,&cool1,&cool2,&ehat1,&ehat2);

#			ifndef NDEBUG
			check1 += rfield.SummedCon[i]*MAX2(cs_tot*rfield.anu[i]-cs1*cool1-cs2*cool2,0.);
#			endif

			/* this accounts for the photons that are fully absorbed by grain */
			phiTilde[i] += gptr->FracPop*rfield.SummedCon[i]*MAX2(gv.bin[nd]->dstab1[i]-cs1,0.);

			/* >>chng 01 oct 10, use bisection search to find ip. On C scale now */

			/* this accounts for photons that eject an electron from the valence band */
			if( cs1 > 0. )
			{
				/* we treat photo-ejection from the valence band and de-excitation
				 * from the conduction band as one simultaneous event */
				long ipLo = 0;
				long ipHi = i;
				float xx = rfield.anu[i] - (float)cool1;
				/* find anu[ipLo]-0.5*widflx[ipLo] <= xx < anu[ipHi]-0.5*widflx[ipHi] */
				while( ipHi-ipLo > 1 )
				{
					long ipMd = (ipLo+ipHi)/2;
					if( gv.anumin[ipMd] > xx )
						ipHi = ipMd;
					else
						ipLo = ipMd;
				}
				phiTilde[ipLo] += gptr->FracPop*rfield.SummedCon[i]*cs1;
			}

			/* no need to account for photons that eject an electron from the conduction band */
			/* >>chng 01 dec 11, cool2 always equal to rfield.anu[i] -> no grain heating */
		}

		*check += gptr->FracPop*check1*EN1RYD*gv.bin[nd]->cnv_H_pCM3;

		if( DEBUG_LOC )
		{
			double integral = 0.;
			for( i=0; i < gv.bin[nd]->qnflux; i++ )
			{
				integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu[i]*EN1RYD;
			}
			sum += gptr->FracPop*check1*EN1RYD*gv.bin[nd]->cnv_H_pCM3;
			printf(" integral test 1: integral %.6e %.6e\n",integral,sum);
		}

		/* add quantum heating due to recombination of electrons, subtract thermionic cooling */

		/* gptr->HeatingRate2 is net heating rate in erg/H/s at standard depl
		 * includes contributions for recombining electrons, autoionizing electrons
		 * subtracted by thermionic emissions here since it is inverse process
		 *
		 * NB - in extreme conditions this rate may be negative (if there
		 * is an intense radiation field leading to very hot grains, but no ionizing
		 * photons, hence very few free electrons). we assume that the photon rates
		 * are high enough under those circumstances to avoid phiTilde becoming negative,
		 * but we will check that in qheat1 anyway. */

		/* >>chng 03 nov 06, check for extremely low HeatingRate and save CPU time, pah_crash.in, PvH */
		if( fabs(gptr->HeatingRate2*gv.bin[nd]->cnv_H_pCM3) > 0.05*CONSERV_TOL*gv.bin[nd]->GrainHeat ) 
		{
			double *RateArr,Sum,ESum,DSum,Delta,E_av2,Corr;
			double fac = BOLTZMANN/EN1RYD*phycon.te;
			/* E0 is barrier that electron needs to overcome, zero for positive grains */
			/* >>chng 03 jan 23, added second term to correctly account for autoionizing states
			 *                   where ThresInfInc is negative, tested in small_grain.in, PvH */
			double E0 = -(MIN2(gptr->PotSurfInc,0.) + MIN2(gptr->ThresInfInc,0.));
			/* >>chng 01 mar 02, this should be energy gap between top electron and infinity, PvH */
			/* >>chng 01 nov 21, use correct barrier: ThresInf[nz] -> ThresInfInc[nz], PvH */
			/* >>chng 03 jan 23, replaced -E0 with MIN2(PotSurfInc[nz],0.), PvH */
			double Einf = gptr->ThresInfInc + MIN2(gptr->PotSurfInc,0.);
			/* this is average energy deposited by one event, in erg
			 * this value is derived from distribution assumed here, and may
			 * not be the same as HeatElectrons/(CollisionRateElectr*eta) !! */
			/* >>chng 01 nov 21, use correct barrier: ThresInf[nz] -> ThresInfInc[nz], PvH */
			/* >>chng 03 jan 23, replaced ThresInfInc[nz] with MAX2(ThresInfInc[nz],0.), PvH */
			double E_av = MAX2(gptr->ThresInfInc,0.)*EN1RYD + 2.*BOLTZMANN*phycon.te;
			/* this is rate in events/H/s at standard depletion */
			double rate = gptr->HeatingRate2/E_av;

			double ylo = -exp(-E0/fac);
			/* this is highest kinetic energy of electron that can be represented in phiTilde */
			/* >>chng 01 nov 29, rfield.nflux -> gv.qnflux, PvH */
			/* >>chng 03 jan 26, gv.qnflux -> gv.bin[nd]->qnflux, PvH */
			double Ehi = gv.anumax[gv.bin[nd]->qnflux-1]-Einf;
			double yhi = ((E0-Ehi)/fac-1.)*exp(-Ehi/fac);
			/* renormalize rate so that integral over phiTilde*anu gives correct total energy */
			rate /= yhi-ylo;

			/* correct for fractional population of this charge state */
			rate *= gptr->FracPop;

			/* >>chng 03 jan 24, add code to correct for discretization errors, hotdust.in, PvH */
			if( (RateArr=(double*)MALLOC((size_t)((unsigned)gv.bin[nd]->qnflux*sizeof(double)))) == NULL )
				BadMalloc();
			Sum = ESum = DSum = 0.;

			/* >>chng 04 jan 21, replaced gv.bin[nd]->qnflux -> gv.bin[nd]->qnflux2, PvH */
			for( i=0; i < gv.bin[nd]->qnflux2; i++ ) 
			{
				Ehi = gv.anumax[i] - Einf;
				if( Ehi >= E0 ) 
				{
					/* Ehi is kinetic energy of electron at infinity */
					yhi = ((E0-Ehi)/fac-1.)*exp(-Ehi/fac);
					/* >>chng 01 mar 24, use MAX2 to protect against roundoff error, PvH */
					RateArr[i] = rate*MAX2(yhi-ylo,0.);
					Sum += RateArr[i];
					ESum += rfield.anu[i]*RateArr[i];
#					ifndef NDEBUG
					DSum += rfield.widflx[i]*RateArr[i];
#					endif
					ylo = yhi;
				}
				else
				{
					RateArr[i] = 0.;
				}
			}
			E_av2 = ESum/Sum*EN1RYD;
			Delta = DSum/Sum*EN1RYD;
			ASSERT( fabs(E_av-E_av2) <= Delta );
			Corr = E_av/E_av2;

			for( i=0; i < gv.bin[nd]->qnflux2; i++ ) 
			{
				phiTilde[i] += RateArr[i]*Corr;
			}

			if( DEBUG_LOC )
			{
				double integral = 0.;
				for( i=0; i < gv.bin[nd]->qnflux; i++ )
				{
					integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu[i]*EN1RYD;
				}
				sum += gptr->FracPop*gptr->HeatingRate2;
				printf(" integral test 2: integral %.6e %.6e\n",integral,sum);
			}

			free( RateArr );
		}
	}

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

	/* add quantum heating due to molecule/ion collisions */

	/* gv.bin[nd]->HeatingRate1 is heating rate in erg/H/s at standard depl
	 * includes contributions from molecules/neutral atoms and recombining ions
	 *
	 * in fully ionized conditions electron heating rates will be much higher
	 * than ion and molecule rates since electrons are so much faster and grains
	 * tend to be positive. in non-ionized conditions the main contribution will
	 * come from neutral atoms and molecules, so it is appropriate to treat both
	 * the same. in fully ionized conditions we don't care since unimportant.
	 *
	 * NB - if grains are hotter than ambient gas, the heating rate may become negative.
	 * if photon rates are not high enough to prevent phiTilde from becoming negative,
	 * we will raise a flag while calculating the quantum heating in qheat1 */

	/* >>chng 03 nov 06, check for extremely low HeatingRate and save CPU time, PvH */
	if( fabs(gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3) > 0.05*CONSERV_TOL*gv.bin[nd]->GrainHeat )
	{
		/* limits for Taylor expansion of (1+x)*exp(-x) */
		/* these choices will assure only 6 digits precision */
		const double LIM2 = pow(3.e-6,1./3.);
		const double LIM3 = pow(8.e-6,1./4.);
		/* if gas temperature is higher than grain temperature we will
		 * consider Maxwell-Boltzmann distribution of incoming particles
		 * and ignore distribution of outgoing particles, if grains
		 * are hotter than ambient gas, we use reverse treatment */
		double fac = BOLTZMANN/EN1RYD*MAX2(phycon.te,gv.bin[nd]->tedust);
		/* this is average energy deposited/extracted by one event, in erg */
		double E_av = 2.*BOLTZMANN*MAX2(phycon.te,gv.bin[nd]->tedust);
		/* this is rate in events/H/s at standard depletion */
		double rate = gv.bin[nd]->HeatingRate1/E_av;

		double ylo = -1.;
		/* this is highest energy of incoming/outgoing particle that can be represented in phiTilde */
		/* >>chng 01 nov 29, rfield.nflux -> gv.qnflux, PvH */
		/* >>chng 03 jan 26, gv.qnflux -> gv.bin[nd]->qnflux, PvH */
		double Ehi = gv.anumax[gv.bin[nd]->qnflux-1];
		double yhi = -(Ehi/fac+1.)*exp(-Ehi/fac);
		/* renormalize rate so that integral over phiTilde*anu gives correct total energy */
		rate /= yhi-ylo;

		for( i=0; i < gv.bin[nd]->qnflux2; i++ ) 
		{
			/* Ehi is kinetic energy of incoming/outgoing particle
			 * we assume that Ehi-E0 is deposited/extracted from grain */
			/* Ehi = gv.anumax[i]; */
			double x = gv.anumax[i]/fac;
			/* (1+x)*exp(-x) = 1 - 1/2*x^2 + 1/3*x^3 - 1/8*x^4 + O(x^5)
			 *               = 1 - Sum_n=2^infty (-x)^n/(n*(n-2)!)      */
			if( x > LIM3 )
				yhi = -(x+1.)*exp(-x);
			else if( x > LIM2 )
				yhi = -(((1./3.)*x - 0.5)*x*x + 1.);
			else
				yhi = -(1. - 0.5*x*x);
			
			/* >>chng 01 mar 24, use MAX2 to protect against roundoff error, PvH */
			phiTilde[i] += rate*MAX2(yhi-ylo,0.);
			ylo = yhi;
		}
	}

	if( DEBUG_LOC )
	{
		double integral = 0.;
		for( i=0; i < gv.bin[nd]->qnflux; i++ )
		{
			integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu[i]*EN1RYD;
		}
		sum += gv.bin[nd]->HeatingRate1;
		printf(" integral test 3: integral %.6e %.6e\n",integral,sum);
	}

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


/*******************************************************************************************
 *
 * GetProbDistr_LowLimit: main routine for calculating non-equilibrium heating of grains
 *
 * This routine implements the algorithm outlined in:
 * >>refer	grain	physics	Guhathakurtha & Draine, 1989, ApJ, 345, 230
 *
 * The original (fortran) version of the code was written by Kevin Volk.
 *
 * Heavily modified and adapted for new style grains by Peter van Hoof.
 *
 *******************************************************************************************/

STATIC void GetProbDistr_LowLimit(long int nd,
				  double rel_tol,
				  double Umax,
				  double fwhm,
				  /*@in@*/ double Phi[],     /* Phi[NQGRID]      */
				  /*@in@*/ double PhiDrv[],  /* PhiDrv[NQGRID]   */
				  /*@out@*/ double qtemp[],  /* qtemp[NQGRID]    */
				  /*@out@*/ double qprob[],  /* qprob[NQGRID]    */
				  /*@out@*/ double dPdlnT[], /* dPdlnT[NQGRID]   */
				  /*@out@*/ long int *qnbin,
				  /*@out@*/ double *new_tmin,
				  long *nWideFail,
				  /*@out@*/ QH_Code *ErrorCode)
{
	int lgAllNegSlope,
	  lgBoundErr;
	long int j,
	  k,
	  l,
	  nbad,
	  nbin,
	  nend=0,
	  nmax,
	  nok,
	  nstart=0,
	  nstart2=0;
	double dCool,
	  delu[NQGRID],
	  dlnLdlnT,
	  dlnpdlnU,
	  fac = 0.,
	  Lambda[NQGRID],
	  maxVal,
	  NextStep,
	  p[NQGRID],
	  /* >>chng 01 aug 22, from array to pointer with MALLOC */
	  qtmin1, 
	  RadCooling,
	  sum,
	  u1[NQGRID],
	  y;


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

	/* sanity checks */
	ASSERT( nd >= 0 && nd < gv.nBin );

	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "   GetProbDistr_LowLimit called for grain %s\n", gv.bin[nd]->chDstLab );
		fprintf( ioQQQ, "    got qtmin1 %.4e\n", gv.bin[nd]->qtmin);
	}

	qtmin1 = gv.bin[nd]->qtmin;
	qtemp[0] = qtmin1;
	/* u1 holds enthalpy in Ryd/grain */
	u1[0] = ufunct(qtemp[0],nd,&lgBoundErr);
	ASSERT( !lgBoundErr ); /* this should never happen */
	/* >>chng 00 mar 22, taken out factor 4, factored in hden and dstAbund
	 * interpolate in dstems array instead of integrating explicitly, by PvH */
	splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,log(qtemp[0]),&y,&lgBoundErr);
	ASSERT( !lgBoundErr ); /* this should never happen */
	/* Lambda holds the radiated energy for grains in this bin, in Ryd/s/grain */
	Lambda[0] = exp(y)*gv.bin[nd]->cnv_H_pGR/EN1RYD;
	/* set up first step of integration */
	/* >>chng 01 nov 14, set to 2.*Lambda[0]/Phi[0] instead of u1[0],
	 * this choice assures that p[1] doesn't make a large jump from p[0], PvH */
	delu[0] = 2.*Lambda[0]/Phi[0];
	p[0] = PROB_CUTOFF_LO;
	dPdlnT[0] = p[0]*qtemp[0]*uderiv(qtemp[0],nd);
	RadCooling = 0.5*p[0]*Lambda[0]*delu[0];
	NextStep = 0.01*Lambda[0]/Phi[0];
	/* >>chng 03 nov 10, added extra safeguard against stepsize too small, PvH */
	if( NextStep < 0.05*DBL_EPSILON*u1[0] )
	{
		*ErrorCode = MAX2(*ErrorCode,QH_STEP_FAIL);

#		ifdef DEBUG_FUN
		fputs( " <->GetProbDistr_LowLimit()\n", debug_fp );
#		endif
		return;
	}
	else if( NextStep < 5.*DBL_EPSILON*u1[0] )
	{
		/* try a last resort... */
		NextStep *= 100.;
	}

	nbad = 0;
	k = 0;

	*qnbin = 0;
	*new_tmin = qtmin1;
	lgAllNegSlope = TRUE;
	maxVal = dPdlnT[0];
	nmax = 0;

	/* this test neglects a negative contribution which is impossible to calculate, so it may
	 * fail to detect cases where the probability distribution starts dropping immediately,
	 * we will use a second test using the variable lgAllNegSlope below to catch those cases, PvH */
	spldrv_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,log(qtemp[0]),&dlnLdlnT,&lgBoundErr);
	ASSERT( !lgBoundErr ); /* this should never happen */
	dlnpdlnU = u1[0]*Phi[0]/Lambda[0] - dlnLdlnT*u1[0]/(qtemp[0]*uderiv(qtemp[0],nd));
	if( dlnpdlnU < 0. )
	{
		/* >>chng 03 nov 06, confirm this by integrating first step..., pah_crash.in, PvH */
		(void)TryDoubleStep(u1,delu,p,qtemp,Lambda,Phi,PhiDrv,NextStep,&dCool,k,nd,&lgBoundErr);
		dPdlnT[2] = p[2]*qtemp[2]*uderiv(qtemp[2],nd);

		if( dPdlnT[2] < dPdlnT[0] ) {
			/* if dPdlnT starts falling immediately, 
			 * qtmin1 was too high and convergence is impossible */
			*ErrorCode = MAX2(*ErrorCode,QH_HIGH_FAIL);

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

	/* NB NB -- every break in this loop should set *ErrorCode (except for regular stop criterion) !! */
	for( l=0; l < MAX_LOOP; ++l )
	{
		double RelError = TryDoubleStep(u1,delu,p,qtemp,Lambda,Phi,PhiDrv,NextStep,&dCool,k,nd,&lgBoundErr);

		/* this happens if the grain temperature in qtemp becomes higher than GRAIN_TMAX
		 * nothing that TryDoubleStep returns can be trusted, so this check should be first */
		if( lgBoundErr )
		{
			nbad += 2;
			*ErrorCode = MAX2(*ErrorCode,QH_THIGH_FAIL);
			break;
		}

		/* estimate new stepsize */
		if( RelError > rel_tol*QHEAT_TOL )
		{
			nbad += 2;

			/* step is rejected, decrease stepsize and try again */
			NextStep *= sqrt(0.9*rel_tol*QHEAT_TOL/RelError);

			/* stepsize too small, this can happen at extreme low temperatures */
			if( NextStep < 5.*DBL_EPSILON*u1[k] )
			{
				*ErrorCode = MAX2(*ErrorCode,QH_STEP_FAIL);
				break;
			}

			continue;
		}
		else
		{
			/* step was OK, adjust stepsize */
			k += 2;

			/* >>chng 03 nov 10, safeguard against division by zero, PvH */
			NextStep *= MIN2(pow(0.9*rel_tol*QHEAT_TOL/MAX2(RelError,1.e-50),1./3.),4.);
			NextStep = MIN2(NextStep,Lambda[k]/Phi[0]);
		}

		dPdlnT[k-1] = p[k-1]*qtemp[k-1]*uderiv(qtemp[k-1],nd);
		dPdlnT[k] = p[k]*qtemp[k]*uderiv(qtemp[k],nd);

		lgAllNegSlope = lgAllNegSlope && ( dPdlnT[k] < dPdlnT[k-2] );

		if( dPdlnT[k-1] > maxVal )
		{
			maxVal = dPdlnT[k-1];
			nmax = k-1;
		}
		if( dPdlnT[k] > maxVal )
		{
			maxVal = dPdlnT[k];
			nmax = k;
		}

		RadCooling += dCool;

/* 		if( nzone >= 870 ) { */
/*  		printf(" k %ld T[k] %.6e U[k] %.6e p[k] %.6e dPdlnT[k] %.6e\n",k-1,qtemp[k-1],u1[k-1],p[k-1],dPdlnT[k-1]); */
/*  		printf(" k %ld T[k] %.6e U[k] %.6e p[k] %.6e dPdlnT[k] %.6e\n",k,qtemp[k],u1[k],p[k],dPdlnT[k]); */
/* 		} */

		/* if qtmin is far too low, p[k] can become extremely large, exceeding
		 * even double precision range. the following check should prevent overflows */
		/* >>chng 01 nov 07, sqrt(DBL_MAX) -> sqrt(DBL_MAX/100.) so that sqrt(p[k]*p[k+1]) is safe */
		if( p[k] > sqrt(DBL_MAX/100.) ) 
		{
			*ErrorCode = MAX2(*ErrorCode,QH_LOW_FAIL);
			break;

#if 0
			/* >>chng 01 apr 21, change DBL_MAX -> DBL_MAX/100. to work around
			 * a bug in the Compaq C compiler V6.3-025 when invoked with -fast, PvH */
			for( j=0; j <= k; j++ ) 
			{
				p[j] /= DBL_MAX/100.;
				dPdlnT[j] /= DBL_MAX/100.;
			}
			RadCooling /= DBL_MAX/100.;
#endif
		}

		/* this may catch a bug in the Compaq C compiler V6.3-025
		 * if this gets triggered, try compiling with -ieee */
		ASSERT( p[k] > 0. && dPdlnT[k] > 0. && RadCooling > 0. );

		/* do a check for negative slope and if there will be enough room to store results */
		/* >>chng 02 jan 30, do this test only every NQTEST steps, PvH */
		/* if( k >= MIN2(NQTEST,NQGRID-2) ) */
		if( k > 0 && k%NQTEST == 0 )
		{
			double wid, avStep, factor;
			/* >>chng 02 jul 31, do additional test for HIGH_FAIL,
			 * first test before main loop doesn't catch everything, PvH */
			if( lgAllNegSlope )
			{
				*ErrorCode = MAX2(*ErrorCode,QH_HIGH_FAIL);
				break;
			}

			/* this is a lower limit for the total width of the probability distr */
			/* >>chng 02 jan 30, corrected calculation of wid and avStep, PvH */
			wid = (sqrt(-log(PROB_CUTOFF_LO)/(4.*LN_TWO)) +
			       sqrt(-log(PROB_CUTOFF_HI)/(4.*LN_TWO)))*fwhm/Umax;
			avStep = log(u1[k]/u1[0])/(double)k;
			/* make an estimate for the number of steps needed */
			/* >>chng 02 jan 30, included factor 1.5 because stepsize increases near peak, PvH */
			/* >>chng 02 aug 06, changed 1.5 to sliding scale because of laqheat2.in test, PvH */
			factor = 1.1 + 3.9*(1.0 - sqrt((double)k/(double)NQGRID));
			if( wid/avStep > factor*(double)NQGRID )
			{
				*ErrorCode = MAX2(*ErrorCode,QH_ARRAY_FAIL);
				break;
			}
		}

		/* if we run out of room to store results, do regular break
		 * the code below will sort out if integration is valid or not */
		if( k >= NQGRID-2 )
		{
			*ErrorCode = MAX2(*ErrorCode,QH_ARRAY_FAIL);
			break;
		}

		/* force thermal equilibrium of the grains */
		fac = RadCooling*gv.bin[nd]->cnv_GR_pCM3*EN1RYD/gv.bin[nd]->GrainHeat;

		/* this is regular stop criterion */
		if( dPdlnT[k] < dPdlnT[k-1] && dPdlnT[k]/fac < PROB_CUTOFF_HI )
		{
			break;
		}
	}

	if( l == MAX_LOOP )
		*ErrorCode = MAX2(*ErrorCode,QH_LOOP_FAIL);

	nok = k;
	nbin = k+1;

	/* there are insufficient bins to attempt rebinning */
	if( *ErrorCode < QH_NO_REBIN && nbin < NQMIN ) 
		*ErrorCode = MAX2(*ErrorCode,QH_NBIN_FAIL);

	/* >>chng 02 aug 07, do some basic checks on the distribution first */
	if( *ErrorCode < QH_NO_REBIN )
		ScanProbDistr(u1,dPdlnT,nbin,maxVal,nmax,&nstart,&nstart2,&nend,nWideFail,ErrorCode);

	if( *ErrorCode >= QH_NO_REBIN )
	{
#		ifdef DEBUG_FUN
		fputs( " <->GetProbDistr_LowLimit()\n", debug_fp );
#		endif
		return;
	}

	for( j=0; j < nbin; j++ )
	{
		p[j] /= fac;
		dPdlnT[j] /= fac;
	}
	RadCooling /= fac;

	/* >>chng 02 aug 08, moved RebinQHeatResults from here further down, this improves new_tmin estimate */
	*new_tmin = 0.;
	for( j=nstart; j < nbin; j++ ) 
	{
		if( dPdlnT[j] < PROB_CUTOFF_LO ) 
		{
			*new_tmin = qtemp[j];
		}
		else 
		{
			if( j == nstart ) 
			{
				if( dPdlnT[j] > SAFETY*PROB_CUTOFF_LO )
					*ErrorCode = MAX2(*ErrorCode,QH_BOUND_FAIL);

				/* >>chng 02 aug 11, use nstart2 for more reliable extrapolation */
				if( dPdlnT[nstart2] < 0.999*dPdlnT[nstart2+NSTARTUP] ) 
				{
					/* >>chng 02 aug 09, new formula for extrapolating new_tmin, PvH */
					/* this assumes that at low temperatures the behaviour
					 * is as follows:   dPdlnT(T) = C1 * exp( -C2/T**3 ) */
					double T1 = qtemp[nstart2];
					double T2 = qtemp[nstart2+NSTARTUP];
					double delta_y = log(dPdlnT[nstart2+NSTARTUP]/dPdlnT[nstart2]);
					double c2 = delta_y/(1./POW3(T1)-1./POW3(T2));
					double help = c2/POW3(T1) + log(dPdlnT[nstart2]/PROB_CUTOFF_LO);
					*new_tmin = pow(c2/help,1./3.);
				}

				/* >>chng 04 nov 09, in case of lower bound failure, assure qtmin is lowered, PvH */ 
				if( dPdlnT[j] > SAFETY*PROB_CUTOFF_LO && *new_tmin >= qtmin1 )
				{
					double delta_x = log(qtemp[nstart2+NSTARTUP]/qtemp[nstart2]);
					double delta_y = log(dPdlnT[nstart2+NSTARTUP]/dPdlnT[nstart2]);
					delta_x *= log(PROB_CUTOFF_LO/dPdlnT[nstart2])/delta_y;
					*new_tmin = qtemp[nstart2]*exp(delta_x);
					if( *new_tmin < qtmin1 )
						/* in general this estimate will be too low -> use geometric mean */
						*new_tmin = sqrt( *new_tmin*qtmin1 );
					else
						/* last resort... */
						*new_tmin = qtmin1/DEF_FAC;
				}
			}
			break;
		}
	}
	*new_tmin = MAX3(*new_tmin,qtmin1/DEF_FAC,GRAIN_TMIN);

	ASSERT( *new_tmin < gv.bin[nd]->tedust );

	/* >>chng 02 jan 30, prevent excessive looping when prob distribution simply won't fit, PvH */
	if( dPdlnT[nbin-1] > SAFETY*PROB_CUTOFF_HI )
	{
		if( *ErrorCode == QH_ARRAY_FAIL || *ErrorCode == QH_LOOP_FAIL )
		{
			++(*nWideFail);

			if( *nWideFail < WIDE_FAIL_MAX )
			{
				/* this indicates that low end was OK, but we ran out of room
				 * to store the high end -> try GetProbDistr_HighLimit instead */
				*ErrorCode = MAX2(*ErrorCode,QH_DELTA_FAIL);
			}
			else
			{
				*ErrorCode = MAX2(*ErrorCode,QH_WIDE_FAIL);
			}
		}
		else
		{
			*ErrorCode = MAX2(*ErrorCode,QH_BOUND_FAIL);
		}
	}

	/* >>chng 01 may 11, rebin the quantum heating results
	 *
	 * for grains in intense radiation fields, the code needs high resolution for stability
	 * and therefore produces lots of small bins, even though the grains never make large
	 * excurions from the equilibrium temperature; adding in the resulting spectra in RT_diffuse
	 * takes up an excessive amount of CPU time where most CPU is spent on grains for which
	 * the quantum treatment matters least, and moreover on temperature bins with very low
	 * probability; rebinning the results on a coarser grid should help reduce the overhead */
	/* >>chng 02 aug 07, use nstart and nend while rebinning */

	nbin = RebinQHeatResults(nd,nstart,nend,p,qtemp,qprob,dPdlnT,u1,delu,Lambda,ErrorCode);

	/* >>chng 01 jul 13, add fail-safe for failure in RebinQHeatResults */
	if( nbin == 0 )
	{
#		ifdef DEBUG_FUN
		fputs( " <->GetProbDistr_LowLimit()\n", debug_fp );
#		endif
		return;
	}

	*qnbin = nbin;

	sum = 0.;
	for( j=0; j < nbin; j++ )
	{
		sum += qprob[j];
	}

	/* the fact that the probability normalization fails probably indicates
	 * that the distribution is too close to a delta function to resolve */
	if( fabs(sum-1.) > PROB_TOL )
		*ErrorCode = MAX2(*ErrorCode,QH_CONV_FAIL);
		
	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ,
			 "    zone %ld %s nbin %ld nok %ld nbad %ld total prob %.4e rel_tol %.3e new_tmin %.4e\n",
			 nzone,gv.bin[nd]->chDstLab,nbin,nok,nbad,sum,rel_tol,*new_tmin );
	}

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


/* try two consecutive integration steps using stepsize "step/2." (yielding p[k]),
 * and also one double integration step using stepsize "step" (yielding p2k).
 * the difference fabs(p2k-p[k])/(3.*p[k]) can be used to estimate the relative
 * accuracy of p[k] and will be used to adapt the stepsize to an optimal value */
STATIC double TryDoubleStep(double u1[],
			    double delu[],
			    double p[],
			    double qtemp[],
			    double Lambda[],
			    double Phi[],
			    double PhiDrv[],
			    double step,
			    /*@out@*/ double *cooling,
			    long index,
			    long nd,
			    /*@out@*/ int *lgBoundFail)
{
	long i,
	  j,
	  jlo,
	  k=-1000;
	double bval_jk,
	  cooling2,
	  p2k,
	  p_max,
	  RelErrCool,
	  RelErrPk,
	  sum,
	  sum2 = -DBL_MAX,
	  trap1,
	  trap12 = -DBL_MAX,
	  trap2,
	  uhi,
	  ulo,
	  umin,
	  y;

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

	/* sanity checks */
	ASSERT( index >= 0 && index < NQGRID-2 && nd >= 0 && nd < gv.nBin && step > 0. );

	ulo = rfield.anu[0];
	/* >>chng 01 nov 29, rfield.nflux -> gv.qnflux, PvH */
	/* >>chng 03 jan 26, gv.qnflux -> gv.bin[nd]->qnflux, PvH */
	uhi = rfield.anu[gv.bin[nd]->qnflux-1];

	/* >>chng 01 nov 21, skip initial bins if they have very low probability */
	p_max = 0.;
	for( i=0; i <= index; i++ )
		p_max = MAX2(p_max,p[i]);
	jlo = 0;
	while( p[jlo] < PROB_CUTOFF_LO*p_max )
		jlo++;

	for( i=1; i <= 2; i++ )
	{
		int lgErr;
		long ipLo = 0;
		long ipHi = gv.bin[nd]->qnflux-1;
		k = index + i;
		delu[k] = step/2.;
		u1[k] = u1[k-1] + delu[k];
		qtemp[k] = inv_ufunct(u1[k],nd,lgBoundFail);
		splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,log(qtemp[k]),&y,&lgErr);
		*lgBoundFail = *lgBoundFail || lgErr;
		Lambda[k] = exp(y)*gv.bin[nd]->cnv_H_pGR/EN1RYD;

		sum = sum2 = 0.;
		trap1 = trap2 = trap12 = 0.;

		/* this loop uses up a large fraction of the total CPU time, it should be well optimized */
		for( j=jlo; j < k; j++ )
		{
			umin = u1[k] - u1[j];

			if( umin >= uhi ) 
			{
				/* for increasing j, umin will decrease monotonically. If ( umin > uhi ) holds for
				 * the current value of j, it must have held for the previous value as well. Hence
				 * both trap1 and trap2 are zero at this point and we would only be adding zero
				 * to sum. Therefore we skip this step to save CPU time */
				continue;
			}
			else if( umin > ulo )
			{
				/* do a bisection search such that rfield.anu[ipLo] <= umin < rfield.anu[ipHi]
				 * ipoint doesn't always give the correct index into anu (it works for AnuOrg)
				 * bisection search is also faster, which is important here to save CPU time.
				 * on the first iteration ipLo equals 0 and the first while loop will be skipped;
				 * after that umin is monotonically decreasing, and ipHi is retained from the
				 * previous iteration since it is a valid upper limit; ipLo will equal ipHi-1 */
				long ipStep = 1;
				/* >>chng 03 feb 03 rjrw: hunt for lower bracket */
				while( rfield.anu[ipLo] > (float)umin )
				{
					ipHi = ipLo;
					ipLo -= ipStep;
					if( ipLo <= 0 ) 
					{
						ipLo = 0;
						break;
					}
					ipStep *= 2;
				}
				/* now do regular bisection search */
				while( ipHi-ipLo > 1 )
				{
					long ipMd = (ipLo+ipHi)/2;
					if( rfield.anu[ipMd] > (float)umin )
						ipHi = ipMd;
					else
						ipLo = ipMd;
				}
				bval_jk = Phi[ipLo] + (umin - rfield.anu[ipLo])*PhiDrv[ipLo];
			}
			else
			{
				bval_jk = Phi[0];
			}

			/* these two quantities are needed to take double step from index -> index+2 */
			trap12 = trap1;
			sum2 = sum;

			/* bval_jk*gv.bin[nd]->cnv_CM3_pGR is the total excitation rate from j to k and
			 * higher due to photon absorptions and particle collisions, it already implements
			 * Eq. 2.17 of Guhathakurtha & Draine, in events/grain/s */
			/* >>chng 00 mar 27, factored in hden (in Phi), by PvH */
			/* >>chng 00 mar 29, add in contribution due to particle collisions, by PvH */
			/* >>chng 01 mar 30, moved multiplication of bval_jk with gv.bin[nd]->cnv_CM3_pGR
			 *   out of loop, PvH */
			trap2 = p[j]*bval_jk;
			/* Trapezoidal rule, multiplication with factor 0.5 is done outside loop */
			sum += (trap1 + trap2)*delu[j];
			trap1 = trap2;
		}

		/* >>chng 00 mar 27, multiplied with delu, by PvH */
		/* >>chng 00 apr 05, taken out Lambda[0], improves convergence at low end dramatically!, by PvH */
		/* qprob[k] = sum*gv.bin[nd]->cnv_CM3_pGR*delu[k]/(Lambda[k] - Lambda[0]); */
		/* this expression includes factor 0.5 from trapezoidal rule above */
		/* p[k] = 0.5*(sum + (trap1 + p[k]*Phi[0])*delu[k])/Lambda[k] */
		p[k] = (sum + trap1*delu[k])/(2.*Lambda[k] - Phi[0]*delu[k]);

		ASSERT( p[k] > 0. );
	}

	/* this is estimate for p[k] using one double step of size "step" */
	p2k = (sum2 + trap12*step)/(2.*Lambda[k] - Phi[0]*step);

	ASSERT( p2k > 0. );

	RelErrPk = fabs(p2k-p[k])/p[k];

	/* this is radiative cooling due to the two probability bins we just added */
	/* simple trapezoidal rule will not do here, RelErrCool would never converge */
	*cooling = log_integral(u1[k-2],p[k-2]*Lambda[k-2],u1[k-1],p[k-1]*Lambda[k-1]);
	*cooling += log_integral(u1[k-1],p[k-1]*Lambda[k-1],u1[k],p[k]*Lambda[k]);

	/* same as cooling, but now for double step of size "step" */
	cooling2 = log_integral(u1[k-2],p[k-2]*Lambda[k-2],u1[k],p2k*Lambda[k]);

	/* p[0] is not reliable, so ignore convergence test on cooling on first step */
	RelErrCool = ( index > 0 ) ? fabs(cooling2-(*cooling))/(*cooling) : 0.;

/*  	printf( " TryDoubleStep p[k-1] %.4e p[k] %.4e p2k %.4e\n",p[k-1],p[k],p2k ); */

#	ifdef DEBUG_FUN
	fputs( " <->TryDoubleStep()\n", debug_fp );
#	endif
	/* error scales as O(step^3), so this is relative accuracy of p[k] or cooling */
	return MAX2(RelErrPk,RelErrCool)/3.;
}


/* calculate logarithmic integral from (xx1,yy1) to (xx2,yy2) */
STATIC double log_integral(double xx1,
			   double yy1,
			   double xx2,
			   double yy2)
{
	double eps,
	  result,
	  xx;

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

	/* sanity checks */
	ASSERT( xx1 > 0. && yy1 > 0. && xx2 > 0. && yy2 > 0. );

	xx = log(xx2/xx1);
	eps = log((xx2/xx1)*(yy2/yy1));
	if( fabs(eps) < 1.e-4 )
	{
		result = xx1*yy1*xx*((eps/6. + 0.5)*eps + 1.);
	}
	else
	{
		result = (xx2*yy2 - xx1*yy1)*xx/eps;
	}

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


/* scan the probability distribution for valid range */
STATIC void ScanProbDistr(double u1[],      /* u1[nbin] */
			  double dPdlnT[],  /* dPdlnT[nbin] */
			  long nbin,
			  double maxVal,
			  long nmax,
			  /*@out@*/long *nstart,
			  /*@out@*/long *nstart2,
			  /*@out@*/long *nend,
			  long *nWideFail,
			  QH_Code *ErrorCode)
{
	int lgSetLo,
	  lgSetHi;
	long i;
	double deriv_max,
	  minVal;
	const double MIN_FAC_LO = 1.e4;
	const double MIN_FAC_HI = 1.e4;

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

	/* sanity checks */
	ASSERT( nbin > 0 && nmax >= 0 && nmax < nbin && maxVal > 0. );

	/* sometimes the probability distribution will start falling before settling on
	 * a rising slope. In such a case nstart will point to the first rising point,
	 * while nstart2 will point to the point with the steepest derivative on the
	 * rising slope. The code will use the distribution from nstart to nend as a
	 * valid probability distribution, but the will use the points near nstart2
	 * to extrapolate a new value for qtmin if needed */
	minVal = maxVal;
	*nstart = nmax;
	for( i=nmax; i >= 0; --i )
	{
		if( dPdlnT[i] < minVal )
		{
			*nstart = i;
			minVal = dPdlnT[i];
		}
	}
	deriv_max = 0.;
	*nstart2 = nmax;
	for( i=nmax; i > *nstart; --i )
	{
		double deriv = log(dPdlnT[i]/dPdlnT[i-1])/log(u1[i]/u1[i-1]);
		if( deriv > deriv_max )
		{
			*nstart2 = i-1;
			deriv_max = deriv;
		}
	}
	*nend = nbin-1;

	/* now do quality control; these checks are more stringent than the ones in GetProbDistr_LowLimit */
	lgSetLo = ( nmax >= *nend || maxVal/dPdlnT[*nend] < MIN_FAC_HI );
	/* >>chng 03 jan 22, prevent problems if both dPdlnT and its derivative are continuously rising,
	 *                   in which case both lgSetLo and lgSetHi are set and QH_WIDE_FAIL is triggered;
	 *                   this can happen if qtmin is far too low; triggered by pahtest.in, PvH */
	if( lgSetLo )
		/* use relaxed test if lgSetLo is already set */
		lgSetHi = ( nmax <= *nstart || maxVal/dPdlnT[*nstart] < MIN_FAC_LO );
	else
		lgSetHi = ( nmax <= *nstart2 || maxVal/dPdlnT[*nstart2] < MIN_FAC_LO );

	if( lgSetLo && lgSetHi )
	{
		++(*nWideFail);

		if( *nWideFail >= WIDE_FAIL_MAX )
			*ErrorCode = MAX2(*ErrorCode,QH_WIDE_FAIL);
	}

	if( lgSetLo )
		*ErrorCode = MAX2(*ErrorCode,QH_LOW_FAIL);

	if( lgSetHi )
		*ErrorCode = MAX2(*ErrorCode,QH_HIGH_FAIL);

	/* there are insufficient bins to attempt rebinning */
	if( *ErrorCode < QH_NO_REBIN && (*nend - *nstart) < NQMIN ) 
		*ErrorCode = MAX2(*ErrorCode,QH_NBIN_FAIL);

	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "    ScanProbDistr nstart %ld nstart2 %ld nend %ld nmax %ld maxVal %.3e",
			 *nstart,*nstart2,*nend,nmax,maxVal );
		fprintf( ioQQQ, " dPdlnT[nstart] %.3e dPdlnT[nstart2] %.3e dPdlnT[nend] %.3e code %d\n",
			 dPdlnT[*nstart],dPdlnT[*nstart2],dPdlnT[*nend],*ErrorCode );
	}

	if( *ErrorCode >= QH_NO_REBIN )
	{
		*nstart = -1;
		*nstart2 = -1;
		*nend = -2;
	}

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


/* rebin the quantum heating results to speed up RT_diffuse */
STATIC long RebinQHeatResults(long nd,
			      long nstart,
			      long nend,
			      double p[],
			      double qtemp[],
			      double qprob[],
			      double dPdlnT[],
			      double u1[],
			      double delu[],
			      double Lambda[],
			      QH_Code *ErrorCode)
{
	long i,
	  newnbin;
	double fac,
	  help,
	  mul_fac,
	  *new_delu/*[NQGRID]*/,
	  *new_dPdlnT/*[NQGRID]*/,
	  *new_Lambda/*[NQGRID]*/,
	  *new_p/*[NQGRID]*/,
	  *new_qprob/*[NQGRID]*/,
	  *new_qtemp/*[NQGRID]*/,
	  *new_u1/*[NQGRID]*/,
	  PP1,
	  PP2,
	  RadCooling,
	  T1,
	  T2,
	  Ucheck,
	  uu1,
	  uu2;

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

	/* sanity checks */
	ASSERT( nd >= 0 && nd < gv.nBin );
	/* >>chng 02 aug 07, changed oldnbin -> nstart..nend */
	ASSERT( nstart >= 0 && nstart < nend && nend < NQGRID );

	/* leading entries may have become very small or zero -> skip */
	for( i=nstart; i <= nend && dPdlnT[i] < PROB_CUTOFF_LO; i++ ) {}

	/* >>chng 04 oct 17, add fail-safe to keep lint happy, but this should never happen... */
	if( i >= NQGRID )
	{
		*ErrorCode = MAX2(*ErrorCode,QH_REBIN_FAIL);

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

	/* >>chng 02 aug 15, use malloc to prevent stack overflows */
	if( ( new_delu = (double*)MALLOC((size_t)(NQGRID*sizeof(double))) ) == NULL )
		BadMalloc();
	if( ( new_dPdlnT = (double*)MALLOC((size_t)(NQGRID*sizeof(double))) ) == NULL )
		BadMalloc();
	if( ( new_Lambda = (double*)MALLOC((size_t)(NQGRID*sizeof(double))) ) == NULL )
		BadMalloc();
	if( ( new_p = (double*)MALLOC((size_t)(NQGRID*sizeof(double))) ) == NULL )
		BadMalloc();
	if( ( new_qprob = (double*)MALLOC((size_t)(NQGRID*sizeof(double))) ) == NULL )
		BadMalloc();
	if( ( new_qtemp = (double*)MALLOC((size_t)(NQGRID*sizeof(double))) ) == NULL )
		BadMalloc();
	if( ( new_u1 = (double*)MALLOC((size_t)(NQGRID*sizeof(double))) ) == NULL )
		BadMalloc();

	newnbin = 0;

	T1 = qtemp[i];
	PP1 = p[i];
	uu1 = u1[i];

	/* >>chng 04 feb 01, change 2.*NQMIN -> 1.5*NQMIN, PvH */
	help = pow(qtemp[nend]/qtemp[i],1./(1.5*NQMIN));
	mul_fac = MIN2(QT_RATIO,help);

	Ucheck = u1[i];
	RadCooling = 0.;

	while( i < nend )
	{
		int lgBoundErr;
		int lgDone= FALSE;
		double s0 = 0.;
		double s1 = 0.;
		double wid = 0.;
		double xx,y;

		T2 = T1*mul_fac;

		do 
		{
			double p1,p2,L1,L2,frac,slope;
			if( qtemp[i] <= T1 && T1 <= qtemp[i+1] )
			{
				/* >>chng 01 nov 15, copy uu2 into uu1 instead, PvH */
				/* uu1 = ufunct(T1,nd); */
				frac = log(T1/qtemp[i]);
				slope = log(p[i+1]/p[i])/log(qtemp[i+1]/qtemp[i]);
				p1 = p[i]*exp(frac*slope);
				slope = log(Lambda[i+1]/Lambda[i])/log(qtemp[i+1]/qtemp[i]);
				L1 = Lambda[i]*exp(frac*slope);
			}
			else
			{
				/* >>chng 01 nov 15, copy uu2 into uu1 instead, PvH */
				/* uu1 = u1[i]; */
				p1 = p[i];
				L1 = Lambda[i];
			}
			if( qtemp[i] <= T2 && T2 <= qtemp[i+1] )
			{
				/* >>chng 02 apr 30, make sure this doesn't point beyond valid range, PvH */
				help = ufunct(T2,nd,&lgBoundErr);
				uu2 = MIN2(help,u1[i+1]);
				ASSERT( !lgBoundErr ); /* this should never be triggered */
				frac = log(T2/qtemp[i]);
				slope = log(p[i+1]/p[i])/log(qtemp[i+1]/qtemp[i]);
				p2 = p[i]*exp(frac*slope);
				slope = log(Lambda[i+1]/Lambda[i])/log(qtemp[i+1]/qtemp[i]);
				L2 = Lambda[i]*exp(frac*slope);
				lgDone = TRUE;
			}
			else
			{
				uu2 = u1[i+1];
				p2 = p[i+1];
				L2 = Lambda[i+1];
				/* >>chng 01 nov 15, this caps the range in p(U) integrated in one bin
				 * it helps avoid spurious QH_BOUND_FAIL's when flank is very steep, PvH */
				if( MAX2(p2,PP1)/MIN2(p2,PP1) > sqrt(SAFETY) )
				{
					lgDone = TRUE;
					T2 = qtemp[i+1];
				}
				++i;
			}
			PP2 = p2;
			wid += uu2 - uu1;
			/* sanity check */
			ASSERT( wid >= 0. );
			s0 += log_integral(uu1,p1,uu2,p2);
			s1 += log_integral(uu1,p1*L1,uu2,p2*L2);
			uu1 = uu2;

		} while( i < nend && ! lgDone );

		/* >>chng 01 nov 14, if T2 == qtemp[oldnbin-1], the code will try another iteration
		 * break here to avoid zero divide, the assert on Ucheck tests if we are really finished */
		/* >>chng 01 dec 04, change ( s0 == 0. ) to ( s0 <= 0. ), PvH */
		if( s0 <= 0. )
		{
			ASSERT( wid == 0. );
			break;
		}

		new_qprob[newnbin] = s0;
		new_Lambda[newnbin] = s1/s0;
		xx = log(new_Lambda[newnbin]*EN1RYD*gv.bin[nd]->cnv_GR_pH);
		splint_safe(gv.bin[nd]->dstems,gv.dsttmp,gv.bin[nd]->dstslp,NDEMS,xx,&y,&lgBoundErr);
		ASSERT( !lgBoundErr ); /* this should never be triggered */
		new_qtemp[newnbin] = exp(y);
		new_u1[newnbin] = ufunct(new_qtemp[newnbin],nd,&lgBoundErr);
		ASSERT( !lgBoundErr ); /* this should never be triggered */
		new_delu[newnbin] = wid;
		new_p[newnbin] = new_qprob[newnbin]/new_delu[newnbin];
		new_dPdlnT[newnbin] = new_p[newnbin]*new_qtemp[newnbin]*uderiv(new_qtemp[newnbin],nd);

		Ucheck += wid;
		RadCooling += new_qprob[newnbin]*new_Lambda[newnbin];

		T1 = T2;
		PP1 = PP2;
		++newnbin;
	}

	/* >>chng 01 jul 13, add fail-safe */
	if( newnbin < NQMIN )
	{
		*ErrorCode = MAX2(*ErrorCode,QH_REBIN_FAIL);

		free(new_delu);
		free(new_dPdlnT);
		free(new_Lambda);
		free(new_p);
		free(new_qprob);
		free(new_qtemp);
		free(new_u1);

#		ifdef DEBUG_FUN
		fputs( " <->RebinQHeatResults()\n", debug_fp );
#		endif
		return 0;
	}
		
	fac = RadCooling*EN1RYD*gv.bin[nd]->cnv_GR_pCM3/gv.bin[nd]->GrainHeat;

	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "     RebinQHeatResults found tol1 %.4e tol2 %.4e\n",
			 fabs(u1[nend]/Ucheck-1.), fabs(fac-1.) );
	}

	/* do quality control */
	/* >>chng 02 apr 30, tighten up check, PvH */
	ASSERT( fabs(u1[nend]/Ucheck-1.) < 10.*sqrt((double)(nend-nstart+newnbin))*DBL_EPSILON );

	if( fabs(fac-1.) > CONSERV_TOL )
		*ErrorCode = MAX2(*ErrorCode,QH_CONV_FAIL);

	for( i=0; i < newnbin; i++ )
	{
		/* renormalize the distribution to assure energy conservation */
		p[i] = new_p[i]/fac;
		qtemp[i] = new_qtemp[i];
		qprob[i] = new_qprob[i]/fac;
		dPdlnT[i] = new_dPdlnT[i]/fac;
		u1[i] = new_u1[i];
		delu[i] = new_delu[i];
		Lambda[i] = new_Lambda[i];

		/* sanity checks */
		ASSERT( qtemp[i] > 0. && qprob[i] > 0. );

/*  		printf(" rk %ld T[k] %.6e U[k] %.6e p[k] %.6e dPdlnT[k] %.6e\n",i,qtemp[i],u1[i],p[i],dPdlnT[i]); */
	}

	free(new_delu);
	free(new_dPdlnT);
	free(new_Lambda);
	free(new_p);
	free(new_qprob);
	free(new_qtemp);
	free(new_u1);

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


/* calculate approximate probability distribution in high intensity limit */
STATIC void GetProbDistr_HighLimit(long nd,
				   double TolFac,
				   double Umax,
				   double fwhm,
				   /*@out@*/double qtemp[],
				   /*@out@*/double qprob[],
				   /*@out@*/double dPdlnT[],
				   /*@out@*/double *tol,
				   /*@out@*/long *qnbin,
				   /*@out@*/double *new_tmin,
				   /*@out@*/QH_Code *ErrorCode)
{
	int lgBoundErr,
	  lgErr;
	long i,
	  nbin;
	double c1,
	  c2,
	  delu[NQGRID],
	  fac,
	  fac1,
	  fac2,
	  help1,
	  help2,
	  L1,
	  L2,
	  Lambda[NQGRID],
	  mul_fac,
	  p[NQGRID],
	  p1,
	  p2,
	  RadCooling,
	  sum,
	  T1,
	  T2,
	  Tlo,
	  Thi,
	  Ulo,
	  Uhi,
	  uu1,
	  uu2,
	  xx,
	  y;

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

	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "   GetProbDistr_HighLimit called for grain %s\n", gv.bin[nd]->chDstLab );
	}

	c1 = sqrt(4.*LN_TWO/PI)/fwhm*exp(-POW2(fwhm/Umax)/(16.*LN_TWO));
	c2 = 4.*LN_TWO*POW2(Umax/fwhm);

	fac1 = fwhm/Umax*sqrt(-log(PROB_CUTOFF_LO)/(4.*LN_TWO));
	/* >>chng 03 nov 10, safeguard against underflow, PvH */
	help1 = Umax*exp(-fac1);
	help2 = exp(gv.bin[nd]->DustEnth[0]);
	Ulo = MAX2(help1,help2);
	/* >>chng 03 jan 28, ignore lgBoundErr on lower boundary, low-T grains have negigible emission, PvH */
	Tlo = inv_ufunct(Ulo,nd,&lgBoundErr);

	fac2 = fwhm/Umax*sqrt(-log(PROB_CUTOFF_HI)/(4.*LN_TWO));
	/* >>chng 03 nov 10, safeguard against overflow, PvH */
	if( fac2 > log(DBL_MAX/10.) )
	{
		*ErrorCode = MAX2(*ErrorCode,QH_WIDE_FAIL);

#		ifdef DEBUG_FUN
		fputs( " <->GetProbDistr_HighLimit()\n", debug_fp );
#		endif
		return;
	}
	Uhi = Umax*exp(fac2);
	Thi = inv_ufunct(Uhi,nd,&lgBoundErr);

	nbin = 0;

	T1 = Tlo;
	uu1 = ufunct(T1,nd,&lgErr);
	lgBoundErr = lgBoundErr || lgErr;
	help1 = log(uu1/Umax);
	p1 = c1*exp(-c2*POW2(help1));
	splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,log(T1),&y,&lgErr);
	lgBoundErr = lgBoundErr || lgErr;
	L1 = exp(y)*gv.bin[nd]->cnv_H_pGR/EN1RYD;

	/* >>chng 03 nov 10, safeguard against underflow, PvH */
	if( uu1*p1*L1 < 1.e5*DBL_MIN )
	{
		*ErrorCode = MAX2(*ErrorCode,QH_WIDE_FAIL);

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

	/* >>chng 04 feb 01, change 2.*NQMIN -> 1.2*NQMIN, PvH */
	help1 = pow(Thi/Tlo,1./(1.2*NQMIN));
	mul_fac = MIN2(QT_RATIO,help1);

	sum = 0.;
	RadCooling = 0.;

	do 
	{
		double s0,s1,wid;

		T2 = T1*mul_fac;
		uu2 = ufunct(T2,nd,&lgErr);
		lgBoundErr = lgBoundErr || lgErr;
		help1 = log(uu2/Umax);
		p2 = c1*exp(-c2*POW2(help1));
		splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,log(T2),&y,&lgErr);
		lgBoundErr = lgBoundErr || lgErr;
		L2 = exp(y)*gv.bin[nd]->cnv_H_pGR/EN1RYD;

		wid = uu2 - uu1;
		s0 = log_integral(uu1,p1,uu2,p2);
		s1 = log_integral(uu1,p1*L1,uu2,p2*L2);

		qprob[nbin] = s0;
		Lambda[nbin] = s1/s0;
		xx = log(Lambda[nbin]*EN1RYD*gv.bin[nd]->cnv_GR_pH);
		splint_safe(gv.bin[nd]->dstems,gv.dsttmp,gv.bin[nd]->dstslp,NDEMS,xx,&y,&lgErr);
		lgBoundErr = lgBoundErr || lgErr;
		qtemp[nbin] = exp(y);
		delu[nbin] = wid;
		p[nbin] = qprob[nbin]/delu[nbin];
		dPdlnT[nbin] = p[nbin]*qtemp[nbin]*uderiv(qtemp[nbin],nd);

		sum += qprob[nbin];
		RadCooling += qprob[nbin]*Lambda[nbin];

		T1 = T2;
		uu1 = uu2;
		p1 = p2;
		L1 = L2;

		++nbin;

	} while( T2 < Thi && nbin < NQGRID );

	fac = RadCooling*EN1RYD*gv.bin[nd]->cnv_GR_pCM3/gv.bin[nd]->GrainHeat;

	for( i=0; i < nbin; ++i )
	{
		qprob[i] /= fac;
		dPdlnT[i] /= fac;
	}

	*tol = sum/fac;
	*qnbin = nbin;
	*new_tmin = qtemp[0];
	*ErrorCode = MAX2(*ErrorCode,QH_ANALYTIC);

	/* do quality control */
	if( TolFac > STRICT )
		*ErrorCode = MAX2(*ErrorCode,QH_ANALYTIC_RELAX);

	if( lgBoundErr )
		*ErrorCode = MAX2(*ErrorCode,QH_THIGH_FAIL);

	if( fabs(sum/fac-1.) > PROB_TOL )
		*ErrorCode = MAX2(*ErrorCode,QH_CONV_FAIL);

	if( dPdlnT[0] > SAFETY*PROB_CUTOFF_LO || dPdlnT[nbin-1] > SAFETY*PROB_CUTOFF_HI )
		*ErrorCode = MAX2(*ErrorCode,QH_BOUND_FAIL);
		
	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "     GetProbDistr_HighLimit found tol1 %.4e tol2 %.4e\n",
			 fabs(sum-1.), fabs(sum/fac-1.) );
		fprintf( ioQQQ, "    zone %ld %s nbin %ld total prob %.4e new_tmin %.4e\n",
			 nzone,gv.bin[nd]->chDstLab,nbin,sum/fac,*new_tmin );
	}

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


/* calculate derivative of the enthalpy function dU/dT (aka the specific heat) at a given temperature, in Ryd/K */
STATIC double uderiv(double temp, 
		     long int nd)
{
	enth_type ecase;
	long int i,
	  j;
	double N_C,
	  N_H;
	double deriv = 0.,
	  hok[3] = {1275., 1670., 4359.},
	  numer,
	  dnumer,
	  denom,
	  ddenom,
	  x;


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

	if( temp <= 0. ) 
	{
		fprintf( ioQQQ, " uderiv called with non-positive temperature: %.6e\n" , temp );
		puts( "[Stop in uderiv]" );
		cdEXIT(EXIT_FAILURE);
	}
	ASSERT( nd >= 0 && nd < gv.nBin );

	ecase = gv.which_enth[gv.bin[nd]->matType];
	switch( ecase )
	{
	case ENTH_CAR:
		numer = (4.15e-22/EN1RYD)*pow(temp,3.3);
		dnumer = (3.3*4.15e-22/EN1RYD)*pow(temp,2.3);
		denom = 1. + 6.51e-03*temp + 1.5e-06*temp*temp + 8.3e-07*pow(temp,2.3);
		ddenom = 6.51e-03 + 2.*1.5e-06*temp + 2.3*8.3e-07*pow(temp,1.3);
		/* dU/dT for pah/graphitic grains in Ryd/K, derived from:
		 * >>refer	grain	physics	Guhathakurta & Draine, 1989, ApJ, 345, 230 */
		deriv = (dnumer*denom - numer*ddenom)/POW2(denom);
		break;
	case ENTH_CAR2:
		/* dU/dT for graphite grains in Ryd/K, using eq 9 of */
		/* >>refer	grain	physics	Draine B.T., and Li A., 2001, ApJ, 551, 807 */
		deriv = (DebyeDeriv(temp/863.,2) + 2.*DebyeDeriv(temp/2504.,2))*BOLTZMANN/EN1RYD;
		break;
	case ENTH_SIL:
		/* dU/dT for silicate grains (and grey grains) in Ryd/K */
		/* limits of tlim set above, 0 and DBL_MAX, so always OK */
		/* >>refer	grain	physics	Guhathakurta & Draine, 1989, ApJ, 345, 230 */
		for( j = 0; j < 4; j++ )
		{
			if( temp > tlim[j] && temp <= tlim[j+1] )
			{
				deriv = cval[j]*pow(temp,power[j]);
				break;
			}
		}
		break;
	case ENTH_SIL2:
		/* dU/dT for silicate grains in Ryd/K, using eq 11 of */
		/* >>refer	grain	physics	Draine B.T., and Li A., 2001, ApJ, 551, 807 */
		deriv = (2.*DebyeDeriv(temp/500.,2) + DebyeDeriv(temp/1500.,3))*BOLTZMANN/EN1RYD;
		break;
	case ENTH_PAH:
		/* dU/dT for PAH grains in Ryd/K, using eq A.4 of */
		/* >>refer	grain	physics	Dwek E., Arendt R.G., Fixsen D.J. et al., 1997, ApJ, 475, 565 */
		/* this expression is only valid upto 2000K */
		x = log10(MIN2(temp,2000.));
		deriv = pow(10.,-21.26+3.1688*x-0.401894*POW2(x))/EN1RYD;
		break;
	case ENTH_PAH2:
		/* dU/dT for PAH grains in Ryd/K, approximately using eq 33 of */
		/* >>refer	grain	physics	Draine B.T., and Li A., 2001, ApJ, 551, 807 */
		/* N_C and N_H should actually be NINT(N_C) and NINT(N_H),
		 * but this can lead to FP overflow for very large grains */
		N_C = NO_ATOMS(nd);
		if( N_C <= 25. )
		{
			N_H = 0.5*N_C;
		}
		else if( N_C <= 100. )
		{
			N_H = 2.5*sqrt(N_C);
		}
		else
		{
			N_H = 0.25*N_C;
		}
		deriv = 0.;
		for( i=0; i < 3; i++ )
		{
			double help1 = hok[i]/temp;
			if( help1 < 300. )
			{
				double help2 = exp(help1);
				double help3 = ( help1 < 1.e-7 ) ? help1*(1.+0.5*help1) : help2-1.;
				deriv += N_H/(N_C-2.)*POW2(help1)*help2/POW2(help3)*BOLTZMANN/EN1RYD;
			}
		}
		deriv += (DebyeDeriv(temp/863.,2) + 2.*DebyeDeriv(temp/2504.,2))*BOLTZMANN/EN1RYD;
		break;
	default:
		fprintf( ioQQQ, " uderiv called with unknown type for enthalpy function: %d\n", ecase );
		puts( "[Stop in uderiv]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* >>chng 00 mar 23, use formula 3.1 of Guhathakurtha & Draine, by PvH */
	/* >>chng 03 jan 17, use MAX2(..,1) to prevent crash for extremely small grains, PvH */
	deriv *= MAX2(NO_ATOMS(nd)-2.,1.);

	if( deriv <= 0. ) 
	{
		fprintf( ioQQQ, " uderiv finds non-positive derivative: %.6e, what's up?\n" , deriv );
		puts( "[Stop in uderiv]" );
		cdEXIT(EXIT_FAILURE);
	}

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


/* calculate the enthalpy of a grain at a given temperature, in Ryd */
STATIC double ufunct(double temp, 
		     long int nd,
		     /*@out@*/ int *lgBoundErr)
{
	double enthalpy,
	  y;

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

	if( temp <= 0. ) 
	{
		fprintf( ioQQQ, " ufunct called with non-positive temperature: %.6e\n" , temp );
		puts( "[Stop in ufunct]" );
		cdEXIT(1);
	}
	ASSERT( nd >= 0 && nd < gv.nBin );

	/* >>chng 02 apr 22, interpolate in DustEnth array to get enthalpy, by PvH */
	splint_safe(gv.dsttmp,gv.bin[nd]->DustEnth,gv.bin[nd]->EnthSlp,NDEMS,log(temp),&y,lgBoundErr);
	enthalpy = exp(y);

	ASSERT( enthalpy > 0. );

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


/* this is the inverse of ufunct: determine grain temperature as a function of enthalpy */
STATIC double inv_ufunct(double enthalpy, 
			 long int nd,
			 /*@out@*/ int *lgBoundErr)
{
	double temp,
	  y;

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

	if( enthalpy <= 0. ) 
	{
		fprintf( ioQQQ, " inv_ufunct called with non-positive enthalpy: %.6e\n" , enthalpy );
		puts( "[Stop in ufunct]" );
		cdEXIT(EXIT_FAILURE);
	}
	ASSERT( nd >= 0 && nd < gv.nBin );

	/* >>chng 02 apr 22, interpolate in DustEnth array to get temperature, by PvH */
	splint_safe(gv.bin[nd]->DustEnth,gv.dsttmp,gv.bin[nd]->EnthSlp2,NDEMS,log(enthalpy),&y,lgBoundErr);
	temp = exp(y);

	ASSERT( temp > 0. );

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


/* initialize interpolation arrys for grain enthalpy */
void InitEnthalpy()
{
	double C_V1,
	  C_V2,
	  tdust1,
	  tdust2;
	long int i,
	  j,
	  nd;

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

	for( nd=0; nd < gv.nBin; nd++ )
	{
		tdust2 = GRAIN_TMIN;
		C_V2 = uderiv(tdust2,nd);
		/* at low temps, C_V = C*T^3 -> U = C*T^4/4 = C_V*T/4 */
		gv.bin[nd]->DustEnth[0] = C_V2*tdust2/4.;
		tdust1 = tdust2;
		C_V1 = C_V2;

		for( i=1; i < NDEMS; i++ )
		{
			double tmid,Cmid;
			tdust2 = exp(gv.dsttmp[i]);
			C_V2 = uderiv(tdust2,nd);
			tmid = sqrt(tdust1*tdust2);
			/* this ensures accuracy for silicate enthalpy */
			for( j=1; j < 4; j++ )
			{
				if( tdust1 < tlim[j] && tlim[j] < tdust2 )
				{
					tmid = tlim[j];
					break;
				}
			}
			Cmid = uderiv(tmid,nd);
			gv.bin[nd]->DustEnth[i] = gv.bin[nd]->DustEnth[i-1] +
				log_integral(tdust1,C_V1,tmid,Cmid) +
				log_integral(tmid,Cmid,tdust2,C_V2);
			tdust1 = tdust2;
			C_V1 = C_V2;
		}
	}

	/* conversion for logarithmic interpolation */
	for( nd=0; nd < gv.nBin; nd++ )
	{
		for( i=0; i < NDEMS; i++ )
		{
			gv.bin[nd]->DustEnth[i] = log(gv.bin[nd]->DustEnth[i]);
		}
	}

	/* now find slopes from spline fit */
	for( nd=0; nd < gv.nBin; nd++ )
	{
		/* set up coefficients for splint */
		spline(gv.dsttmp,gv.bin[nd]->DustEnth,NDEMS,2e31,2e31,gv.bin[nd]->EnthSlp);
		spline(gv.bin[nd]->DustEnth,gv.dsttmp,NDEMS,2e31,2e31,gv.bin[nd]->EnthSlp2);
	}

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


/* helper function for calculating specific heat, uses Debye theory */
STATIC double DebyeDeriv(double x,
			 long n)
{
	long i,
	  nn;
	double res,
	  *xx,
	  *rr,
	  *aa,
	  *ww;

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

	ASSERT( x > 0. );
	ASSERT( n == 2 || n == 3 );

	if( x < 0.001 )
	{
		/* for general n this is Gamma(n+2)*zeta(n+1)*powi(x,n) */
		if( n == 2 )
		{
			res = 6.*1.202056903159594*POW2(x);
		}
		else if( n == 3 )
		{
			res = 24.*1.082323233711138*POW3(x);
		}
		else
			/* added to keep lint happy - note that assert above confimred that n is 2 or 3,
			 * but lint flagged possible flow without defn of res */
			TotalInsanity();
	}
	else
	{
		nn = 4*MAX2(4,2*(long)(0.05/x));
		xx = (double *)MALLOC(sizeof(double)*(unsigned)nn);
		rr = (double *)MALLOC(sizeof(double)*(unsigned)nn);
		aa = (double *)MALLOC(sizeof(double)*(unsigned)nn);
		ww = (double *)MALLOC(sizeof(double)*(unsigned)nn);
		gauss_legendre(nn,xx,aa);
		gauss_init(nn,0.,1.,xx,aa,rr,ww);

		res = 0.;
		for( i=0; i < nn; i++ )
		{
			double help1 = rr[i]/x;
			if( help1 < 300. )
			{
				double help2 = exp(help1);
				double help3 = ( help1 < 1.e-7 ) ? help1*(1.+0.5*help1) : help2-1.;
				res += ww[i]*powi(rr[i],n+1)*help2/POW2(help3);
			}
		}
		res /= POW2(x);

		free(ww);
		free(aa);
		free(rr);
		free(xx);
	}

#	ifdef DEBUG_FUN
	fputs( " <->DebyeDeriv()\n", debug_fp );
#	endif
	return (double)n*res;
}
/*lint +e817 Conceivably negative subscript*/
/*lint +e661 possible access of out-of-bounds pointer */
