/* 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 */
/*grain main routine to converge grains thermal solution */
#include "cddefines.h"
#include "physconst.h"
#include "rfield.h"
#include "hmi.h"
#include "trace.h"
#include "converge.h"
#include "ionbal.h"
#include "thermal.h"
#include "phycon.h"
#include "doppvel.h"
#include "taulines.h"
#include "co.h"
#include "heavy.h"
#include "negdrg.h"
#include "interpolate.h"
#include "dense.h"
#include "ipoint.h"
#include "grainvar.h"
#include "grains.h"
#include "elementnames.h"

#ifdef STATIC
#undef STATIC
#endif

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

/* option to revert to old grain potential algorithm; uncomment to activate */
/* #define OLD_POT_ALGO 1 */

/* the next three defines are for debugging purposes only, uncomment to activate */
/*  #define WD_TEST2 1 */
/*  #define IGNORE_GRAIN_ION_COLLISIONS 1 */
/*  #define IGNORE_THERMIONIC 1 */

#define NINT(X) ((long)((X) < 0. ? (X)-0.5 : (X)+0.5))

/* no parentheses around PTR needed since it needs to be an lvalue */
#define FREE_CHECK(PTR) { ASSERT( PTR != NULL ); free( PTR ); PTR = NULL; }
#define FREE_SAFE(PTR) { if( PTR != NULL ) free( PTR ); PTR = NULL; }

#define POT2CHRG(X) ((X)*EVRYD/ELEM_CHARGE*gv.bin[nd]->Capacity - 1.)
#define CHRG2POT(X) (((X)+1.)*ELEM_CHARGE/EVRYD/gv.bin[nd]->Capacity)
#define ONE_ELEC   (ELEM_CHARGE/EVRYD/gv.bin[nd]->Capacity)

/*#define EPSM(Z_0,Z) ( ((Z) <= (Z_0)-1) ? 1. : (((Z) >= (Z_0)) ? 0. : (double)((Z_0)-(Z))) )*/
#define EPSP(Z_0,Z) ( ((Z) <= (Z_0)) ? 0. : (((Z) >= (Z_0)+1) ? 1. : (double)((Z)-(Z_0))) )

static const int INCL_TUNNEL = TRUE;
static const int NO_TUNNEL = FALSE;

static const int ALL_STAGES = TRUE;
/* static const int NONZERO_STAGES = FALSE; */

/* counts how many times GrainDrive has been called,
 * set to zero in GrainZero */
static long int nCalledGrainDrive;

/*================================================================================*/
/* these are used for setting up grain emissivities in InitEmissivities() */

/* NTOP is number of bins for temps between GRAIN_TMID and GRAIN_TMAX */
static const long NTOP = NDEMS/5;

/*================================================================================*/
/* miscellaneous grain physics */

/* constant for thermionic emissions, 7.501e20 e/cm^2/s/K^2 */
static const double THERMCONST = PI4*ELECTRON_MASS*POW2(BOLTZMANN)/POW3(HPLANCK);

/* sticking probabilities */
static const double STICK_ELEC = 0.5;
static const double STICK_ION = 1.0;

/* mean free path for electron penetrating grain, in cm */
static const double MEAN_PATH = 1.e-7;

/*================================================================================*/
/* these are used when iterating the grain charge in GrainCharge() */
static const double TOLER = CONSERV_TOL/10.;
#ifdef OLD_POT_ALGO
static const double BIG_POTENTIAL = 5.;

/* >>chng 04 jan 24, lowered BRACKET_STEP from 0.2 -> 0.001, tests show that in 85%
 *      of the cases this bracket is large enough to contain the new pot value, PvH */
static const double BRACKET_STEP = 0.001;
#endif
static const long BRACKET_MAX = 50L;

/* this is the largest number of charge bins that may be in use at any one time; the
 * remaining NCHS-NCHU charge bins are used for backing up data for possible later use */
#define NCHU (NCHS/3)

/* maximum number of tries to converge charge/temperature in GrainChargeTemp() */
static const long CT_LOOP_MAX = 10L;

/* maximum number of tries to converge grain temperature in GrainChargeTemp() */
static const long T_LOOP_MAX = 30L;

/* these will become the new tolerance levels used throughout the code */
static double HEAT_TOLER = DBL_MAX;
static double HEAT_TOLER_BIN = DBL_MAX;
static double CHRG_TOLER = DBL_MAX;
static double CHRG_TOLER_BIN = DBL_MAX;

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

static int lgGvInitialized = FALSE;

/* a_0 thru a_2 constants for calculating IP_V and EA, in cm */
static const double AC0 = 3.e-9;
static const double AC1G = 4.e-8;
static const double AC2G = 7.e-8;

/* >>chng 01 sep 13, create structure for dynamically allocating backup store, PvH */
/* the following are placeholders for intermediate results that depend on grain type,
 * however they are only used inside the main loop over nd in GrainChargeTemp(),
 * so it is OK to reuse the same memory for each grain bin separately. */
/* >>chng 01 dec 19, added entry for Auger rate, PvH */
/* >>chng 04 jan 17, created substructure for individual charge states, PvH */
/* >>chng 04 jan 20, moved cache data into gv data structure, PvH */

/* initialize grain emissivities */
STATIC void InitEmissivities(void);
/* PlanckIntegral compute total radiative cooling due to large grains */
STATIC double PlanckIntegral(double,long,long);
/* invalidate charge dependent data from previous iteration */
STATIC void NewChargeData(long);
/* iterate grain charge and temperature */
STATIC void GrainChargeTemp(void);
#ifdef OLD_POT_ALGO
/* GrainCharge compute grains charge */
STATIC void GrainCharge(long,/*@out@*/double*,double);
/* bracket the solution for the grain potential */
STATIC void BracketGrainPot(long,/*@in@*/double*,/*@out@*/double*,/*@out@*/double*,/*@out@*/double*,
			    /*@in@*/double*,/*@out@*/double*,/*@out@*/double*,/*@out@*/double*,
			    /*@in@*/double*,/*@out@*/double*,/*@out@*/double*,/*@out@*/double*);
/* GrainElecRecomb compute electron recombination onto grain surface */
STATIC double GrainElecRecomb(long);
/* GrainElecEmis computes electron loss rate for the grain species
 * due to photoelectric effect, positive ion recombinations on the
 * grain surface and thermionic emissions */
STATIC double GrainElecEmis(long,/*@out@*/double*);
#else
/* GrainCharge compute grains charge */
STATIC void GrainCharge(long,/*@out@*/double*);
#endif
/* grain electron recombination rates for single charge state */
STATIC double GrainElecRecomb1(long,long,/*@out@*/double*,/*@out@*/double*);
/* grain electron emission rates for single charge state */
STATIC double GrainElecEmis1(long,long,/*@out@*/double*,/*@out@*/double*,
			     /*@out@*/double*,/*@out@*/double*,/*@out@*/double*);
/* correction factors for grain charge screening (including image potential
 * to correct for polarization of the grain as charged particle approaches). */
STATIC void GrainScreen(long,long,long,double*,double*);
/* helper function for GrainScreen */
STATIC double ThetaNu(double);
#ifdef OLD_POT_ALGO
/* update items that depend on grain potential */
STATIC void UpdatePot(long,long);
/* calculate charge state populations */
STATIC void GetFracPop(long,double,long,/*@out@*/long*);
#else
/* update items that depend on grain potential */
STATIC void UpdatePot(long,long,long,/*@out@*/double[],/*@out@*/double[]);
/* calculate charge state populations */
STATIC void GetFracPop(long,long,/*@in@*/double[],/*@in@*/double[],/*@out@*/long*);
#endif
/* this routine updates all quantities that depend only on grain charge and radius */
STATIC void UpdatePot1(long,long,long,long);
/* this routine updates all quantities that depend on grain charge, radius and temperature */
STATIC void UpdatePot2(long,long);
/* find highest ionization stage with non-zero population */
STATIC long HighestIonStage(void);
/* determine charge Z0 ion recombines to upon impact on grain */
STATIC void UpdateRecomZ0(long,long,int);
/* helper routine for UpdatePot */
STATIC void GetPotValues(long,long,/*@out@*/double*,/*@out@*/double*,/*@out@*/double*,
			 /*@out@*/double*,/*@out@*/double*,int);
/* given grain nd in charge state nz, and incoming ion (nelem,ion),
 * detemine outgoing ion (nelem,Z0) and chemical energy ChEn released
 * ChemEn is net contribution of ion recombination to grain heating */
STATIC void GrainIonColl(long,long,long,long,const double[],const double[],/*@out@*/long*,
			 /*@out@*/float*,/*@out@*/float*);
/* initialize ion recombination rates on grain species nd */
STATIC void GrainChrgTransferRates(long);
/* this routine updates all grain quantities that depend on radius, except gv.dstab and gv.dstsc */
STATIC void GrainUpdateRadius1(void);
/* this routine adds all the grain opacities in gv.dstab and gv.dstsc */
STATIC void GrainUpdateRadius2(int);
/* GrainTemperature computes grains temperature, and gas cooling */
STATIC void GrainTemperature(long,/*@out@*/float*,/*@out@*/double*,/*@out@*/double*,
			     /*@out@*/double*);
/* GrainCollHeating computes grains collisional heating cooling */
STATIC void GrainCollHeating(long,/*@out@*/float*,/*@out@*/float*,/*@out@*/float[]);
#if 0
/* GrainCollHeating compute grains collisional heating cooling dur to ions/electrons */
STATIC void GrainCollHeating1(long,long,/*@out@*/double*,/*@out@*/double*,/*@out@*/double*,
			      /*@out@*/double*,/*@out@*/double*);
/* GrainCollHeating2 compute grains collisional heating/cooling due to molecules */
STATIC void GrainCollHeating2(long,/*@out@*/double*,/*@out@*/double*,/*@out@*/double*);
#endif
/*GrnVryDpth set grains abundance as a function of depth into cloud*/
STATIC double GrnVryDpth(long);

/*lint -e662 creation of out of bounds pointer */

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

/* this routine is called by zero(), so it should contain initializations
 * that need to be done every time before the input lines are parsed */
void GrainZero(void)
{
	long int nelem, ion, ion_to;

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

	gv.lgAnyDustVary = FALSE;
	gv.TotalEden = 0.;
	gv.dHeatdT = 0.;
	gv.lgQHeatAll = FALSE;
	/* gv.lgGrainElectrons - should grain electron source/sink be included in overall electron sum?
	 * default is true, set false with no grain electrons command */
	gv.lgGrainElectrons = TRUE;
	gv.lgQHeatOn = TRUE;
	gv.lgDHetOn = TRUE;
	gv.lgDColOn = TRUE;
	gv.GrainMetal = 1.;
	gv.lgBakes = FALSE;
	gv.lgLeidenBakesPAH_heat = FALSE;
	gv.nChrgRequested = NCHRG_DEFAULT;
	gv.ReadPtr = 0L;
	/* by default grains always reevaluated - command grains reevaluate off sets to false */
	gv.lgReevaluate = TRUE;
	/* set true with "set grain leiden hack", use simple grain extinction */
	gv.lgLeidenGrnOpacHack = FALSE;

	/* counts how many times GrainDrive has been called */
	nCalledGrainDrive = 0;

	/* this is option to turn off all grain physics while leaving
	 * the opacity in, set false with no grain physics command */
	gv.lgGrainPhysicsOn = TRUE;

	/* the following entries define the physical behavior of each type of grains
	 * (entropy function, expression for Zmin and ionization potential, etc) */
	/* >>chng 02 sep 18, defined which_ial for all material types, needed for special rfi files, PvH */
	gv.which_enth[MAT_CAR] = ENTH_CAR;
	gv.which_zmin[MAT_CAR] = ZMIN_CAR;
	gv.which_pot[MAT_CAR] = POT_CAR;
	gv.which_ial[MAT_CAR] = IAL_CAR;
	gv.which_pe[MAT_CAR] = PE_CAR;
	gv.which_strg[MAT_CAR] = STRG_CAR;
	gv.which_H2distr[MAT_CAR] = H2_CAR;

	gv.which_enth[MAT_SIL] = ENTH_SIL;
	gv.which_zmin[MAT_SIL] = ZMIN_SIL;
	gv.which_pot[MAT_SIL] = POT_SIL;
	gv.which_ial[MAT_SIL] = IAL_SIL;
	gv.which_pe[MAT_SIL] = PE_SIL;
	gv.which_strg[MAT_SIL] = STRG_SIL;
	gv.which_H2distr[MAT_SIL] = H2_SIL;

	gv.which_enth[MAT_PAH] = ENTH_PAH;
	gv.which_zmin[MAT_PAH] = ZMIN_CAR;
	gv.which_pot[MAT_PAH] = POT_CAR;
	gv.which_ial[MAT_PAH] = IAL_CAR;
	gv.which_pe[MAT_PAH] = PE_CAR;
	gv.which_strg[MAT_PAH] = STRG_CAR;
	gv.which_H2distr[MAT_PAH] = H2_CAR;

	gv.which_enth[MAT_CAR2] = ENTH_CAR2;
	gv.which_zmin[MAT_CAR2] = ZMIN_CAR;
	gv.which_pot[MAT_CAR2] = POT_CAR;
	gv.which_ial[MAT_CAR2] = IAL_CAR;
	gv.which_pe[MAT_CAR2] = PE_CAR;
	gv.which_strg[MAT_CAR2] = STRG_CAR;
	gv.which_H2distr[MAT_CAR2] = H2_CAR;

	gv.which_enth[MAT_SIL2] = ENTH_SIL2;
	gv.which_zmin[MAT_SIL2] = ZMIN_SIL;
	gv.which_pot[MAT_SIL2] = POT_SIL;
	gv.which_ial[MAT_SIL2] = IAL_SIL;
	gv.which_pe[MAT_SIL2] = PE_SIL;
	gv.which_strg[MAT_SIL2] = STRG_SIL;
	gv.which_H2distr[MAT_SIL2] = H2_SIL;

	gv.which_enth[MAT_PAH2] = ENTH_PAH2;
	gv.which_zmin[MAT_PAH2] = ZMIN_CAR;
	gv.which_pot[MAT_PAH2] = POT_CAR;
	gv.which_ial[MAT_PAH2] = IAL_CAR;
	gv.which_pe[MAT_PAH2] = PE_CAR;
	gv.which_strg[MAT_PAH2] = STRG_CAR;
	gv.which_H2distr[MAT_PAH2] = H2_CAR;

	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		for( ion=0; ion <= nelem+1; ion++ )
		{
			for( ion_to=0; ion_to <= nelem+1; ion_to++ )
			{
				gv.GrainChTrRate[nelem][ion][ion_to] = 0.f;
			}
		}
	}

	/* >>>chng 01 may 08, return memory possibly allocated in previous calls to cloudy(), PvH
	 * this routine MUST be called before ParseCommands() so that grain commands find a clean slate */
	ReturnGrainBins();

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


/* this routine is called by IterStart(), so anything that needs to be reset before each
 * iteration starts should be put here; typically variables that are integrated over radius */
void GrainStartIter(void)
{
	long nd;

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

	if( gv.lgDustOn && gv.lgGrainPhysicsOn )
	{
		for( nd=0; nd < gv.nBin; nd++ )
		{
			/* >>chng 97 jul 5, save and reset this
			 * save grain potential */
			gv.bin[nd]->dstpotsav = gv.bin[nd]->dstpot;
			gv.bin[nd]->qtmin = ( gv.bin[nd]->qtmin_zone1 > 0. ) ?
				gv.bin[nd]->qtmin_zone1 : DBL_MAX;
			gv.bin[nd]->avdust = 0.;
			gv.bin[nd]->avdpot = 0.;
			gv.bin[nd]->avdft = 0.;
			gv.bin[nd]->avDGRatio = 0.;
			gv.bin[nd]->TeGrainMax = -1.f;
			gv.bin[nd]->qpres = QHEAT_INIT_RES;
			gv.bin[nd]->lgEverQHeat = FALSE;
			gv.bin[nd]->QHeatFailures = 0L;
			gv.bin[nd]->lgQHTooWide = FALSE;
			gv.bin[nd]->lgPAHsInIonizedRegion = FALSE;
			gv.bin[nd]->nChrgOrg = gv.bin[nd]->nChrg;
		}
	}

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


/* this routine is called by IterRestart(), so anything that needs to be
 * reset or saved after an iteration is finished should be put here */
void GrainRestartIter(void)
{
	long nd;

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

	if( gv.lgDustOn && gv.lgGrainPhysicsOn )
	{
		for( nd=0; nd < gv.nBin; nd++ )
		{
			/* >>chng 97 jul 5, reset grain potential
			 * reset grain to pential to initial value from previous iteration */
			gv.bin[nd]->dstpot = gv.bin[nd]->dstpotsav;
			gv.bin[nd]->nChrg = gv.bin[nd]->nChrgOrg;
		}
	}

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


/* this routine is called by ParseSet() */
void SetNChrgStates(long nChrg)
{
#	ifdef DEBUG_FUN
	fputs( "<+>SetNChrgStates()\n", debug_fp );
#	endif

	ASSERT( nChrg >= 2 && nChrg <= NCHU );
	gv.nChrgRequested = nChrg;

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


long NewGrainBin(void)
{
	long nd, nz;

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

	ASSERT( lgGvInitialized );

	if( gv.nBin >= NDUST ) 
	{
		fprintf( ioQQQ, " The code has run out of grain bins; increase NDUST and recompile.\n" );
		puts( "[Stop in NewGrainBin]" );
		cdEXIT(EXIT_FAILURE);
	}
	nd = gv.nBin;

	ASSERT( gv.bin[nd] == NULL ); /* prevent memory leaks */
	if( ( gv.bin[nd] = (GrainBin *)MALLOC(sizeof(GrainBin)) ) == NULL ) 
		BadMalloc();

	/* the first 4 are allocated in mie_read_opc, the rest in InitGrains */
	gv.bin[nd]->dstab1 = NULL;
	gv.bin[nd]->pure_sc1 = NULL;
	gv.bin[nd]->asym = NULL;
	gv.bin[nd]->inv_att_len = NULL;
	gv.bin[nd]->y1 = NULL;
	for( nz=0; nz < NCHS; nz++ )
		gv.bin[nd]->chrg[nz] = NULL;

	gv.lgDustOn = TRUE;
	gv.bin[nd]->lgQHeat = FALSE;
	gv.bin[nd]->qnflux = LONG_MAX;
	gv.bin[nd]->lgDustVary = FALSE;
	gv.bin[nd]->DustDftVel = FLT_MAX;
	gv.bin[nd]->TeGrainMax = FLT_MAX;
	/* NB - this number should not be larger than NCHU */
	gv.bin[nd]->nChrg = gv.nChrgRequested;
	/* this must be zero for the first solutions to be able to converge */
	/* >>chng 00 jun 19, tedust has to be greater than zero
	 * to prevent division by zero in GrainElecEmis and GrainCollHeating, PvH */
	gv.bin[nd]->tedust = 1.f;
	gv.bin[nd]->GrainHeat = DBL_MAX/10.;
	gv.bin[nd]->GrainGasCool = DBL_MAX/10.;
	/* used to check that energy scale in grains opacity files is same as
	 * current cloudy scale */
	gv.bin[nd]->EnergyCheck = 0.;
	gv.bin[nd]->dstAbund = -FLT_MAX;
	gv.bin[nd]->dstfactor = 1.f;
	gv.bin[nd]->cnv_H_pGR = -DBL_MAX;
	gv.bin[nd]->cnv_H_pCM3 = -DBL_MAX;
	gv.bin[nd]->cnv_CM3_pGR = -DBL_MAX;
	gv.bin[nd]->cnv_CM3_pH = -DBL_MAX;
	gv.bin[nd]->cnv_GR_pH = -DBL_MAX;
	gv.bin[nd]->cnv_GR_pCM3 = -DBL_MAX;
	/* >>chng 04 feb 05, zero this rate in case "no molecules" is set, will.in, PvH */
	gv.bin[nd]->rate_h2_form_grains_used = 0.;

	gv.nBin++;

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


void ReturnGrainBins(void)
{
	long nd, nz;

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

	if( lgGvInitialized )
	{
		/* >>chng 01 sep 12, allocate/free [rfield.nupper] arrays dynamically */
		for( nd=0; nd < gv.nBin; nd++ ) 
		{
			ASSERT( gv.bin[nd] != NULL );

			FREE_SAFE( gv.bin[nd]->dstab1 );
			FREE_SAFE( gv.bin[nd]->pure_sc1 );
			FREE_SAFE( gv.bin[nd]->asym );
			FREE_SAFE( gv.bin[nd]->y1 );
			FREE_SAFE( gv.bin[nd]->inv_att_len );

			for( nz=0; nz < NCHS; nz++ )
			{
				if( gv.bin[nd]->chrg[nz] != NULL )
				{
					FREE_SAFE( gv.bin[nd]->chrg[nz]->yhat );
					FREE_SAFE( gv.bin[nd]->chrg[nz]->cs_pdt );
					FREE_SAFE( gv.bin[nd]->chrg[nz]->fac1 );
					FREE_SAFE( gv.bin[nd]->chrg[nz]->fac2 );

					FREE_CHECK( gv.bin[nd]->chrg[nz] );
				}
			}

			FREE_CHECK( gv.bin[nd] );
		}

		FREE_SAFE( gv.anumin );
		FREE_SAFE( gv.anumax );
		FREE_SAFE( gv.dstab );
		FREE_SAFE( gv.dstsc );
		FREE_SAFE( gv.GrainEmission );
		FREE_SAFE( gv.GraphiteEmission );
		FREE_SAFE( gv.SilicateEmission );
	}
	else
	{
		/* >>chng 01 sep 12, moved initialization of data from NewGrainBin to here, PvH */
		/* >>chng 01 may 08, make sure bin pointers are properly initialized, PvH */
		for( nd=0; nd < NDUST; nd++ )
		{
			gv.bin[nd] = NULL;
		}

		gv.anumin = NULL;
		gv.anumax = NULL;
		gv.dstab = NULL;
		gv.dstsc = NULL;
		gv.GrainEmission = NULL;
		gv.GraphiteEmission = NULL;
		gv.SilicateEmission = NULL;

		lgGvInitialized = TRUE;
	}

	gv.lgDustOn = FALSE;
	gv.nBin = 0;

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


/*InitGrains, called one time by opacitycreateall at initialization of calculation, 
 * called after commands have been parsed,
 * not after every iteration or every model */
void InitGrains(void)
{
	long int i,
	  nelem,
	  nd,
	  nz;

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

	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, " InitGrains called.\n" );
	}


	/* >>chng 01 sep 12, allocate/free [rfield.nupper] arrays dynamically */
	ASSERT( gv.anumin == NULL ); /* prevent memory leaks */
	if( ( gv.anumin = (float*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(float))) ) == NULL ) 
		BadMalloc();
	ASSERT( gv.anumax == NULL );
	if( ( gv.anumax = (float*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(float))) ) == NULL ) 
		BadMalloc();
	ASSERT( gv.dstab == NULL );
	if( ( gv.dstab = (double*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(double))) ) == NULL ) 
		BadMalloc();
	ASSERT( gv.dstsc == NULL );
	if( ( gv.dstsc = (double*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(double))) ) == NULL ) 
		BadMalloc();
	ASSERT( gv.GrainEmission == NULL );
	if( ( gv.GrainEmission = (float*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(float))) ) == NULL ) 
		BadMalloc();
	ASSERT( gv.GraphiteEmission == NULL );
	if( ( gv.GraphiteEmission = (float*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(float))) ) == NULL ) 
		BadMalloc();
	ASSERT( gv.SilicateEmission == NULL );
	if( ( gv.SilicateEmission = (float*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(float))) ) == NULL ) 
		BadMalloc();
	
	/* sanity check */
	ASSERT( gv.nBin >= 0 && gv.nBin < NDUST );
	for( nd=gv.nBin; nd < NDUST; nd++ ) 
	{
		ASSERT( gv.bin[nd] == NULL );
	}

	/* >>chng 02 jan 15, initialize to zero in case grains are not used, needed in IonIron(), PvH */
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		gv.elmSumAbund[nelem] = 0.f;
	}

	for( i=0; i < rfield.nupper; i++ )
	{
		gv.dstab[i] = 0.;
		gv.dstsc[i] = 0.;
		/* >>chng 01 sep 12, moved next three initializations from GrainZero(), PvH */
		gv.GrainEmission[i] = 0.;
		gv.SilicateEmission[i] = 0.;
		gv.GraphiteEmission[i] = 0.;
	}

	if( !gv.lgDustOn )
	{
		/* grains are not on, set all heating/cooling agents to zero */
		gv.GrainHeatInc = 0.;
		gv.GrainHeatDif = 0.;
		gv.GrainHeatLya = 0.;
		gv.GrainHeatCollSum = 0.;
		gv.GrainHeatSum = 0.;
		gv.GasCoolColl = 0.;
		thermal.heating[0][13] = 0.;
		thermal.heating[0][14] = 0.;
		thermal.heating[0][25] = 0.;

		if( trace.lgTrace && trace.lgDustBug )
		{
			fprintf( ioQQQ, " InitGrains exits.\n" );
		}

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

	HEAT_TOLER = conv.HeatCoolRelErrorAllowed / 10.;
	HEAT_TOLER_BIN = HEAT_TOLER / sqrt(gv.nBin);
	CHRG_TOLER = conv.EdenErrorAllowed / 10.;
	CHRG_TOLER_BIN = CHRG_TOLER / sqrt(gv.nBin);

	for( nd=0; nd < gv.nBin; nd++ )
	{
		double help,atoms,p_rad,ThresInf,dum[4];
		long low1,low2,low3,lowm;

		/* sanity check */
		ASSERT( gv.bin[nd] != NULL );

		/* this is QHEAT ALL command */
		if( gv.lgQHeatAll )
		{
			gv.bin[nd]->lgQHeat = TRUE;
		}

		/* this is NO GRAIN QHEAT command, always takes precedence */
		if( !gv.lgQHeatOn ) 
		{
			gv.bin[nd]->lgQHeat = FALSE;
		}

		/* >>chng 04 jun 01, disable quantum heating when constant grain temperature is used, PvH */
		if( thermal.ConstGrainTemp > 0. )
		{
			gv.bin[nd]->lgQHeat = FALSE;
		}

#ifndef IGNORE_QUANTUM_HEATING
		gv.bin[nd]->lgQHTooWide = FALSE;
		gv.bin[nd]->qtmin = DBL_MAX;
#endif

		/* this is total grain depletion factor, GrainMetal is depletion
		 * set with metals command, normally one */
		gv.bin[nd]->dstfactor *= gv.GrainMetal;

		/* >>chng 01 nov 21, grain abundance may depend on radius,
		 * invalidate for now; GrainUpdateRadius1() will set correct value */
		gv.bin[nd]->dstAbund = -FLT_MAX;

		gv.bin[nd]->qtmin_zone1 = -1.;

		if( gv.bin[nd]->DustWorkFcn < rfield.anu[0] || gv.bin[nd]->DustWorkFcn > rfield.anu[rfield.nupper] )
		{
			fprintf( ioQQQ, " Grain work function for %s has insane value: %.4e\n",
				 gv.bin[nd]->chDstLab,gv.bin[nd]->DustWorkFcn );
			puts( "[Stop in InitGrains]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* >>chng 01 sep 12, allocate/free [rfield.nupper] arrays dynamically */
		ASSERT( gv.bin[nd]->y1 == NULL ); /* prevent memory leaks */
		if( (gv.bin[nd]->y1 = (float*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(float)))) == NULL )
			BadMalloc();
		for( nz=0; nz < NCHS; nz++ )
		{
			ASSERT( gv.bin[nd]->chrg[nz] == NULL );
			if( ( gv.bin[nd]->chrg[nz] = (ChargeBin *)MALLOC(sizeof(ChargeBin)) ) == NULL ) 
				BadMalloc();
			gv.bin[nd]->chrg[nz]->yhat = NULL;
			gv.bin[nd]->chrg[nz]->cs_pdt = NULL;
			gv.bin[nd]->chrg[nz]->fac1 = NULL;
			gv.bin[nd]->chrg[nz]->fac2 = NULL;
			gv.bin[nd]->chrg[nz]->nfill = 0;

			gv.bin[nd]->chrg[nz]->DustZ = LONG_MIN;
			gv.bin[nd]->chrg[nz]->FracPop = -DBL_MAX;
			gv.bin[nd]->chrg[nz]->tedust = 1.f;
		}

		for( i=0; i < rfield.nupper; i++ )
		{
			double alpha,beta,af,bf;

			beta = gv.bin[nd]->AvRadius*gv.bin[nd]->inv_att_len[i];
			if( beta > 1.e-4 ) 
			{
				bf = POW2(beta) - 2.*beta + 2. - 2.*exp(-beta);
			}
			else 
			{
				bf = POW3(beta)/3.;
			}
			alpha = beta + gv.bin[nd]->AvRadius/MEAN_PATH;
			if( alpha > 1.e-4 ) 
			{
				af = POW2(alpha) - 2.*alpha + 2. - 2.*exp(-alpha);
			}
			else 
			{
				af = POW3(alpha)/3.;
			}

			/* this is size-dependent geometrical yield enhancement
			 * defined in Weingartner & Draine, 2000 */
			gv.bin[nd]->y1[i] = (float)(POW2(beta/alpha)*af/bf);
		}

		ASSERT( gv.bin[nd]->nChrg >= 2 && gv.bin[nd]->nChrg <= NCHU );

		/* >>chng 00 jun 19, this value is absolute lower limit for the grain
		 * potential, electrons cannot be bound for lower values..., PvH */
		if( gv.lgBakes )
		{
			/* this corresponds to
			 * >>refer	grain	physics	Bakes & Tielens, 1994, ApJ, 427, 822 */
			help = ceil(POT2CHRG(gv.bin[nd]->BandGap-gv.bin[nd]->DustWorkFcn+ONE_ELEC/2.));
			low1 = NINT(help);
		}
		else
		{
			zmin_type zcase = gv.which_zmin[gv.bin[nd]->matType];
			/* >>chng 01 jan 18, the following expressions are taken from Weingartner & Draine, 2001 */
			switch( zcase )
			{
			case ZMIN_CAR:
				help = gv.bin[nd]->AvRadius*1.e7;
				help = ceil(-(1.2*POW2(help)+3.9*help+0.2)/1.44);
				low1 = NINT(help);
/*  			help = POT2CHRG(-0.2866-8.82e5*gv.bin[nd]->AvRadius-1.47e-9/gv.bin[nd]->AvRadius); */
/*  			help = ceil(help) + 1.; */
				break;
			case ZMIN_SIL:
				help = gv.bin[nd]->AvRadius*1.e7;
				help = ceil(-(0.7*POW2(help)+2.5*help+0.8)/1.44);
				low1 = NINT(help);
/*  			help = POT2CHRG(-0.1837-5.15e5*gv.bin[nd]->AvRadius-5.88e-9/gv.bin[nd]->AvRadius); */
/*  			help = ceil(help) + 1.; */
				break;
			case ZMIN_BAKES:
				/* >>refer	grain	physics	Bakes & Tielens, 1994, ApJ, 427, 822 */
				help = ceil(POT2CHRG(gv.bin[nd]->BandGap-gv.bin[nd]->DustWorkFcn+ONE_ELEC/2.));
				low1 = NINT(help);
				break;
			default:
				fprintf( ioQQQ, " InitGrains detected unknown Zmin type: %d\n" , zcase );
				puts( "[Stop in InitGrains]" );
				cdEXIT(1);
			}
		}

		/* this is to assure that gv.bin[nd]->LowestZg > LONG_MIN */
		ASSERT( help > (double)(LONG_MIN+1) );

		/* >>chng 01 apr 20, iterate to get LowestPot such that the exponent in the thermionic
		 * rate never becomes positive; the value can be derived by equating ThresInf >= 0;
		 * the new expression for Emin (see GetPotValues) cannot be inverted analytically,
		 * hence it is necessary to iterate for LowestPot. this also automatically assures that
		 * the expressions for ThresInf and LowestPot are consistent with each other, PvH */
		low2 = low1;
		GetPotValues(nd,low2,&ThresInf,&dum[0],&dum[1],&dum[2],&dum[3],INCL_TUNNEL);
		if( ThresInf < 0. )
		{
			low3 = 0;
			/* do a bisection search for the lowest charge such that
			 * ThresInf >= 0, the end result will eventually be in low3 */
			while( low3-low2 > 1 )
			{
				lowm = (low2+low3)/2;
				GetPotValues(nd,lowm,&ThresInf,&dum[0],&dum[1],&dum[2],&dum[3],INCL_TUNNEL);
				if( ThresInf < 0. )
					low2 = lowm;
				else
					low3 = lowm;
			}
			low2 = low3;
		}

		/* the first term implements the minimum charge due to autoionization
		 * the second term assures that the exponent in the thermionic rate never
		 * becomes positive; the expression was derived by equating ThresInf >= 0 */
		gv.bin[nd]->LowestZg = MAX2(low1,low2);
		gv.bin[nd]->LowestPot = CHRG2POT(gv.bin[nd]->LowestZg);

		/* >>chng 00 jul 13, new sticking probability for electrons */
		/* the second term is chance that electron passes through grain,
		 * 1-p_rad is chance that electron is ejected before grain settles
		 * see discussion in 
		 * >>refer	grain	physics	Weingartner & Draine, 2001, ApJS, 134, 263 */
		gv.bin[nd]->StickElecPos = STICK_ELEC*(1. - exp(-gv.bin[nd]->AvRadius/MEAN_PATH));
		atoms = gv.bin[nd]->AvVol*gv.bin[nd]->dustp[0]/ATOMIC_MASS_UNIT/gv.bin[nd]->atomWeight;
		p_rad = 1./(1.+exp(20.-atoms));
		gv.bin[nd]->StickElecNeg = gv.bin[nd]->StickElecPos*p_rad;

		/* >>chng 02 feb 15, these quantities depend on radius and are normally set
		 * in GrainUpdateRadius1(), however, it is necessary to initialize them here
		 * as well so that they are valid the first time hmole is called. */
		gv.bin[nd]->dstAbund = (float)(gv.bin[nd]->dstfactor*GrnVryDpth(nd));
		ASSERT( gv.bin[nd]->dstAbund > 0.f );
		/* grain unit conversion, <unit>/H (default depl) -> <unit>/cm^3 (actual depl) */
		gv.bin[nd]->cnv_H_pCM3 = dense.gas_phase[ipHYDROGEN]*gv.bin[nd]->dstAbund;
		gv.bin[nd]->cnv_CM3_pH = 1./gv.bin[nd]->cnv_H_pCM3;
		/* grain unit conversion, <unit>/cm^3 (actual depl) -> <unit>/grain */
		gv.bin[nd]->cnv_CM3_pGR = gv.bin[nd]->cnv_H_pGR/gv.bin[nd]->cnv_H_pCM3;
		gv.bin[nd]->cnv_GR_pCM3 = 1./gv.bin[nd]->cnv_CM3_pGR;
	}

	/* >>chng 02 dec 19, these quantities depend on radius and are normally set
	 * in GrainUpdateRadius1(), however, it is necessary to initialize them here
	 * as well so that they are valid for the initial printout in Cloudy, PvH */
	/* calculate the summed grain abundances, these are valid at the inner radius;
	 * these numbers depend on radius and are updated in GrainUpdateRadius1() */
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		gv.elmSumAbund[nelem] = 0.f;
	}

	for( nd=0; nd < gv.nBin; nd++ )
	{
		for( nelem=0; nelem < LIMELM; nelem++ )
		{
			gv.elmSumAbund[nelem] += gv.bin[nd]->elmAbund[nelem]*(float)gv.bin[nd]->cnv_H_pCM3;
		}
	}

	gv.nfill = -1;
	gv.nzone = -1;
	gv.lgAnyNegCharge = FALSE;
	gv.GrnRecomTe = -1.f;

	/* >>chng 01 nov 21, total grain opacities depend on charge and therefore on radius,
	 *                   invalidate for now; GrainUpdateRadius2() will set correct values */
	for( i=0; i < rfield.nupper; i++ )
	{
		/* these are total absorption and scattering cross sections,
		 * the latter should contain the asymmetry factor (1-g) */
		gv.dstab[i] = -DBL_MAX;
		gv.dstsc[i] = -DBL_MAX;
	}

	for( i=0; i < rfield.nupper; i++ )
	{
		gv.anumin[i] = rfield.anu[i] - 0.5f*rfield.widflx[i];
		gv.anumax[i] = rfield.anu[i] + 0.5f*rfield.widflx[i];
	}

	InitEmissivities();
	InitEnthalpy();

	if( trace.lgDustBug && trace.lgTrace )
	{
		fprintf( ioQQQ, "     There are %ld grain types turned on.\n", gv.nBin );

		fprintf( ioQQQ, "     grain depletion factors, dstfactor=" );
		for( nd=0; nd < gv.nBin; nd++ )
		{
			fprintf( ioQQQ, "%10.2e", gv.bin[nd]->dstfactor );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "     nChrg =" );
		for( nd=0; nd < gv.nBin; nd++ )
		{
			fprintf( ioQQQ, " %ld", gv.bin[nd]->nChrg );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "     lowest charge (e) =" );
		for( nd=0; nd < gv.nBin; nd++ )
		{
			fprintf( ioQQQ, "%10.2e", POT2CHRG(gv.bin[nd]->LowestPot) );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "     lgDustVary flag for depth dependence:" );
		for( nd=0; nd < gv.nBin; nd++ )
		{
			fprintf( ioQQQ, "%2c", TorF(gv.bin[nd]->lgDustVary) );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "     Quantum heating flag:" );
		for( nd=0; nd < gv.nBin; nd++ )
		{
			fprintf( ioQQQ, "%2c", TorF(gv.bin[nd]->lgQHeat) );
		}
		fprintf( ioQQQ, "\n" );

		/* >>chng 01 nov 21, removed total abs and sct cross sections, they are invalid */
		fprintf( ioQQQ, "     NU(Ryd), Abs cross sec per proton\n" );

		fprintf( ioQQQ, "    Ryd   " );
		for( nd=0; nd < gv.nBin; nd++ )
		{
			fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
		}
		fprintf( ioQQQ, "\n" );

		for( i=0; i < rfield.nupper; i += 40 )
		{
			fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
			for( nd=0; nd < gv.nBin; nd++ )
			{
				fprintf( ioQQQ, " %10.2e  ", gv.bin[nd]->dstab1[i] );
			}
			fprintf( ioQQQ, "\n" );
		}

		fprintf( ioQQQ, "     NU(Ryd), Sct cross sec per proton\n" );

		fprintf( ioQQQ, "    Ryd   " );
		for( nd=0; nd < gv.nBin; nd++ )
		{
			fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
		}
		fprintf( ioQQQ, "\n" );

		for( i=0; i < rfield.nupper; i += 40 )
		{
			fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
			for( nd=0; nd < gv.nBin; nd++ )
			{
				fprintf( ioQQQ, " %10.2e  ", gv.bin[nd]->pure_sc1[i] );
			}
			fprintf( ioQQQ, "\n" );
		}

		fprintf( ioQQQ, "     NU(Ryd), Q abs\n" );

		fprintf( ioQQQ, "    Ryd   " );
		for( nd=0; nd < gv.nBin; nd++ )
		{
			fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
		}
		fprintf( ioQQQ, "\n" );

		for( i=0; i < rfield.nupper; i += 40 )
		{
			fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
			for( nd=0; nd < gv.nBin; nd++ )
			{
				fprintf( ioQQQ, " %10.2e  ", gv.bin[nd]->dstab1[i]*4./gv.bin[nd]->IntArea );
			}
			fprintf( ioQQQ, "\n" );
		}

		fprintf( ioQQQ, "     NU(Ryd), Q sct\n" );

		fprintf( ioQQQ, "    Ryd   " );
		for( nd=0; nd < gv.nBin; nd++ )
		{
			fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
		}
		fprintf( ioQQQ, "\n" );

		for( i=0; i < rfield.nupper; i += 40 )
		{
			fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
			for( nd=0; nd < gv.nBin; nd++ )
			{
				fprintf( ioQQQ, " %10.2e  ", gv.bin[nd]->pure_sc1[i]*4./gv.bin[nd]->IntArea );
			}
			fprintf( ioQQQ, "\n" );
		}

		fprintf( ioQQQ, "     NU(Ryd), asymmetry factor\n" );

		fprintf( ioQQQ, "    Ryd   " );
		for( nd=0; nd < gv.nBin; nd++ )
		{
			fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
		}
		fprintf( ioQQQ, "\n" );

		for( i=0; i < rfield.nupper; i += 40 )
		{
			fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
			for( nd=0; nd < gv.nBin; nd++ )
			{
				fprintf( ioQQQ, " %10.2e  ", gv.bin[nd]->asym[i] );
			}
			fprintf( ioQQQ, "\n" );
		}

		fprintf( ioQQQ, "     NU(Ryd), yield enhancement\n" );

		fprintf( ioQQQ, "    Ryd   " );
		for( nd=0; nd < gv.nBin; nd++ )
		{
			fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
		}
		fprintf( ioQQQ, "\n" );

		for( i=0; i < rfield.nupper; i += 40 )
		{
			fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
			for( nd=0; nd < gv.nBin; nd++ )
			{
				fprintf( ioQQQ, " %10.2e  ", gv.bin[nd]->y1[i] );
			}
			fprintf( ioQQQ, "\n" );
		}

		fprintf( ioQQQ, " InitGrains exits.\n" );
	}

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

STATIC void InitEmissivities()
{
	double fac,
	  fac2,
	  mul,
	  tdust;
	long int i,
	  nd;

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

	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "  InitEmissivities starts\n" );
		fprintf( ioQQQ, "    ND    Tdust       Emis       BB Check   4pi*a^2*<Q>\n" );
	}

	ASSERT( NTOP >= 2 && NDEMS > 2*NTOP );
	fac = exp(log(GRAIN_TMID/GRAIN_TMIN)/(double)(NDEMS-NTOP));
	tdust = GRAIN_TMIN;
	for( i=0; i < NDEMS-NTOP; i++ )
	{
		gv.dsttmp[i] = log(tdust);
		for( nd=0; nd < gv.nBin; nd++ )
		{
			gv.bin[nd]->dstems[i] = log(PlanckIntegral(tdust,nd,i));
		}
		tdust *= fac;
	}

	/* temperatures above GRAIN_TMID are unrealistic -> make grid gradually coarser */
	fac2 = exp(log(GRAIN_TMAX/GRAIN_TMID/powi(fac,NTOP-1))/(double)((NTOP-1)*NTOP/2));
	for( i=NDEMS-NTOP; i < NDEMS; i++ )
	{
		gv.dsttmp[i] = log(tdust);
		for( nd=0; nd < gv.nBin; nd++ )
		{
			gv.bin[nd]->dstems[i] = log(PlanckIntegral(tdust,nd,i));
		}
		fac *= fac2;
		tdust *= fac;
	}

	/* sanity checks */
	mul = 1.;
	ASSERT( fabs(gv.dsttmp[0] - log(GRAIN_TMIN)) < 10.*mul*DBL_EPSILON );
	mul = sqrt((double)(NDEMS-NTOP));
	ASSERT( fabs(gv.dsttmp[NDEMS-NTOP] - log(GRAIN_TMID)) < 10.*mul*DBL_EPSILON );
	mul = (double)NTOP + sqrt((double)NDEMS);
	ASSERT( fabs(gv.dsttmp[NDEMS-1] - log(GRAIN_TMAX)) < 10.*mul*DBL_EPSILON );

	/* now find slopes form spline fit */
	for( nd=0; nd < gv.nBin; nd++ )
	{
		/* set up coefficients for spline */
		spline(gv.bin[nd]->dstems,gv.dsttmp,NDEMS,2e31,2e31,gv.bin[nd]->dstslp);
		spline(gv.dsttmp,gv.bin[nd]->dstems,NDEMS,2e31,2e31,gv.bin[nd]->dstslp2);
	}

#	if 0
	/* test the dstems interpolation */
	nd = NINT(fudge(0));
	ASSERT( nd >= 0 && nd < gv.nBin );
	for( i=0; i < 2000; i++ )
	{
		double x,y,z;
		z = pow(10.,-40. + (double)i/50.);
		splint(gv.bin[nd]->dstems,gv.dsttmp,gv.bin[nd]->dstslp,NDEMS,log(z),&y);
		if( exp(y) > GRAIN_TMIN && exp(y) < GRAIN_TMAX )
		{
			x = PlanckIntegral(exp(y),nd,3);
			printf(" input %.6e temp %.6e output %.6e rel. diff. %.6e\n",z,exp(y),x,(x-z)/z);
		}
	}
	cdEXIT(EXIT_SUCCESS);
#	endif

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


/* PlanckIntegral compute total radiative cooling due to grains */
STATIC double PlanckIntegral(double tdust, 
			     long int nd, 
			     long int ip)
{
	long int i;

	double arg,
	  ExpM1,
	  integral1 = 0.,  /* integral(Planck) */
	  integral2 = 0.,  /* integral(Planck*abs_cs) */
	  Planck1,
	  Planck2,
	  TDustRyg, 
	  x;

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

	/******************************************************************
	 *
	 * >>>chng 99 mar 12, this sub rewritten following Peter van Hoof
	 * comments.  Original coding was in single precision, and for
	 * very low temperature the exponential was set to zero.  As 
	 * a result Q was far too large for grain temperatures below 10K
	 *
	 ******************************************************************/

	/* Boltzmann factors for Planck integration */
	TDustRyg = TE1RYD/tdust;

	x = 0.999*log(DBL_MAX);

	for( i=0; i < rfield.nupper; i++ )
	{
		/* this is hnu/kT for grain at this temp and photon energy */
		arg = TDustRyg*rfield.anu[i];

		/* want the number exp(hnu/kT) - 1, two expansions */
		if( arg < 1.e-5 )
		{
			/* for small arg expand exp(hnu/kT) - 1 to second order */
			ExpM1 = arg*(1. + arg/2.);
		}
		else
		{
			/* for large arg, evaluate the full expansion */
			ExpM1 = exp(MIN2(x,arg)) - 1.;
		}

		Planck1 = PI4*2.*HPLANCK/POW2(SPEEDLIGHT)*POW2(FR1RYD)*POW2(FR1RYD)*
			rfield.anu3[i]/ExpM1*rfield.widflx[i];
		Planck2 = Planck1*gv.bin[nd]->dstab1[i];

		/* add integral over RJ tail, maybe useful for extreme low temps */
		if( i == 0 ) 
		{
			integral1 = Planck1/rfield.widflx[0]*rfield.anu[0]/3.;
			integral2 = Planck2/rfield.widflx[0]*rfield.anu[0]/5.;
		}
		/* if we are in the Wien tail - exit */
		if( Planck1/integral1 < DBL_EPSILON && Planck2/integral2 < DBL_EPSILON )
			break;

		integral1 += Planck1;
		integral2 += Planck2;
	}

	/* this is an option to print out every few steps, when 'trace grains' is set */
	if( trace.lgDustBug && trace.lgTrace && ip%10 == 0 )
	{
		fprintf( ioQQQ, "  %4ld %11.4e %11.4e %11.4e %11.4e\n", nd, tdust, 
		  integral2, integral1/4./5.67051e-5/powi(tdust,4), integral2*4./integral1 );
	}

	ASSERT( integral2 > 0. );

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


/* invalidate charge dependent data from previous iteration */
STATIC void NewChargeData(long nd)
{
	long nz;

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

	for( nz=0; nz < NCHS; nz++ )
	{
		gv.bin[nd]->chrg[nz]->RSum1 = -DBL_MAX;
		gv.bin[nd]->chrg[nz]->RSum2 = -DBL_MAX;
		gv.bin[nd]->chrg[nz]->ESum1a = -DBL_MAX;
		gv.bin[nd]->chrg[nz]->ESum1b = -DBL_MAX;
		gv.bin[nd]->chrg[nz]->ESum1c = -DBL_MAX;
		gv.bin[nd]->chrg[nz]->ESum2 = -DBL_MAX;

		/*TODO  should any of the following 3 statements be removed ???? */
		gv.bin[nd]->chrg[nz]->ThermRate = -DBL_MAX;
		gv.bin[nd]->chrg[nz]->HeatingRate2 = -DBL_MAX;
		gv.bin[nd]->chrg[nz]->GrainHeat = -DBL_MAX;

		gv.bin[nd]->chrg[nz]->hots1 = -DBL_MAX;
		gv.bin[nd]->chrg[nz]->bolflux1 = -DBL_MAX;
		gv.bin[nd]->chrg[nz]->pe1 = -DBL_MAX;
	}

	if( fabs(phycon.te-gv.GrnRecomTe) > 10.f*FLT_EPSILON*gv.GrnRecomTe )
	{
		for( nz=0; nz < NCHS; nz++ )
		{
			memset( gv.bin[nd]->chrg[nz]->eta, 0, (LIMELM+2)*sizeof(double) );
			memset( gv.bin[nd]->chrg[nz]->xi, 0, (LIMELM+2)*sizeof(double) );
		}
	}

	if( nzone != gv.nzone )
	{
		for( nz=0; nz < NCHS; nz++ )
		{
			gv.bin[nd]->chrg[nz]->hcon1 = -DBL_MAX;
		}
	}

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


/* this is the main routine that drives the grain physics */
void GrainDrive(void)
{
#	ifdef DEBUG_FUN
	fputs( "<+>GrainDrive()\n", debug_fp );
#	endif

	/* gv.lgGrainPhysicsOn set false with no grain physics command */
	if( gv.lgDustOn && gv.lgGrainPhysicsOn )
	{
		static float tesave=-1.f;
		static long int nzonesave=-1;

		/* option to only reevaluate grain physics if something has changed.  
		 * gv.lgReevaluate is set false with keyword no reevaluate on grains command 
		 * option to force constant reevaluation of grain physics - 
		 * by default is TRUE 
		 * usually reevaluate grains at all times, but NO REEVALUATE will
		 * save some time but may affect stability */
		if( gv.lgReevaluate || conv.lgSearch || nzonesave != nzone || 
			/* need to reevaluate the grains when temp changes since */
			/*lint -e777 float test equality */
			phycon.te !=tesave || 
			/* >>chng 03 dec 30, check that electrons locked in grains are not important,
			 * if they are, then reevaluate */
			 fabs(gv.TotalEden)/dense.eden > conv.EdenErrorAllowed/5. ||
			 /* >>chng 04 aug 06, always reevaluate when thermal effects of grains are important,
			  * first is collisional energy exchange with gas, second is grain photoionization */
			 (fabs( gv.GasCoolColl ) + fabs( thermal.heating[0][13] ))/SDIV(thermal.ctot)>0.1 )
			/*lint +e777 float test equality */
			/*>>chng 03 oct 12, from not exactly same to 1% */
			/*fabs(phycon.te-tesave)/phycon.te>0.01 )*/
		{
			nzonesave = nzone;
			tesave = phycon.te;

			if( trace.lgTrConvg >= 5 )
			{
				fprintf( ioQQQ, "     grain5 calling GrainChargeTemp and GrainDrift\n");
			}
			/* find dust charge and temperature - this must be called at least once per zone
			 * since grain abundances, set here, may change with depth */
			GrainChargeTemp();

			/* >>chng 04 jan 31, moved call to GrainDrift to ConvPresTempEdenIoniz(), PvH */
		}
	}
	else if( gv.lgDustOn && !gv.lgGrainPhysicsOn )
	{
		/* very minimalistic treatment of grains; only extinction of continuum is considered
		 * however, the absorbed energy is not reradiated, so this creates thermal imbalance! */
		if( nCalledGrainDrive == 0 )
		{
			long nelem, ion, ion_to, nd;

			/* when not doing grain physics still want some exported quantities
			 * to be reasonable, grain temperature used for H2 formation */
			for( nd=0; nd < gv.nBin; nd++ )
			{
				long nz;

				/* this disables warnings about PAHs in the ionized region */
				gv.bin[nd]->lgPAHsInIonizedRegion = FALSE;

				for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
				{
					gv.bin[nd]->chrg[nz]->DustZ = nz;
					gv.bin[nd]->chrg[nz]->FracPop = ( nz == 0 ) ? 1. : 0.;
					gv.bin[nd]->chrg[nz]->nfill = 0;
					gv.bin[nd]->chrg[nz]->tedust = 100.f;
				}

				gv.bin[nd]->AveDustZ = 0.;
				gv.bin[nd]->dstpot = CHRG2POT(0.);

				gv.bin[nd]->tedust = 100.f;
				gv.bin[nd]->TeGrainMax = 100.;

				/* set all heating/cooling agents to zero */
				gv.bin[nd]->BolFlux = 0.;
				gv.bin[nd]->GrainCoolTherm = 0.;
				gv.bin[nd]->GasHeatPhotoEl = 0.;
				gv.bin[nd]->GrainHeat = 0.;
				gv.bin[nd]->GrainHeatColl = 0.;
				gv.bin[nd]->ChemEn = 0.;
				gv.bin[nd]->ChemEnH2 = 0.;
				gv.bin[nd]->thermionic = 0.;

				gv.bin[nd]->lgUseQHeat = FALSE;
				gv.bin[nd]->lgEverQHeat = FALSE;
				gv.bin[nd]->QHeatFailures = 0;

				gv.bin[nd]->DustDftVel = 0.;

				gv.bin[nd]->avdust = gv.bin[nd]->tedust;
				gv.bin[nd]->avdft = 0.f;
				gv.bin[nd]->avdpot = (float)(gv.bin[nd]->dstpot*EVRYD);
				gv.bin[nd]->avDGRatio = -1.f;
			}

			gv.lgAnyNegCharge = FALSE;

			gv.TotalEden = 0.;
			gv.GrnElecDonateMax = 0.f;
			gv.GrnElecHoldMax = 0.f;

			for( nelem=0; nelem < LIMELM; nelem++ )
			{
				for( ion=0; ion <= nelem+1; ion++ )
				{
					for( ion_to=0; ion_to <= nelem+1; ion_to++ )
					{
						gv.GrainChTrRate[nelem][ion][ion_to] = 0.f;
					}
				}
			}

			/* set all heating/cooling agents to zero */
			gv.GrainHeatInc = 0.;
			gv.GrainHeatDif = 0.;
			gv.GrainHeatLya = 0.;
			gv.GrainHeatCollSum = 0.;
			gv.GrainHeatSum = 0.;
			gv.GrainHeatChem = 0.;
			gv.GasCoolColl = 0.;
			gv.TotalDustHeat = 0.f;
			gv.dphmax = 0.f;
			gv.dclmax = 0.f;

			thermal.heating[0][13] = 0.;
			thermal.heating[0][14] = 0.;
			thermal.heating[0][25] = 0.;
		}

		if( nCalledGrainDrive == 0 || gv.lgAnyDustVary )
		{
			GrainUpdateRadius1();
			GrainUpdateRadius2(FALSE);
		}
	}

	++nCalledGrainDrive;

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

/* iterate grain charge and temperature */
STATIC void GrainChargeTemp(void)
{
	int lgAnyNegCharge = FALSE;
	long int i,
	  ion,
	  ion_to,
	  nelem,
	  nd,
	  nz;
	float dccool = FLT_MAX;
	double delta,
	  GasHeatNet,
	  hcon = DBL_MAX,
	  hla = DBL_MAX,
	  hots = DBL_MAX,
	  oldtemp,
	  oldTotalEden,
	  oldGasHeatNet,
	  ratio,
	  ThermRatio;

	static long int oldZone = -1;
	static double oldTe = -DBL_MAX,
	  oldHeat = -DBL_MAX;

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

	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "\n GrainChargeTemp called lgSearch%2c\n\n", TorF(conv.lgSearch) );
	}

	oldTotalEden = gv.TotalEden;
	oldGasHeatNet = gv.GasHeatNet;

	/* these will sum heating agents over grain populations */
	gv.GrainHeatInc = 0.;
	gv.GrainHeatDif = 0.;
	gv.GrainHeatLya = 0.;
	gv.GrainHeatCollSum = 0.;
	gv.GrainHeatSum = 0.;
	gv.GrainHeatChem = 0.;

	gv.GasCoolColl = 0.;
	gv.GasHeatPhotoEl = 0.;
	gv.GasHeatTherm = 0.;

	gv.TotalEden = 0.;

	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		for( ion=0; ion <= nelem+1; ion++ )
		{
			for( ion_to=0; ion_to <= nelem+1; ion_to++ )
			{
				gv.GrainChTrRate[nelem][ion][ion_to] = 0.f;
			}
		}
	}

	gv.HighestIon = HighestIonStage();

	/* this sets dstAbund and conversion factors, but not gv.dstab and gv.dstsc! */
	GrainUpdateRadius1();

	for( nd=0; nd < gv.nBin; nd++ )
	{
		double one;
		long relax = ( conv.lgSearch ) ? 3 : 1;
#ifdef OLD_POT_ALGO
		double tol1 = TOLER;
#endif

		/* >>chng 02 nov 11, added test for the presence of PAHs in the ionized region, PvH */
		if( gv.bin[nd]->matType == MAT_PAH || gv.bin[nd]->matType == MAT_PAH2 )
		{
			if( dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] > 0.50 )
			{
				gv.bin[nd]->lgPAHsInIonizedRegion = TRUE;
			}
		}

		/* >>chng 01 sep 13, dynamically allocate backup store, remove ncell dependence, PvH */
		/* allocate data inside loop to avoid accidental spillover to next iteration */
		/* >>chng 04 jan 18, no longer delete and reallocate data inside loop to speed up the code, PvH */
		NewChargeData(nd);

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

		delta = 2.*TOLER;
		/* >>chng 01 nov 29, relax max no. of iterations during initial search */
		for( i=0; i < relax*CT_LOOP_MAX && delta > TOLER; ++i )
		{
			long j;
			double TdBracketLo = 0., TdBracketHi = -DBL_MAX;
			double ThresEst = 0.;
			oldtemp = gv.bin[nd]->tedust;

			/* solve for charge using previous estimate for grain temp
			 * grain temp only influences thermionic emissions
			 * Thermratio is fraction thermionic emissions contribute
			 * to the total electron loss rate of the grain */
#ifdef OLD_POT_ALGO
			GrainCharge(nd,&ThermRatio,tol1);
#else
			GrainCharge(nd,&ThermRatio);
#endif

			ASSERT( gv.bin[nd]->GrainHeat > 0. );
			ASSERT( gv.bin[nd]->tedust >= GRAIN_TMIN && gv.bin[nd]->tedust <= GRAIN_TMAX );

			/* >>chng 04 may 31, in conditions where collisions become an important
			 * heating/cooling source (e.g. gas that is predominantly heated by cosmic
			 * rays), the heating rate depends strongly on the assumed dust temperature.
			 * hence it is necessary to iterate for the dust temperature. PvH */
			gv.bin[nd]->lgTdustConverged = FALSE;
			for( j=0; j < relax*T_LOOP_MAX; ++j )
			{
				double oldTemp2 = gv.bin[nd]->tedust;
				double oldHeat = gv.bin[nd]->GrainHeat;
				double oldCool = gv.bin[nd]->GrainGasCool;

				/* now solve grain temp using new value for grain potential */
				GrainTemperature(nd,&dccool,&hcon,&hots,&hla);

				gv.bin[nd]->GrainGasCool = dccool;

				if( trace.lgTrace && trace.lgDustBug )
				{
					fprintf( ioQQQ, "  >>loop %ld BracketLo %.6e BracketHi %.6e",
						 j, TdBracketLo, TdBracketHi );
				}

				/* this test assures that convergence can only happen if GrainHeat > 0
				 * and therefore the value of tedust is guaranteed to be valid as well */
				/* >>chng 04 aug 05, test that gas cooling is converged as well,
				 * in deep PDRs gas cooling depends critically on grain temperature, PvH */
				if( fabs(gv.bin[nd]->GrainHeat-oldHeat) < HEAT_TOLER*gv.bin[nd]->GrainHeat &&
				    fabs(gv.bin[nd]->GrainGasCool-oldCool) < HEAT_TOLER_BIN*thermal.ctot )
				{
					gv.bin[nd]->lgTdustConverged = TRUE;
					if( trace.lgTrace && trace.lgDustBug )
						fprintf( ioQQQ, " converged\n" );
					break;
				}

				/* update the bracket for the solution */
				if( gv.bin[nd]->tedust < oldTemp2 )
					TdBracketHi = oldTemp2;
				else
					TdBracketLo = oldTemp2;

				/* GrainTemperature yields a new estimate for tedust, and initially
				 * that estimate will be used. In most zones this will converge quickly.
				 * However, sometimes the solution will oscillate and converge very
				 * slowly. So, as soon as j >= 2 and the bracket is set up, we will
				 * force convergence by using a bisection search within the bracket */
				/*TODO   this algorithm might be more efficient with Brent */

				/* this test assures that TdBracketHi is initialized */
				if( TdBracketHi > TdBracketLo )
				{
					/* if j >= 2, the solution is converging too slowly
					 * so force convergence by doing a bisection search */
					if( ( j >= 2 && TdBracketLo > 0. ) ||
					    gv.bin[nd]->tedust <= TdBracketLo ||
					    gv.bin[nd]->tedust >= TdBracketHi )
					{
						gv.bin[nd]->tedust = (float)(0.5*(TdBracketLo + TdBracketHi));
						if( trace.lgTrace && trace.lgDustBug )
							fprintf( ioQQQ, " bisection\n" );
					}
					else
					{
						if( trace.lgTrace && trace.lgDustBug )
							fprintf( ioQQQ, " iteration\n" );
					}
				}
				else
				{
					if( trace.lgTrace && trace.lgDustBug )
						fprintf( ioQQQ, " iteration\n" );
				}

				ASSERT( gv.bin[nd]->tedust >= GRAIN_TMIN && gv.bin[nd]->tedust <= GRAIN_TMAX );
			}

			if( !gv.bin[nd]->lgTdustConverged )
			{
				int lgBoundErr;
				double y, x = log(gv.bin[nd]->tedust);
				/* make sure GrainHeat is consistent with value of tedust */
				splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,x,&y,&lgBoundErr);
				gv.bin[nd]->GrainHeat = exp(y)*gv.bin[nd]->cnv_H_pCM3;

				fprintf( ioQQQ," PROBLEM  temperature of grain species %s not converged\n",
					 gv.bin[nd]->chDstLab );
				ConvFail("grai","");
			}

			ASSERT( gv.bin[nd]->GrainHeat > 0. );
			ASSERT( gv.bin[nd]->tedust >= GRAIN_TMIN && gv.bin[nd]->tedust <= GRAIN_TMAX );

			/* delta estimates relative change in electron emission rate
			 * due to the update in the grain temperature, if it is small
			 * we won't bother to iterate (which is usually the case)
			 * the formula assumes that thermionic emission is the only
			 * process that depends on grain temperature */
			ratio = gv.bin[nd]->tedust/oldtemp;
			for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
			{
				ThresEst += gv.bin[nd]->chrg[nz]->FracPop*gv.bin[nd]->chrg[nz]->ThresInf;
			}
			delta = ThresEst*TE1RYD/gv.bin[nd]->tedust*(ratio - 1.);
			/*TODO use something like log(ThermRatio) + log(delta) ???? */
			delta = ( delta < 0.9*log(DBL_MAX) ) ?
				ThermRatio*fabs(POW2(ratio)*exp(delta)-1.) : DBL_MAX;
#ifdef OLD_POT_ALGO
			/* >>chng 01 jan 19, increase accuracy of dstpot solution if convergence fails.
			 * when grain charge is close to LowestPot, thermionic rates become important
			 * even for low grain temperatures. under those conditions the rate is a very
			 * sensitive function of the grain potential. if tol1 is too large, ThermRatio
			 * will be inaccurate and delta will blow up, PvH */
			if( i > 0 )
				tol1 /= 10.;
#endif

			if( trace.lgTrace && trace.lgDustBug )
			{
				fprintf( ioQQQ, " >>GrainChargeTemp finds delta=%.4e, ",delta );
				fprintf( ioQQQ, " old/new temp=%.5e %.5e, ",oldtemp,gv.bin[nd]->tedust );
				if( delta > TOLER ) 
				{
					fprintf( ioQQQ, "doing another iteration\n" );
				}
				else 
				{
					fprintf( ioQQQ, "converged\n" );
				}
			}
		}
		if( delta > TOLER )
		{
			fprintf( ioQQQ, " PROBLEM  charge/temperature not converged for %s zone %.2f\n",
				 gv.bin[nd]->chDstLab , fnzone );
			ConvFail("grai","");
		}

		for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
		{
			if( gv.bin[nd]->chrg[nz]->DustZ <= -1 )
				lgAnyNegCharge = TRUE;
		}

		/* add in ion recombination rates on this grain species */
		/* ionbal.lgGrainIonRecom is 1 by default, set to 0 with
		 * no grain neutralization command */
		if( ionbal.lgGrainIonRecom )
			GrainChrgTransferRates(nd);

		/* >>chng 04 jan 31, moved call to UpdateRadius2 outside loop, PvH */

		/* following used to keep track of heating agents in printout
		 * no physics done with GrainHeatInc
		 * dust heating by incident continuum, and elec friction before ejection */
		gv.GrainHeatInc += hcon;
		/* remember total heating by diffuse fields, for printout (includes Lya) */
		gv.GrainHeatDif += hots;
		/* GrainHeatLya - total heating by LA in this zone, erg cm-3 s-1, only here
		 * for eventual printout, hots is total ots line heating */
		gv.GrainHeatLya += hla;

		/* this will be total collisional heating, for printing in lines */
		gv.GrainHeatCollSum += gv.bin[nd]->GrainHeatColl;

		/* GrainHeatSum is total heating of all grain types in this zone,
		 * will be carried by total cooling, only used in lines to print tot heat
		 * printed as entry "GraT    0 " */
		gv.GrainHeatSum += gv.bin[nd]->GrainHeat;

		/* net amount of chemical energy donated by recombining ions and molecule formation */
		gv.GrainHeatChem += gv.bin[nd]->ChemEn + gv.bin[nd]->ChemEnH2;

		/* dccool is gas cooling due to collisions with grains */
		gv.GasCoolColl += dccool;
		gv.GasHeatPhotoEl += gv.bin[nd]->GasHeatPhotoEl;
		gv.GasHeatTherm += gv.bin[nd]->thermionic;

		/* this is grain charge in e/cm^3, positive number means grain supplied free electrons */
		/* >>chng 01 mar 24, changed DustZ+1 to DustZ, PvH */
		one = 0.;
		for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
		{
			one += gv.bin[nd]->chrg[nz]->FracPop*(double)gv.bin[nd]->chrg[nz]->DustZ*
				gv.bin[nd]->cnv_GR_pCM3;
		}
		gv.TotalEden += one;
		{
			/*@-redef@*/
			enum {DEBUG_LOC=FALSE};
			/*@+redef@*/
			if( DEBUG_LOC )
			{
				fprintf(ioQQQ," DEBUG grn chr nz\t%.2f\teden\t%.3e\tnd\t%li",
					fnzone,
					dense.eden,
					nd);
				fprintf(ioQQQ,"\tne\t%.2e\tAveDustZ\t%.2e\t%.2e\t%.2e\t%.2e",
					one,
					gv.bin[nd]->AveDustZ,
					gv.bin[nd]->chrg[0]->FracPop,(double)gv.bin[nd]->chrg[0]->DustZ,
					gv.bin[nd]->cnv_GR_pCM3);
				fprintf(ioQQQ,"\n");
			}
		}

		if( trace.lgTrace && trace.lgDustBug )
		{
			fprintf(ioQQQ,"     %s Pot %.5e Thermal %.5e GasCoolColl %.5e" , 
				gv.bin[nd]->chDstLab, gv.bin[nd]->dstpot, gv.bin[nd]->GrainHeat, dccool );
			fprintf(ioQQQ," GasPEHeat %.5e GasThermHeat %.5e ChemHeat %.5e\n\n" , 
				gv.bin[nd]->GasHeatPhotoEl, gv.bin[nd]->thermionic, gv.bin[nd]->ChemEn );
		}
	}

	/*fprintf(ioQQQ,"DEBUG eden grn %.2f %e %e %e %e %e\n",
		fnzone ,
		dense.eden,
		phycon.te,
		gv.TotalEden,
		gv.bin[20]->dstpot,
		gv.GrainChTrRate[ipSILICON][1][0] );*/

	/* >>chng 04 aug 06, added test of convergence of the net gas heating/cooling, PvH */
	GasHeatNet = gv.GasHeatPhotoEl + gv.GasHeatTherm - gv.GasCoolColl;

	if( fabs(phycon.te-gv.GrnRecomTe) > 10.f*FLT_EPSILON*gv.GrnRecomTe )
	{
		oldZone = gv.nzone;
		oldTe = gv.GrnRecomTe;
		oldHeat = gv.GasHeatNet;
	}

	/* >>chng 04 aug 07, added estimate for heating derivative, PvH */
	if( nzone == oldZone && fabs(phycon.te-oldTe) > 10.f*FLT_EPSILON*phycon.te )
	{
		gv.dHeatdT = (GasHeatNet-oldHeat)/(phycon.te-oldTe);
	}

	/* >>chng 04 sep 15, add test for convergence of gv.TotalEden, PvH */
	if( nzone != gv.nzone || fabs(phycon.te-gv.GrnRecomTe) > 10.f*FLT_EPSILON*gv.GrnRecomTe ||
		fabs(gv.GasHeatNet-oldGasHeatNet) > HEAT_TOLER*thermal.ctot ||
		fabs(gv.TotalEden-oldTotalEden) > CHRG_TOLER*dense.eden )
	{
		/* >>chng 04 aug 07, add test whether eden on grain converged */
		/* flag that change in eden was too large */
		/*conv.lgConvEden = FALSE;*/
		conv.lgConvIoniz = FALSE;
		if( fabs(gv.TotalEden-oldTotalEden) > CHRG_TOLER*dense.eden )
		{
			strcpy( conv.chConvIoniz, "grn eden chg" );
			conv.BadConvIoniz[0] = oldTotalEden;
			conv.BadConvIoniz[1] = gv.TotalEden;
		}
		else if( fabs(gv.GasHeatNet-oldGasHeatNet) > HEAT_TOLER*thermal.ctot )
		{
			strcpy( conv.chConvIoniz, "grn het chg" );
			conv.BadConvIoniz[0] = oldGasHeatNet;
			conv.BadConvIoniz[1] = gv.GasHeatNet;
		}
		else if( fabs(phycon.te-gv.GrnRecomTe) > 10.f*FLT_EPSILON*gv.GrnRecomTe )
		{
			strcpy( conv.chConvIoniz, "grn ter chg" );
			conv.BadConvIoniz[0] = gv.GrnRecomTe;
			conv.BadConvIoniz[1] = phycon.te;
		}
		else if( nzone != gv.nzone )
		{
			strcpy( conv.chConvIoniz, "grn zon chg" );
			conv.BadConvIoniz[0] = gv.nzone;
			conv.BadConvIoniz[1] = nzone;
		}
		else
			TotalInsanity();
	}

	/*printf( "DEBUG GasHeatNet %.6e -> %.6e TotalEden %e -> %e conv.lgConvIoniz %c\n",
		gv.GasHeatNet, GasHeatNet, gv.TotalEden, oldTotalEden, TorF(conv.lgConvIoniz) );*/

/* 	printf( "DEBUG %.2f %e %e\n", fnzone, phycon.te, dense.eden ); */

	gv.nzone = nzone;
	gv.GrnRecomTe = phycon.te;
	gv.GasHeatNet = GasHeatNet;

	/* update total grain opacities in gv.dstab and gv.dstsc,
	 * they depend on grain charge and may depend on depth
	 * also add in the photo-dissociation cs in gv.dstab */
	GrainUpdateRadius2(lgAnyNegCharge);

#ifdef WD_TEST2
	printf("wd test: proton fraction %.5e Total DustZ %.6f heating/cooling rate %.5e %.5e\n",
	       dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN],
	       gv.bin[0]->AveDustZ,gv.GasHeatPhotoEl/dense.gas_phase[ipHYDROGEN]/fudge(0),
	       gv.GasCoolColl/dense.gas_phase[ipHYDROGEN]/fudge(0));
#endif

	if( trace.lgTrace )
	{
		/*@-redef@*/
		enum {DEBUG_LOC=TRUE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			fprintf( ioQQQ, "     %.2f Grain surface charge transfer rates\n", fnzone );

			for( nelem=0; nelem < LIMELM; nelem++ )
			{
				if( dense.lgElmtOn[nelem] )
				{
					printf( "      %s:", elementnames.chElementSym[nelem] );
					for( ion=dense.IonLow[nelem]; ion <= dense.IonHigh[nelem]; ++ion )
					{
						for( ion_to=0; ion_to <= nelem+1; ion_to++ )
						{
							if( gv.GrainChTrRate[nelem][ion][ion_to] > 0.f )
							{
								printf( "  %ld->%ld %.2e", ion, ion_to,
									gv.GrainChTrRate[nelem][ion][ion_to] );
							}
						}
					}
					printf( "\n" );
				}
			}
		}

		fprintf( ioQQQ, "     %.2f Grain contribution to electron density %.2e\n", 
			fnzone , gv.TotalEden );

		fprintf( ioQQQ, "     Grain electons: " );
		for( nd=0; nd < gv.nBin; nd++ )
		{
			double sum = 0.;
			for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
			{
				sum += gv.bin[nd]->chrg[nz]->FracPop*(double)gv.bin[nd]->chrg[nz]->DustZ*
					gv.bin[nd]->cnv_GR_pCM3;
			}
			fprintf( ioQQQ, " %.2e", sum );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "     Grain potentials:" );
		for( nd=0; nd < gv.nBin; nd++ )
		{
			fprintf( ioQQQ, " %.2e", gv.bin[nd]->dstpot );
		}
		fprintf( ioQQQ, "\n" );
		
		fprintf( ioQQQ, "     Grain temperatures:" );
		for( nd=0; nd < gv.nBin; nd++ )
		{
			fprintf( ioQQQ, " %.2e", gv.bin[nd]->tedust );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "     GrainCollCool: %.6e\n", gv.GasCoolColl );
	}

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


#ifdef OLD_POT_ALGO
#define SORT2(X1,X2,F1,F2,T1,T2,N1,N2) \
	if( X1 > X2 )                  \
	{                              \
		double TMP = X1;       \
		X1 = X2;               \
		X2 = TMP;              \
		TMP = F1;              \
		F1 = F2;               \
		F2 = TMP;              \
		TMP = T1;              \
		T1 = T2;               \
		T2 = TMP;              \
		TMP = N1;              \
		N1 = N2;               \
		N2 = TMP;              \
	}


STATIC void GrainCharge(long int nd,
			/*@out@*/double *ThermRatio, /* ratio of thermionic to total rate */
			double toler)
{
	int lgBigError;
	long int lopbig,
	  luplim,
	  nz;
	double dummy,
	  emission=0.,
	  fhi,
	  flo,
	  fmid,
	  help,
	  nhi=0.,
	  nlo=0.,
	  nmid=0.,
	  old_tol,
	  recombination=0.,
	  renorm,
	  thhi,
	  thlo,
	  thmid,
	  tol=DBL_MAX,
	  xhi,
	  xlo,
	  xmid;

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


	/* find dust charge */
	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "  Starting Charge loop for %s, search%2c\n",
			 gv.bin[nd]->chDstLab, TorF(conv.lgSearch) );
	}

	gv.bin[nd]->lgChrgConverged = FALSE;

	if( conv.lgSearch )
	{
		/* this is limit to number of iterations for getting the grain charge */
		luplim = 100;
		/* set up initial estimates for grain potential, will be checked by BracketGrainPot */
		xlo = gv.bin[nd]->LowestPot;
		/* >>chng 03 jan 21, old expression gave problems for extremely small grains, PvH */
		/* xmid = 0.; */
		/* xhi = BIG_POTENTIAL; */
		xmid = gv.bin[nd]->LowestPot + 0.5*BIG_POTENTIAL;
		xhi = gv.bin[nd]->LowestPot + BIG_POTENTIAL;
	}
	else
	{
		/* >>chng 00 jul 19, convergence around zero pot can be tricky; luplim 10 -> 20, PvH */
		/* >>chng 01 jun 02, due to increased demands on precision, change luplim 20 -> 40, PvH */
		/* >>chng 03 jan 24, test case will.in prompted a further increase of luplim 40 -> 50, PvH */
		luplim = 50;
		xlo = MAX2(gv.bin[nd]->LowestPot,gv.bin[nd]->dstpot - BRACKET_STEP);
		xmid = gv.bin[nd]->dstpot;
		xhi = gv.bin[nd]->dstpot + BRACKET_STEP;
	}

	/* >>chng 00 jul 18, keep solution bracketed, this is needed since
	 * bandgap introduces discontinuities in emission and recombination rate;
	 * solution can keep oscillating around zero, PvH */
	/* >>chng 01 jan 08, validate bracket and alter if
	 * necessary, it sets first new step in gv.bin[nd]->dstpot 
	 * (usually the value from the previous zone) */
	BracketGrainPot(nd,&xlo,&flo,&thlo,&nlo,&xmid,&fmid,&thmid,&nmid,&xhi,&fhi,&thhi,&nhi);

	/* if BracketGrainPot accidentally stumbled on the
	 * correct solution, it will be in xmid */

	if( trace.lgTrace && trace.lgDustBug )
	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			long i;
			double backup = gv.bin[nd]->dstpot;
			/* switch off trace output temporarily */
			trace.lgDustBug = FALSE;
			for( i=0; i<=100; i++ )
			{
				gv.bin[nd]->dstpot = xlo + (double)i/100.*(xhi-xlo);
				UpdatePot(nd,1);
				emission = GrainElecEmis(nd,&dummy);
				recombination = GrainElecRecomb(nd);
				printf(" GrainPotTest: %.6e %.6e %.6e\n",
				       POT2CHRG(gv.bin[nd]->dstpot),emission,recombination );
			}
			/* reset to old values */
			gv.bin[nd]->dstpot = backup;
			UpdatePot(nd,1);
			trace.lgDustBug = TRUE;
		}
	}

	lgBigError = TRUE;

	help = fabs(gv.bin[nd]->LowestPot);
	renorm = MAX2(help,1.);

	/* >>chng 01 jan 09, completely rewritten optimization loop. it now uses
	 * Brent's algorithm instead of Newton-Raphson which is a more stable
	 * choice if the derivative is not known explicitly, PvH */
	for( lopbig=0; lopbig < luplim; lopbig++ )
	{
		int lgSlow;
		double delta,
		  fnew,
		  minstep,
		  nnew,
		  thnew,
		  tol1,
		  tol2,
		  xnew;

		char* method;

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

		old_tol = tol;
		tol = ( flo*fmid > 0. ) ? xhi-xmid : xlo-xmid;
		tol2 = ( nmid > 0. ) ? fmid/nmid : 0.;
		lgSlow = fabs(tol/old_tol) > 0.5;

		/* >>chng 01 jan 17, using this instead of tol1 = toler assures that the potential
		 * is very accurate near lowest potential where the emission rate is very steep, PvH */
		/* >>chng 01 may 31, change 10.*DBL_EPSILON -> 10.*DBL_EPSILON*xmid, PvH */
		tol1 = MAX2(toler*(xmid-gv.bin[nd]->LowestPot)/renorm,10.*DBL_EPSILON*xmid);

		if( DEBUG_LOC )
		{
			int lgTest = ( fabs(tol) < tol1 );
			printf("   tol %.15e tol1 %.15e -> %c\n",tol,tol1,TorF(lgTest));
			lgTest = ( fabs(tol2) < toler );
			printf("   tol2 %.15e toler %.15e -> %c\n",tol2,toler,TorF(lgTest));
			lgTest = ( fmid == 0. );
			printf("   fmid %.15e nmid %.15e -> %c\n",fmid,nmid,TorF(lgTest));
			lgTest = ( fmid == 0. || ( fabs(tol) < tol1 && fabs(tol2) < toler ));
			printf("   converged ? --> %c\n",TorF(lgTest));
		}

		/* >>chng 01 may 31, put in additional check to assure that emission equals recombination
		 * this is important for bandgap materials with average charges between -1 and 0; the
		 * emission rates can be strongly non-linear and lead to significant discrepancies, PvH */
		if( fmid == 0. || ( fabs(tol) < tol1 && fabs(tol2) < toler ) )
		{
			lgBigError = FALSE;
			break;
		}

		if( ! lgSlow )
		{
			double r = fmid/fhi;
			double s = fmid/flo;
			double t = flo/fhi;
			double p = s*((1.-r)*(xmid-xlo) - t*(r-t)*(xhi-xmid));
			double q = (1.-t)*(r-1.)*(s-1.);
			/* we are only interested in the ratio p/q, make sure q is positive */
			if( q < 0. )
			{
				p = -p;
				q = -q;
			}
			/* inverse quadratic interpolation is within bounds */
			if( p > q*(xlo-xmid) && p < q*(xhi-xmid) )
			{
				delta = p/q;
				method = "InvQuadr";
			}
			/* inverse quadratic interpolation is out of bounds, use bisection */
			else
			{
				delta = 0.5*tol;
				method = "OBBisect";
			}
		}
		/* progress for inverse quadratic interpolation is too slow, use bisection */
		else
		{
			delta = 0.5*tol;
			method = "SLBisect";
		}

		/* >>chng 00 jul 14, make sure change is not infinitesimally small, PvH */
		/* >>chng 01 may 31, prevent stepsize from becoming too small */
		help = fabs(xmid);
		minstep = 5.*DBL_EPSILON*MAX2(help,1.);
		delta = ( delta < 0. ) ? MIN2(delta,-minstep) : MAX2(delta,minstep);

		/* >>chng 01 jan 08, removed upper limit on change, code above handles this, PvH */
		xnew = MAX2(xmid + delta,gv.bin[nd]->LowestPot);

		if( DEBUG_LOC )
		{
			printf(" xlo %.15e xmid %.15e xhi %.15e delta %.15e xnew %.15e\n",xlo,xmid,xhi,delta,xnew);
			printf(" flo %.15e fmid %.15e fhi %.15e\n",flo,fmid,fhi);
		}

		if( xnew < xlo || xnew > xhi )
		{
			fprintf( ioQQQ, " GrainCharge failed to keep solution bracketed.\n" );
			lgBigError = TRUE;
			break;
		}

		if( trace.lgTrace && trace.lgDustBug )
		{
			fprintf( ioQQQ, "  Charge loop %ld %s Vmid %.4f (%.1e) Vnext %.4f Vlo %.4f (%.1e)", 
			  lopbig, gv.bin[nd]->chDstLab, xmid, fmid, xnew, xlo, flo );
			fprintf( ioQQQ, " Vhi %.4f (%.1e) %s tol %.4f Em/Rc %.3e %.3e\n", 
			  xhi, fhi, method, fabs(tol), emission, recombination );
		}

		gv.bin[nd]->dstpot = xnew;
		UpdatePot(nd,1);

		/* function GrainElecEmis returns electron loss rate for the grain species
		 * due to photoelectric effect, positive ion recombinations on the
		 * grain surface and thermionic emissions */
		emission = GrainElecEmis(nd,&thnew);

		/* the function GrainElecRecomb returns the recombination rate
		 * for this grain and charge */
		recombination = GrainElecRecomb(nd);

		/* this is the difference between emission and recombination
		 * rates - we want this to be small */
		fnew = emission - recombination;
		nnew = 0.5*(emission+recombination);

		if( DEBUG_LOC )
			printf(" xnew %.15e fnew %.15e nnew %.15e\n",xnew,fnew,nnew);

		/* update bracket for the correct solution */
		if( fnew == 0. )
		{
			/* this sometimes happens, really ! */
			xlo = xmid = xhi = xnew;
			flo = fmid = fhi = fnew;
			thlo = thmid = thhi = thnew;
			nlo = nmid = nhi = nnew;
		}
		else if( fnew*fmid > 0. )
		{
			/* there are three points on one side of zero point -> remove farthest value */
			if( flo*fnew > 0. )
			{
				xlo = xnew;
				flo = fnew;
				thlo = thnew;
				nlo = nnew;
				SORT2(xlo,xmid,flo,fmid,thlo,thmid,nlo,nmid);
			}
			else
			{
				xhi = xnew;
				fhi = fnew;
				thhi = thnew;
				nhi = nnew;
				SORT2(xmid,xhi,fmid,fhi,thmid,thhi,nmid,nhi);
			}
		}
		else
		{
			/* there are two points on either side of zero point -> minimize bracket size */
			if( MIN2(xnew,xmid)-xlo > xhi-MAX2(xnew,xmid) )
			{
				xlo = xnew;
				flo = fnew;
				thlo = thnew;
				nlo = nnew;
				SORT2(xlo,xmid,flo,fmid,thlo,thmid,nlo,nmid);
			}
			else
			{
				xhi = xnew;
				fhi = fnew;
				thhi = thnew;
				nhi = nnew;
				SORT2(xmid,xhi,fmid,fhi,thmid,thhi,nmid,nhi);
			}
		}

		if( DEBUG_LOC )
		{
			printf(" xlo %.15e xmid %.15e xhi %.15e delta %.15e xnew %.15e\n",xlo,xmid,xhi,delta,xnew);
			printf(" flo %.15e fmid %.15e fhi %.15e\n",flo,fmid,fhi);
		}

		/* sanity check */
		/* >>chng 01 jul 26, added "fmid == 0." clause to make it run with gcc 2.96 on Linux-PC's */
		ASSERT( fmid == 0. || ( xlo <= xmid && xmid <= xhi && flo*fhi <= 0. ) );
	}

	*ThermRatio = thmid;
	gv.bin[nd]->dstpot = xmid;
	UpdatePot(nd,1);

	gv.bin[nd]->lgChrgConverged = !lgBigError;

	if( lgBigError && !conv.lgSearch )
	{
		fprintf( ioQQQ," GrainCharge did not converge ionization of grain species %s\n",
			 gv.bin[nd]->chDstLab );
		ConvFail("grai","");
	}

	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "  >Grain potential:%12.12s %.6fV, charge %.5e (e)", 
		  gv.bin[nd]->chDstLab, gv.bin[nd]->dstpot*EVRYD, POT2CHRG(gv.bin[nd]->dstpot) );
		fprintf( ioQQQ, " emis - recom: %.4e\n", fmid );
		for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
		{
			fprintf( ioQQQ, "    Thres[%ld]: %.4f V ThresVal[%ld]: %.4f V", 
				 gv.bin[nd]->chrg[nz]->DustZ, gv.bin[nd]->chrg[nz]->ThresInf*EVRYD,
				 gv.bin[nd]->chrg[nz]->DustZ, gv.bin[nd]->chrg[nz]->ThresInfVal*EVRYD );
		}
		fprintf( ioQQQ, "\n" );
	}


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

/* bracket the solution for the grain potential */
STATIC void BracketGrainPot(long nd,
			    /*@in@*/double* xlo,
			    /*@out@*/double* flo,
			    /*@out@*/double* thlo,
			    /*@out@*/double* nlo,
			    /*@in@*/double* xmid,
			    /*@out@*/double* fmid,
			    /*@out@*/double* thmid,
			    /*@out@*/double* nmid,
			    /*@in@*/double* xhi,
			    /*@out@*/double* fhi,
			    /*@out@*/double* thhi,
			    /*@out@*/double* nhi)
{
	long loop;
	double emission,
	  recombination,
	  stephi,
	  steplo;

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

	steplo = *xmid - *xlo;
	stephi = *xhi - *xmid;

	*flo = -1.;
	*thlo = 0.;
	*nlo = 0.;
	*fhi = 1.;
	*thhi = 0.;
	*nhi = 0.;

	for( loop = 0; *flo < 0. && loop < BRACKET_MAX; ++loop )
	{
		gv.bin[nd]->dstpot = *xlo;
		UpdatePot(nd,1);
		emission = GrainElecEmis(nd,thlo);
		recombination = GrainElecRecomb(nd);
		*flo = emission - recombination;
		*nlo = 0.5*(emission+recombination);
		if( *flo < 0. )
		{
			*xhi = *xlo;
			*fhi = *flo;
			*thhi = *thlo;
			*nhi = *nlo;
			*xlo = MAX2(*xlo-steplo,gv.bin[nd]->LowestPot);
			steplo *= 2.;
			/* value from previous zone is useless now */
			*xmid = (*xlo + *xhi)/2.;
		}
	}
	if( *flo == 0. )
	{
		*xmid = *xlo;
		*fmid = *flo;
		*thmid = *thlo;
		*nmid = *nlo;

#		ifdef DEBUG_FUN
		fputs( " <->BracketGrainPot()\n", debug_fp );
#		endif
		return;
	}
	else if( *flo < 0. )
	{
		/*lint -e777 float test equality */
		if( *xlo == gv.bin[nd]->LowestPot )
		/*lint +e777 float test equality */
		{
			fprintf( ioQQQ, " insanity: grain potential below lowest allowed value for %s\n" ,
				 gv.bin[nd]->chDstLab );
			ShowMe();
		}
		else
		{
			fprintf( ioQQQ, " could not bracket grain potential for %s\n" , gv.bin[nd]->chDstLab );
		}
		puts( "[Stop in BracketGrainPot]" );
		cdEXIT(EXIT_FAILURE);
	}

	for( loop = 0; *fhi > 0. && loop < BRACKET_MAX; ++loop )
	{
		gv.bin[nd]->dstpot = *xhi;
		UpdatePot(nd,1);
		emission = GrainElecEmis(nd,thhi);
		recombination = GrainElecRecomb(nd);
		*fhi = emission - recombination;
		*nhi = 0.5*(emission+recombination);
		if( *fhi > 0. )
		{
			*xlo = *xhi;
			*flo = *fhi;
			*thlo = *thhi;
			*nlo = *nhi;
			*xhi += stephi;
			stephi *= 2.;
			/* value from previous zone is useless now */
			*xmid = (*xlo + *xhi)/2.;
		}
	}
	if( *fhi == 0. )
	{
		*xmid = *xhi;
		*fmid = *fhi;
		*thmid = *thhi;
		*nmid = *nhi;

#		ifdef DEBUG_FUN
		fputs( " <->BracketGrainPot()\n", debug_fp );
#		endif
		return;
	}
	else if( *fhi > 0. )
	{
		fprintf( ioQQQ, " could not bracket grain potential for %s\n" , gv.bin[nd]->chDstLab );
		puts( "[Stop in BracketGrainPot]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* set initial estimate for grain potential */
	gv.bin[nd]->dstpot = *xmid;
	UpdatePot(nd,1);
	emission = GrainElecEmis(nd,thmid);
	recombination = GrainElecRecomb(nd);
	*fmid = emission - recombination;
	*nmid = 0.5*(emission+recombination);

	/* sanity checks */
	ASSERT( (*flo)*(*fhi) <= 0. );
	ASSERT( gv.bin[nd]->LowestPot <= *xlo && *xlo < *xhi );
	ASSERT( *xlo <= gv.bin[nd]->dstpot && gv.bin[nd]->dstpot <= *xhi );

#	ifdef DEBUG_FUN
	fputs( " <->BracketGrainPot()\n", debug_fp );
#	endif
	return;
}
#else
STATIC void GrainCharge(long int nd,
			/*@out@*/double *ThermRatio) /* ratio of thermionic to total rate */
{
	int lgBigError,
	  lgInitial;
	long backup,
	  i,
	  loopMax,
	  newZlo,
	  nz,
	  power,
	  stride,
	  stride0,
	  Zlo;
	double crate,
	  csum1,
	  csum1a,
	  csum1b,
	  csum1c,
	  csum2,
	  csum3,
	  netloss0 = -DBL_MAX,
	  netloss1 = -DBL_MAX,
	  rate_up[NCHU],
	  rate_dn[NCHU],
	  step;

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

	/* find dust charge */
	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "    Charge loop, search %c,", TorF(conv.lgSearch) );
	}

	ASSERT( nd >= 0 && nd < gv.nBin );

	for( nz=0; nz < NCHS; nz++ )
	{
		gv.bin[nd]->chrg[nz]->FracPop = -DBL_MAX;
	}

	/* this algorithm determines the value of Zlo and the charge state populations
	 * in the n-charge state model as described in:
	 *
	 * >>refer van Hoof P.A.M., Weingartner J.C., et al., 2003, MNRAS, submitted
	 *
	 * the algorithm first uses the n charge states to bracket the solution by
	 * separating the charge states with a stride that is an integral power of
	 * n-1, e.g. like this for an n == 4 calculation:
	 *
	 *  +gain    +gain    -gain    -gain
	 *    |        |        |        |
	 *   -8        1        10       19
	 *
	 * for each of the charge states the total electron emission and recombination
	 * rates are calculated. the solution is considered properly bracketed if the
	 * sign of the net gain rate (emission-recombination) is different for the first
	 * and the last of the n charge states. then the algorithm searches for two
	 * adjacent charge states where the net gain rate changes sign and divides
	 * that space into n-1 equal parts, like here
	 *
	 *   net gain  +  +  +  -
	 *             |  |  |  |
	 *        Zg   1  4  7  10
	 *
	 * note that the first and the last charge state can be retained from the
	 * previous iteration and do not need to be recalculated (UpdatePot as well
	 * as GrainElecEmis1 and GrainElecRecomb1 have mechanisms for re-using
	 * previously calculated data, so GrainCharge doesn't need to worry about
	 * this). The dividing algorithm is repeated until the stride equals 1, like
	 * here
	 *
	 *        net gain   +---
	 *                   ||||
	 *             Zg    7  10
	 *
	 * finally, the bracket may need to be shifted in order for the criterion
	 * J1 x J2 <= 0 to be fulfilled (see the paper quoted above for a detailed
	 * explanation). in the example here one shift might be necessary:
	 *
	 *        net gain  ++--
	 *                  ||||
	 *             Zg   6  9
	 *
	 * for n == 2, the algorithm would have to be slightly different. In order
	 * to avoid complicating the code, we force the code to use n == 3 in the
	 * first two steps of the process, reverting back to n == 2 upon the last
	 * step. this should not produce any noticeable additional CPU overhead */

	lgInitial = ( gv.bin[nd]->chrg[0]->DustZ == LONG_MIN );

	backup = gv.bin[nd]->nChrg;
	gv.bin[nd]->nChrg = MAX2(gv.bin[nd]->nChrg,3);

	stride0 = gv.bin[nd]->nChrg-1;

	/* set up initial bracket for grain charge, will be checked below */
	if( lgInitial )
	{
		double xxx;
		step = MAX2((double)(-gv.bin[nd]->LowestZg),1.);
		power = (int)(log(step)/log((double)stride0));
		power = MAX2(power,0);
		xxx = powi((double)stride0,power);
		stride = NINT(xxx);
		Zlo = gv.bin[nd]->LowestZg;
	}
	else
	{
		/* the previous solution is the best choice here */
		stride = 1;
		Zlo = gv.bin[nd]->chrg[0]->DustZ;
	}
	UpdatePot( nd, Zlo, stride, rate_up, rate_dn );

	/* check if the grain charge is correctly bracketed */
	for( i=0; i < BRACKET_MAX; i++ )
	{
		netloss0 = rate_up[0] - rate_dn[0];
		netloss1 = rate_up[gv.bin[nd]->nChrg-1] - rate_dn[gv.bin[nd]->nChrg-1];

		if( netloss0*netloss1 <= 0. )
			break;

		if( netloss1 > 0. )
			Zlo += (gv.bin[nd]->nChrg-1)*stride;

		if( i > 0 )
			stride *= stride0;

		if( netloss1 < 0. )
			Zlo -= (gv.bin[nd]->nChrg-1)*stride;

		Zlo = MAX2(Zlo,gv.bin[nd]->LowestZg);
		UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
	}

	if( netloss0*netloss1 > 0. ) {
		fprintf( ioQQQ, " insanity: could not bracket grain charge for %s\n", gv.bin[nd]->chDstLab );
		ShowMe();
		puts( "[Stop in GrainCharge]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* home in on the charge */
	while( stride > 1 )
	{
		stride /= stride0;

		netloss0 = rate_up[0] - rate_dn[0];
		for( nz=0; nz < gv.bin[nd]->nChrg-1; nz++ )
		{
			netloss1 = rate_up[nz+1] - rate_dn[nz+1];

			if( netloss0*netloss1 <= 0. )
			{
				Zlo = gv.bin[nd]->chrg[nz]->DustZ;
				break;
			}

			netloss0 = netloss1;
		}
		UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
	}

	ASSERT( netloss0*netloss1 <= 0. );

	gv.bin[nd]->nChrg = backup;

	/* >>chng 04 feb 15, relax upper limit on initial search when nChrg is much too large, PvH */ 
	loopMax = ( lgInitial ) ? 4*gv.bin[nd]->nChrg : 2*gv.bin[nd]->nChrg;

	lgBigError = TRUE;
	for( i=0; i < loopMax; i++ )
	{
		GetFracPop( nd, Zlo, rate_up, rate_dn, &newZlo );

		if( newZlo == Zlo )
		{
			lgBigError = FALSE;
			break;
		}

		Zlo = newZlo;
		UpdatePot( nd, Zlo, 1, rate_up, rate_dn );
	}

	if( lgBigError ) {
		fprintf( ioQQQ, " insanity: could not converge grain charge for %s\n", gv.bin[nd]->chDstLab );
		ShowMe();
		puts( "[Stop in GrainCharge]" );
		cdEXIT(EXIT_FAILURE);
	}

	/*TODO  remove gv.bin[nd]->lgChrgConverged, gv.bin[nd]->LowestPot, gv.bin[nd]->dstpotsav
	 *      gv.bin[nd]->RateUp, gv.bin[nd]->RateDn; also  gv.HighestIon??, HighestIonStage()?? */
	gv.bin[nd]->lgChrgConverged = TRUE;

	gv.bin[nd]->AveDustZ = 0.;
	crate = csum3 = 0.;
	for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
	{
		double d[5];

		gv.bin[nd]->AveDustZ += gv.bin[nd]->chrg[nz]->FracPop*gv.bin[nd]->chrg[nz]->DustZ;

		crate += gv.bin[nd]->chrg[nz]->FracPop*GrainElecEmis1(nd,nz,&d[0],&d[1],&d[2],&d[3],&d[4]);
		csum3 += gv.bin[nd]->chrg[nz]->FracPop*d[4];
	}
	gv.bin[nd]->dstpot = CHRG2POT(gv.bin[nd]->AveDustZ);
	*ThermRatio = ( crate > 0. ) ? csum3/crate : 0.;

	ASSERT( *ThermRatio >= 0. );

	if( trace.lgTrace && trace.lgDustBug )
	{
		double d[5];

		fprintf( ioQQQ, "\n" );

		crate = csum1a = csum1b = csum1c = csum2 = csum3 = 0.;
		for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
		{
			crate += gv.bin[nd]->chrg[nz]->FracPop*GrainElecEmis1(nd,nz,&d[0],&d[1],&d[2],&d[3],&d[4]);
			csum1a += gv.bin[nd]->chrg[nz]->FracPop*d[0];
			csum1b += gv.bin[nd]->chrg[nz]->FracPop*d[1];
			csum1c += gv.bin[nd]->chrg[nz]->FracPop*d[2];
			csum2 += gv.bin[nd]->chrg[nz]->FracPop*d[3];
			csum3 += gv.bin[nd]->chrg[nz]->FracPop*d[4];
		}

		fprintf( ioQQQ, "    ElecEm  rate1a=%.4e, rate1b=%.4e, rate1c=%.4e, ", csum1a, csum1b, csum1c );
		fprintf( ioQQQ, "rate2=%.4e, rate3=%.4e, sum=%.4e\n", csum2, csum3, crate );
		if( crate > 0. ) 
		{
			fprintf( ioQQQ, "      rate1a/sum=%.4e, rate1b/sum=%.4e, rate1c/sum=%.4e, ",
				 csum1a/crate, csum1b/crate, csum1c/crate );
			fprintf( ioQQQ, "rate2/sum=%.4e, rate3/sum=%.4e\n", csum2/crate, csum3/crate );
		}

		crate = csum1 = csum2 = 0.;
		for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
		{
			crate += gv.bin[nd]->chrg[nz]->FracPop*GrainElecRecomb1(nd,nz,&d[0],&d[1]);
			csum1 += gv.bin[nd]->chrg[nz]->FracPop*d[0];
			csum2 += gv.bin[nd]->chrg[nz]->FracPop*d[1];
		}

		fprintf( ioQQQ, "    ElecRc  rate1=%.4e, rate2=%.4e, sum=%.4e\n", csum1, csum2, crate );
		if( crate > 0. ) 
		{
			fprintf( ioQQQ, "      rate1/sum=%.4e, rate2/sum=%.4e\n", csum1/crate, csum2/crate );
		}

		fprintf( ioQQQ, "    Charging rates:" );
		for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
		{
			fprintf( ioQQQ, "    Zg %ld up %.4e dn %.4e",
				 gv.bin[nd]->chrg[nz]->DustZ, rate_up[nz], rate_dn[nz] );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "    FracPop: fnzone %.2f nd %ld AveZg %.5e", fnzone, nd, gv.bin[nd]->AveDustZ );
		for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
		{
			fprintf( ioQQQ, "    Zg %ld %.5f", Zlo+nz, gv.bin[nd]->chrg[nz]->FracPop );
		}
		fprintf( ioQQQ, "\n" );

		fprintf( ioQQQ, "  >Grain potential:%12.12s %.6fV", 
			 gv.bin[nd]->chDstLab, gv.bin[nd]->dstpot*EVRYD );
		for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
		{
			fprintf( ioQQQ, "    Thres[%ld]: %.4f V ThresVal[%ld]: %.4f V", 
				 gv.bin[nd]->chrg[nz]->DustZ, gv.bin[nd]->chrg[nz]->ThresInf*EVRYD,
				 gv.bin[nd]->chrg[nz]->DustZ, gv.bin[nd]->chrg[nz]->ThresInfVal*EVRYD );
		}
		fprintf( ioQQQ, "\n" );
	}

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


#ifdef OLD_POT_ALGO
STATIC double GrainElecRecomb(long int nd)
{
	long nz;
	double crate,
	  sum1[NCHS],
	  csum1,
	  sum2[NCHS],
	  csum2;

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

	/* >>chng 01 may 07, this routine now completely supports the hybrid grain
	 * charge model, and the average charge state is not used anywhere anymore, PvH */

	crate = 0.;
	csum1 = 0.;
	csum2 = 0.;

	for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
	{
		crate += gv.bin[nd]->chrg[nz]->FracPop*GrainElecRecomb1(nd,nz,&sum1[nz],&sum2[nz]);
		csum1 += gv.bin[nd]->chrg[nz]->FracPop*sum1[nz];
		csum2 += gv.bin[nd]->chrg[nz]->FracPop*sum2[nz];
	}

	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "    GrainElecRecomb finds rate1=%.4e, rate2=%.4e, sum=%.4e\n",
			 csum1, csum2, crate );
		if( crate > 0. ) 
		{
			fprintf( ioQQQ, "    rel. fractions rate1/sum=%.4e, rate2/sum=%.4e\n",
				 csum1/crate, csum2/crate );
		}
	}


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


/* grain electron recombination rates for single charge state */
STATIC double GrainElecRecomb1(long nd,
			       long nz,
			       /*@out@*/ double *sum1,
			       /*@out@*/ double *sum2)
{
	long ion,
	  nelem;
	double eta,
	  rate,
	  Stick,
	  ve,
	  xi;

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

	ASSERT( nd >= 0 && nd < gv.nBin );
	ASSERT( nz >= 0 && nz < gv.bin[nd]->nChrg );

	/* >>chng 01 may 31, try to find cached results first */
	/* >>chng 04 feb 11, moved cached data to ChargeBin, PvH */
	if( gv.bin[nd]->chrg[nz]->RSum1 >= 0. )
	{
		*sum1 = gv.bin[nd]->chrg[nz]->RSum1;
		*sum2 = gv.bin[nd]->chrg[nz]->RSum2;
		rate = *sum1 + *sum2;

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

	/* -1 makes psi correct for impact by electrons */
	ion = -1;
	/* VE is mean (not RMS) electron velocity */
	/*ve = TePowers.sqrte*6.2124e5;*/
	ve = sqrt(8.*BOLTZMANN/PI/ELECTRON_MASS*phycon.te);

	Stick = ( gv.bin[nd]->chrg[nz]->DustZ <= -1 ) ? gv.bin[nd]->StickElecNeg : gv.bin[nd]->StickElecPos;
	/* >>chng 00 jul 19, replace classical results with results including image potential
	 * to correct for polarization of the grain as charged particle approaches. */
	GrainScreen(ion,nd,nz,&eta,&xi);
	/* this is grain surface recomb rate for electrons */
	*sum1 = EPSP(gv.bin[nd]->LowestZg,gv.bin[nd]->chrg[nz]->DustZ)*Stick*dense.eden*ve*eta;

	/* >>chng 00 jul 13, add in gain rate from atoms and ions, PvH */
	*sum2 = 0.;

#ifndef IGNORE_GRAIN_ION_COLLISIONS
	for( ion=0; ion <= LIMELM; ion++ )
	{
		double CollisionRateAll = 0.;

		for( nelem=MAX2(ion-1,0); nelem < LIMELM; nelem++ )
		{
			if( dense.lgElmtOn[nelem] && dense.xIonDense[nelem][ion] > 0. &&
			    gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion] > ion )
			{
				/* this is rate with which charged ion strikes grain */
				CollisionRateAll += STICK_ION*dense.xIonDense[nelem][ion]*DoppVel.AveVel[nelem]*
					(double)(gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion]-ion);
			}
		}

		if( CollisionRateAll > 0. )
		{
			/* >>chng 00 jul 19, replace classical results with results
			 * including image potential to correct for polarization 
			 * of the grain as charged particle approaches. */
			GrainScreen(ion,nd,nz,&eta,&xi);
			*sum2 += CollisionRateAll*eta;
		}
	}
#endif

	rate = *sum1 + *sum2;

	/* >>chng 01 may 31, store results so that they may be used agian */
	gv.bin[nd]->chrg[nz]->RSum1 = *sum1;
	gv.bin[nd]->chrg[nz]->RSum2 = *sum2;

	ASSERT( *sum1 >= 0. && *sum2 >= 0. );

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


#ifdef OLD_POT_ALGO
/* GrainElecEmis computes electron loss rate for this grain species due
 * to photoelectric effect, positive ion recombinations on the
 * grain surface and thermionic emissions */
STATIC double GrainElecEmis(long int nd,
			    double *ThermRatio) /* ratio of thermionic to total emissions */
{
	long int nz;
	double crate,
	  sum1a[NCHS],
	  csum1a,
	  sum1b[NCHS],
	  csum1b,
	  sum1c[NCHS],
	  csum1c,
	  sum2[NCHS],
	  csum2,
	  sum3[NCHS],
	  csum3;

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

	/* >>chng 01 may 07, this routine now completely supports the hybrid grain
	 * charge model, and the average charge state is not used anywhere anymore, PvH */

	/* this is the loss rate due to photo-electric effect */
	crate = 0.;
	csum1a = 0.;
	csum1b = 0.;
	csum1c = 0.;
	csum2 = 0.;
	csum3 = 0.;
	for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
	{
		crate += gv.bin[nd]->chrg[nz]->FracPop*
			GrainElecEmis1(nd,nz,&sum1a[nz],&sum1b[nz],&sum1c[nz],&sum2[nz],&sum3[nz]);
		csum1a += gv.bin[nd]->chrg[nz]->FracPop*sum1a[nz];
		csum1b += gv.bin[nd]->chrg[nz]->FracPop*sum1b[nz];
		csum1c += gv.bin[nd]->chrg[nz]->FracPop*sum1c[nz];
		csum2 += gv.bin[nd]->chrg[nz]->FracPop*sum2[nz];
		csum3 += gv.bin[nd]->chrg[nz]->FracPop*sum3[nz];
	}

	*ThermRatio = ( crate > 0. ) ? csum3/crate : 0.;
		
	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "    GrainElecEmis finds rate1a=%.4e, rate1b=%.4e, ", csum1a, csum1b );
		fprintf( ioQQQ, "rate1c=%.4e, rate2=%.4e, rate3=%.4e, sum=%.4e\n", csum1c, csum2, csum3, crate );
		if( crate > 0. ) 
		{
			fprintf( ioQQQ, "    rate1a/sum=%.4e, rate1b/sum=%.4e, rate1c/sum=%.4e, ",
				 csum1a/crate, csum1b/crate, csum1c/crate );
			fprintf( ioQQQ, "rate2/sum=%.4e, rate3/sum=%.4e\n",
				 csum2/crate, csum3/crate );
		}
	}

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


/* grain electron emission rates for single charge state */
STATIC double GrainElecEmis1(long nd,
			     long nz,
			     /*@out@*/ double *sum1a,
			     /*@out@*/ double *sum1b,
			     /*@out@*/ double *sum1c,
			     /*@out@*/ double *sum2,
			     /*@out@*/ double *sum3)
{
	long int i,
	  ion,
	  ipLo,
	  nelem;
	double eta,
	  rate,
	  xi;

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

	ASSERT( nd >= 0 && nd < gv.nBin );
	ASSERT( nz >= 0 && nz < gv.bin[nd]->nChrg );

	/* >>chng 01 may 31, try to find cached results first */
	/* >>chng 04 feb 11, moved cached data to ChargeBin, PvH */
	if( gv.bin[nd]->chrg[nz]->ESum1a >= 0. )
	{
		*sum1a = gv.bin[nd]->chrg[nz]->ESum1a;
		*sum1b = gv.bin[nd]->chrg[nz]->ESum1b;
		*sum1c = gv.bin[nd]->chrg[nz]->ESum1c;
		*sum2 = gv.bin[nd]->chrg[nz]->ESum2;
		/* don't cache thermionic rates as they depend on grain temp */
		*sum3 = 4.*gv.bin[nd]->chrg[nz]->ThermRate;
		rate = *sum1a + *sum1b + *sum1c + *sum2 + *sum3;

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

	/* this is the loss rate due to photo-electric effect */
	*sum1a = 0.;
	ipLo = gv.bin[nd]->chrg[nz]->ipThresInfVal;
	for( i=ipLo; i < rfield.nflux; i++ )
	{
#		ifdef WD_TEST2
		*sum1a += rfield.flux[i]*gv.bin[nd]->dstab1[i]*gv.bin[nd]->chrg[nz]->yhat[i-ipLo];
#		else
		*sum1a += rfield.SummedCon[i]*gv.bin[nd]->dstab1[i]*gv.bin[nd]->chrg[nz]->yhat[i-ipLo];
#		endif
	}
	/* normalize to rates per cm^2 of projected grain area */
	*sum1a /= gv.bin[nd]->IntArea/4.;

	/* >>chng 01 dec 18, added code for modeling secondary electron emissions, PvH */
	/* this code does a crude correction for the Auger effect in grains,
	 * it is roughly valid for neutral and negative grains, but overestimates
	 * the effect for positively charged grains. Many of the Auger electrons have
	 * rather low energies and will not make it out of the potential well for
	 * high grain potentials typical of AGN conditions, see Table 4.1 & 4.2 of
	 * >>refer	grain	physics	Dwek E. & Smith R.K., 1996, ApJ, 459, 686 */

	/*TODO	2	note that the number of primary electrons is given by yhat,
	 * which may not be one, so this is not necessarily consistent */
	/*TODO	2	avAuger depends on grain charge, this should be treated explicitly here */

	*sum1b = 0.; /* this will be the secondary electron rate, in e/s/H at default depl */
#	if 0
	broken(); /* the Auger/Compton-recoil code is severely flawed, see hotdust.in case */

	ion = 0; /* do Auger effect only for neutral atoms */

	/* this can be turned off with the NO AUGER command, since that command sets
	 * the array yield.n_elec_eject to zero */

	/* no Auger electrons for first three elements, but there are compton
	 * recoil ionizations for those elements, so start with Hydrogen */
	for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
	{
		long ns, nelec, nej;

		if( gv.bin[nd]->elmAbund[nelem] > 0. )
		{
			double eject_rate;
			/* loop over all shells, except the valence shell which doesn't produce Auger electrons
			 * and is properly part of the low energy grain opacity */
			for( ns=0; ns < Heavy.nsShells[nelem][ion]-1; ns++ )
			{
				/* this will become the average number of Auger electrons freed */
				double avAuger = 0.;
				/* following is most number of electrons freed */
				nelec = yield.n_elec_eject[nelem][ion][ns];

				/* following loop is over number of electrons that come out of shell
				 * with multiple electron ejection - use (nej-1) so that primary
				 * photoelectron is not counted, which is already done above. */
				for( nej=2; nej <= nelec; nej++ )
				{
					avAuger += yield.frac_elec_eject[nelem][ion][ns][nej-1]*(double)(nej-1);
				}
				/* add on Compton recoil of K electrons, with yield */
				if( ns==0 )
				{
					/* first term is Auger shakeoff electrons */
					eject_rate = ionbal.PhotoRate_Shell[nelem][ion][ns][0]*avAuger +
						/* this is compton recoil ionization term, which adds
						 * primary plus Auger electrons */
						/* >>chng 02 jan 11, add compton recoil ionization */
						/* >>chng 02 mar 24, recoil rate was added to outer
						 * shell of atom in photoionize, but not to inner shell rate */
						/* add ionization due to compton recoil - 
						 * this is turned off with the NO RECOIL command */
						ionbal.CompRecoilIonRate[nelem][ion]*(avAuger+1.);
				}
				else
				{
					/* not K shell, only include explicit shell */
					eject_rate = ionbal.PhotoRate_Shell[nelem][ion][ns][0]*avAuger;
				}
				/* ionbal.PhotoRate_Shell is the photo-ionization rate per atom per second,
				 * it is updated after grain() is called but before RT_OTS_Update */
				*sum1b  += eject_rate * gv.bin[nd]->elmAbund[nelem];
			}
		}
	}
#	endif

#	ifdef WD_TEST2
	*sum1b = 0.;
#	else
	*sum1b /= gv.bin[nd]->IntArea/4.;
#	endif

	*sum1c = 0.;
	if( gv.bin[nd]->chrg[nz]->DustZ <= -1 )
	{
		ipLo = gv.bin[nd]->chrg[nz]->ipThresInf;
		for( i=ipLo; i < rfield.nflux; i++ )
		{
			/* >>chng 00 jul 17, use description of Weingartner & Draine, 2001 */
#			ifdef WD_TEST2
			*sum1c += rfield.flux[i]*gv.bin[nd]->chrg[nz]->cs_pdt[i-ipLo];
#			else
			*sum1c += rfield.SummedCon[i]*gv.bin[nd]->chrg[nz]->cs_pdt[i-ipLo];
#			endif
		}
		*sum1c /= gv.bin[nd]->IntArea/4.;
	}

	/* >>chng 00 jun 19, add in loss rate due to recombinations with ions, PvH */
	*sum2 = 0.;
#	ifndef IGNORE_GRAIN_ION_COLLISIONS
	for( ion=0; ion <= LIMELM; ion++ )
	{
		double CollisionRateAll = 0.;

		for( nelem=MAX2(ion-1,0); nelem < LIMELM; nelem++ )
		{
			if( dense.lgElmtOn[nelem] && dense.xIonDense[nelem][ion] > 0. &&
			    ion > gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion] )
			{
				/* this is rate with which charged ion strikes grain */
				CollisionRateAll += STICK_ION*dense.xIonDense[nelem][ion]*DoppVel.AveVel[nelem]*
					(double)(ion-gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion]);
			}
		}

		if( CollisionRateAll > 0. )
		{
			/* >>chng 00 jul 19, replace classical results with results
			 * including image potential to correct for polarization 
			 * of the grain as charged particle approaches. */
			GrainScreen(ion,nd,nz,&eta,&xi);
			*sum2 += CollisionRateAll*eta;
		}
	}
#	endif

	/* >>chng 01 may 30, moved calculation of ThermRate to UpdatePot */
	/* >>chng 01 jan 19, multiply by 4 since thermionic emissions scale with total grain
	 * surface area while the above two processes scale with projected grain surface area, PvH */
	*sum3 = 4.*gv.bin[nd]->chrg[nz]->ThermRate;

	rate = *sum1a + *sum1b + *sum1c + *sum2 + *sum3;

	/* >>chng 01 may 31, store results so that they may be used agian */
	gv.bin[nd]->chrg[nz]->ESum1a = *sum1a;
	gv.bin[nd]->chrg[nz]->ESum1b = *sum1b;
	gv.bin[nd]->chrg[nz]->ESum1c = *sum1c;
	gv.bin[nd]->chrg[nz]->ESum2 = *sum2;

	ASSERT( *sum1a >= 0. && *sum1b >= 0. && *sum1c >= 0. && *sum2 >= 0. && *sum3 >= 0. );

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


/* correction factors for grain charge screening (including image potential
 * to correct for polarization of the grain as charged particle approaches). */
STATIC void GrainScreen(long ion,
			long nd,
			long nz,
			/*@out@*/ double *eta,
			/*@out@*/ double *xi)
{
#	ifdef DEBUG_FUN
	fputs( "<+>GrainScreen()\n", debug_fp );
#	endif

	/* >>chng 04 jan 31, start caching eta, xi results, PvH */

	/* add 1 to allow for electron charge ion = -1 */
	long ind = ion+1;
	ASSERT( ind >= 0 && ind < LIMELM+2 );

	if( gv.bin[nd]->chrg[nz]->eta[ind] > 0. )
	{
		*eta = gv.bin[nd]->chrg[nz]->eta[ind];
		*xi = gv.bin[nd]->chrg[nz]->xi[ind];

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

	/* >>refer	grain	physics	Draine & Sutin, 1987, ApJ, 320, 803
	 * eta = J-tilde (eq. 3.3 thru 3.5), xi = Lambda-tilde/2. (eq. 3.8 thru 3.10) */
	if( ion == 0 ) 
	{
		*eta = 1.;
		*xi = 1.;
	}
	else 
	{
		/* >>chng 01 jan 03, assume that grain charge is distributed in two states just below and
		 * above the average charge, instead of the delta function we assume elsewhere. by averaging
		 * over the distribution we smooth out the discontinuities of the the Draine & Sutin expressions
		 * around nu == 0. they were the cause of temperature instabilities in globule.in. PvH */
		/* >>chng 01 may 07, revert back to single charge state since only integral charge states are
		 * fed into this routine now, making the two-charge state approximation obsolete, PvH */
		double nu = (double)gv.bin[nd]->chrg[nz]->DustZ/(double)ion;
		double tau = gv.bin[nd]->Capacity*BOLTZMANN*phycon.te*1.e-7/POW2((double)ion*ELEM_CHARGE);
		if( nu < 0. )
		{
			*eta = (1. - nu/tau)*(1. + sqrt(2./(tau - 2.*nu)));
			*xi = (1. - nu/(2.*tau))*(1. + 1./sqrt(tau - nu));
		}
		else if( nu == 0. ) 
		{
			*eta = 1. + sqrt(PI/(2.*tau));
			*xi = 1. + 0.75*sqrt(PI/(2.*tau));
		}
		else 
		{
			double theta_nu = ThetaNu(nu);
			/* >>>chng 00 jul 27, avoid passing functions to macro so set to local temp var */
			double xxx = 1. + 1./sqrt(4.*tau+3.*nu);
			*eta = POW2(xxx)*exp(-theta_nu/tau);
#			ifdef WD_TEST2
			*xi = (1. + nu/(2.*tau))*(1. + 1./sqrt(3./(2.*tau)+3.*nu))*exp(-theta_nu/tau);
#			else
			/* >>chng 01 jan 24, use new expression for xi which only contains the excess
			 * energy above the potential barrier of the incoming particle (accurate to
			 * 2% or better), and then add in potential barrier separately, PvH */
			xxx = 0.25*pow(nu/tau,0.75)/(pow(nu/tau,0.75) + pow((25.+3.*nu)/5.,0.75)) +
				(1. + 0.75*sqrt(PI/(2.*tau)))/(1. + sqrt(PI/(2.*tau)));
			*xi = (MIN2(xxx,1.) + theta_nu/(2.*tau))*(*eta);
#			endif
		}
		
		ASSERT( *eta >= 0. && *xi >= 0. );
	}

	gv.bin[nd]->chrg[nz]->eta[ind] = *eta;
	gv.bin[nd]->chrg[nz]->xi[ind] = *xi;

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


STATIC double ThetaNu(double nu)
{
	double theta_nu;

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

	if( nu > 0. )
	{
		double xxx;
		const double REL_TOLER = 10.*DBL_EPSILON;

		/* >>chng 01 jan 24, get first estimate for xi_nu and iteratively refine, PvH */
		double xi_nu = 1. + 1./sqrt(3.*nu);
		double xi_nu2 = POW2(xi_nu);
		do 
		{
			double old = xi_nu;
			/* >>chng 04 feb 01, use Newton-Raphson to speed up convergence, PvH */
                        /* xi_nu = sqrt(1.+sqrt((2.*POW2(xi_nu)-1.)/(nu*xi_nu))); */
			double fnu = 2.*xi_nu2 - 1. - nu*xi_nu*POW2(xi_nu2 - 1.);
			double dfdxi = 4.*xi_nu - nu*((5.*xi_nu2 - 6.)*xi_nu2 + 1.);
			xi_nu -= fnu/dfdxi;
			xi_nu2 = POW2(xi_nu);
			xxx = fabs(old-xi_nu);
		} while( xxx > REL_TOLER*xi_nu );

		theta_nu = nu/xi_nu - 1./(2.*xi_nu2*(xi_nu2-1.));
	}
	else
	{
		theta_nu = 0.;
	}

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


#ifdef OLD_POT_ALGO
/* update items that depend on grain potential */
STATIC void UpdatePot(long nd,
		      long stride)
{
	long loop,
	  newZlo,
	  nz,
	  Zg,
	  Zlo;
	double BoltzFac,
	  HighEnergy;
	const long MAXLOOP = 10L;

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

	/* sanity checks */
	ASSERT( gv.bin[nd]->dstpot >= gv.bin[nd]->LowestPot );
	/* temporary */
	ASSERT( stride == 1 );

	/* MAX2 is to protect against roundoff error when dstpot == LowestPot */
	gv.bin[nd]->AveDustZ = MAX2(POT2CHRG(gv.bin[nd]->dstpot),(double)gv.bin[nd]->LowestZg);

	/* invalidate level populations */
	for( nz=0; nz < NCHS; nz++ )
	{
		gv.bin[nd]->chrg[nz]->FracPop = -DBL_MAX;
	}

	/* >>chng 01 may 30, calculation of newZlo is now valid in n-charge state model, unchanged for n = 2 */
	if( gv.bin[nd]->nChrg%2 == 0 )
	{
		double help = gv.bin[nd]->AveDustZ - (double)stride/2.;
		/* integer division is intended in next two statements */
		newZlo = NINT(help) - ((gv.bin[nd]->nChrg-2)*stride)/2;
	}
	else
	{
		newZlo = NINT(gv.bin[nd]->AveDustZ) - ((gv.bin[nd]->nChrg-1)*stride)/2;
	}
	newZlo = MAX2(newZlo,gv.bin[nd]->LowestZg);

	loop = 0;

	do
	{
		Zlo = newZlo;

		/* >>chng 00 jul 17, use description of Weingartner & Draine, 2001 */
		/* >>chng 01 mar 21, assume that grain charge is distributed in two states just below and
		 * above the average charge. */
		/* >>chng 01 may 07, this routine now completely supports the hybrid grain
		 * charge model, and the average charge state is not used anywhere anymore, PvH */
		/* >>chng 01 may 30, reorganize code such that all relevant data can be copied
		 * when a valid set of data is available from a previous call, this saves CPU time, PvH */
		/* >>chng 04 jan 17, reorganized code to use pointers to the charge bins, PvH */

		for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
		{
			long ind, zz;
			ChargeBin *ptr;

			Zg = Zlo + nz*stride;

			/* check if charge state is already present */
			ind = NCHS-1;
			for( zz=0; zz < NCHS-1; zz++ )
			{
				if( gv.bin[nd]->chrg[zz]->DustZ == Zg )
				{
					ind = zz;
					break;
				}
			}

			/* in the code below the old charge bins are shifted to the back in order to assure
			 * that the most recently used ones are upfront; they are more likely to be reused */
			ptr = gv.bin[nd]->chrg[ind];

			/* make room to swap in charge state */
			for( zz=ind-1; zz >= nz; zz-- )
				gv.bin[nd]->chrg[zz+1] = gv.bin[nd]->chrg[zz];

			gv.bin[nd]->chrg[nz] = ptr;

			if( gv.bin[nd]->chrg[nz]->DustZ != Zg )
				UpdatePot1(nd,nz,Zg,0);
			else if( gv.bin[nd]->chrg[nz]->nfill < rfield.nflux )
				UpdatePot1(nd,nz,Zg,gv.bin[nd]->chrg[nz]->nfill);

			UpdatePot2(nd,nz);

			ASSERT( gv.bin[nd]->chrg[nz]->DustZ == Zg );
			ASSERT( gv.bin[nd]->chrg[nz]->nfill >= rfield.nflux );
		}

		GetFracPop(nd,gv.bin[nd]->AveDustZ,Zlo,&newZlo);

		++loop;
	}
	while( loop < MAXLOOP && newZlo != Zlo );

	if( newZlo != Zlo )
	{
		fprintf( ioQQQ, " PROBLEM  UpdatePot did not converge level populations\n" );
		puts( "[Stop in UpdatePot]" );
		cdEXIT(EXIT_FAILURE);
	}
		
	/* determine highest energy to be considered by quantum heating routines.
	 * since the Boltzmann distribution is resolved, the upper limit has to be
	 * high enough that a negligible amount of energy is in the omitted tail */
	/* >>chng 03 jan 26, moved this code from GrainChargeTemp to UpdatePot
	 *                   since the new code depends on grain potential, HTT91.in, PvH */
	BoltzFac = (-log(CONSERV_TOL) + 8.)*BOLTZMANN/EN1RYD;
	HighEnergy = 0.;
	for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
	{
		/* >>chng 04 jan 21, changed phycon.te -> MAX2(phycon.te,gv.bin[nd]->tedust), PvH */
		HighEnergy = MAX2(HighEnergy,
		  MAX2(gv.bin[nd]->chrg[nz]->ThresInfInc,0.) + BoltzFac*MAX2(phycon.te,gv.bin[nd]->tedust));
	}
	HighEnergy = MIN2(HighEnergy,rfield.anu[rfield.nupper-1]);
	gv.bin[nd]->qnflux2 = ipoint(HighEnergy);
	gv.bin[nd]->qnflux = MAX2(rfield.nflux,gv.bin[nd]->qnflux2);

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


/* calculate charge state populations */
STATIC void GetFracPop(long nd,
		       double AveZg,
		       long Zlo,
		       /*@out@*/ long *newZlo)
{
	long nz;

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

	if( gv.bin[nd]->nChrg == 2 )
	{
		for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
		{
			long Zg = Zlo + nz;
			gv.bin[nd]->chrg[nz]->FracPop = 1.-fabs(AveZg-(double)Zg);
		}
		*newZlo = Zlo;
	}
	else
	{
		long i,j,k;
		double avePop[2],d[5],frac,pop[2][NCHS-1],rate_up[NCHS],rate_dn[NCHS];
#		if 0
		double test1,test2;
#		endif
		/* this will stop the lint from generating a warning */
		for( nz=0; nz < NCHS; ++nz )
		{
			rate_up[nz] = -FLT_MAX;
			rate_dn[nz] = -FLT_MAX;
		}

		/* initialize charge transfer rates */
		for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
		{
			rate_up[nz] = GrainElecEmis1(nd,nz,&d[0],&d[1],&d[2],&d[3],&d[4]);
			rate_dn[nz] = GrainElecRecomb1(nd,nz,&d[0],&d[1]);
			/* sanity check */
			ASSERT( rate_up[nz] >= 0. && rate_dn[nz] >= 0. );
		}

		/* solve level populations for levels 0..nChrg-2 (i == 0) and
		 * levels 1..nChrg-1 (i == 1), and determine average charge
		 * for each of those subsystems. Next we demand that
		 *         avePop[0] <= AveZg <= avePop[1]
		 * and determine FracPop by linearly adding the subsystems such that
		 *      AveZg == frac*avePop[0] + (1-frac)*avePop[1]
		 * this assures that all charge state populations are positive and all
		 * emission and recombination rates behave continuously as AveZg changes */
		for( i=0; i < 2; i++ )
		{
			double sum;

			pop[i][0] = 1.;
			sum = pop[i][0];
			for( j=1; j < gv.bin[nd]->nChrg-1; j++ )
			{
				nz = i + j;
				if( rate_dn[nz] > 10.*rate_up[nz-1]/sqrt(DBL_MAX) )
				{
					pop[i][j] = pop[i][j-1]*rate_up[nz-1]/rate_dn[nz];
					sum += pop[i][j];
				}
				else
				{
					for( k=0; k < j; k++ )
					{
						pop[i][k] = 0.;
					}
					pop[i][j] = 1.;
					sum = pop[i][j];
				}
				/* guard against overflow */
				if( pop[i][j] > sqrt(DBL_MAX) )
				{
					for( k=0; k <= j; k++ )
					{
						pop[i][k] /= DBL_MAX/10.;
					}
					sum /= DBL_MAX/10.;
				}
			}
			avePop[i] = 0.;
			for( j=0; j < gv.bin[nd]->nChrg-1; j++ )
			{
				nz = i + j;
				pop[i][j] /= sum;
				avePop[i] += pop[i][j]*(Zlo + nz);
			}
		}

		/* ascertain that the choice of Zlo was correct, this is to ensure positive
		 * level populations and continuous emission and recombination rates */
		if( AveZg < avePop[0] )
		{
			*newZlo = MAX2(Zlo - 1,gv.bin[nd]->LowestZg);
		}
		else if( AveZg > avePop[1] )
		{
			*newZlo = Zlo + 1;
		}
		else
		{
			*newZlo = Zlo;
		}

		frac = (AveZg - avePop[1])/(avePop[0] - avePop[1]);
		gv.bin[nd]->chrg[0]->FracPop = frac*pop[0][0];
		gv.bin[nd]->chrg[gv.bin[nd]->nChrg-1]->FracPop = (1.-frac)*pop[1][gv.bin[nd]->nChrg-2];
		for( nz=1; nz < gv.bin[nd]->nChrg-1; nz++ )
		{
			gv.bin[nd]->chrg[nz]->FracPop = frac*pop[0][nz] + (1.-frac)*pop[1][nz-1];
		}

#		if 0
		test1 = test2 = 0.;
		for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
		{
			test1 += gv.bin[nd]->chrg[nz]->FracPop;
			test2 += gv.bin[nd]->chrg[nz]->FracPop*(double)(Zlo+nz);
		}
		printf(" GetFracPop test 1: %.15e\n",test1);
		printf(" GetFracPop test 2: %.15e %.15e\n",test2,AveZg);
#		endif
	}

	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "      FracPop: fnzone %.2f nd %ld AveZg %.4e", fnzone, nd, AveZg );
		for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
		{
			fprintf( ioQQQ, " Zg %ld %.4e", Zlo+nz, gv.bin[nd]->chrg[nz]->FracPop );
		}
		fprintf( ioQQQ, " newZlo %ld\n", *newZlo );
	}

#	ifdef DEBUG_FUN
	fputs( " <->GetFracPop()\n", debug_fp );
#	endif
	return;
}
#else
/* update items that depend on grain potential */
STATIC void UpdatePot(long nd,
		      long Zlo,
		      long stride,
		      /*@out@*/ double rate_up[], /* rate_up[NCHU] */
		      /*@out@*/ double rate_dn[]) /* rate_dn[NCHU] */
{
	long nz,
	  Zg;
	double BoltzFac,
	  HighEnergy;

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

	ASSERT( nd >= 0 && nd < gv.nBin );
	ASSERT( Zlo >= gv.bin[nd]->LowestZg );
	ASSERT( stride >= 1 );

	/* >>chng 00 jul 17, use description of Weingartner & Draine, 2001 */
	/* >>chng 01 mar 21, assume that grain charge is distributed in two states just below and
	 * above the average charge. */
	/* >>chng 01 may 07, this routine now completely supports the hybrid grain
	 * charge model, and the average charge state is not used anywhere anymore, PvH */
	/* >>chng 01 may 30, reorganize code such that all relevant data can be copied
	 * when a valid set of data is available from a previous call, this saves CPU time, PvH */
	/* >>chng 04 jan 17, reorganized code to use pointers to the charge bins, PvH */

	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, " %ld/%ld", Zlo, stride );
	}

	for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
	{
		long ind, zz;
		double d[5];
		ChargeBin *ptr;

		Zg = Zlo + nz*stride;

		/* check if charge state is already present */
		ind = NCHS-1;
		for( zz=0; zz < NCHS-1; zz++ )
		{
			if( gv.bin[nd]->chrg[zz]->DustZ == Zg )
			{
				ind = zz;
				break;
			}
		}

		/* in the code below the old charge bins are shifted to the back in order to assure
		 * that the most recently used ones are upfront; they are more likely to be reused */
		ptr = gv.bin[nd]->chrg[ind];

		/* make room to swap in charge state */
		for( zz=ind-1; zz >= nz; zz-- )
			gv.bin[nd]->chrg[zz+1] = gv.bin[nd]->chrg[zz];

		gv.bin[nd]->chrg[nz] = ptr;

		if( gv.bin[nd]->chrg[nz]->DustZ != Zg )
			UpdatePot1(nd,nz,Zg,0);
		else if( gv.bin[nd]->chrg[nz]->nfill < rfield.nflux )
			UpdatePot1(nd,nz,Zg,gv.bin[nd]->chrg[nz]->nfill);

		UpdatePot2(nd,nz);

		rate_up[nz] = GrainElecEmis1(nd,nz,&d[0],&d[1],&d[2],&d[3],&d[4]);
		rate_dn[nz] = GrainElecRecomb1(nd,nz,&d[0],&d[1]);

		/* sanity checks */
		ASSERT( gv.bin[nd]->chrg[nz]->DustZ == Zg );
		ASSERT( gv.bin[nd]->chrg[nz]->nfill >= rfield.nflux );
		ASSERT( rate_up[nz] >= 0. && rate_dn[nz] >= 0. );
	}
		
	/* determine highest energy to be considered by quantum heating routines.
	 * since the Boltzmann distribution is resolved, the upper limit has to be
	 * high enough that a negligible amount of energy is in the omitted tail */
	/* >>chng 03 jan 26, moved this code from GrainChargeTemp to UpdatePot
	 *                   since the new code depends on grain potential, HTT91.in, PvH */
	BoltzFac = (-log(CONSERV_TOL) + 8.)*BOLTZMANN/EN1RYD;
	HighEnergy = 0.;
	for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
	{
		/* >>chng 04 jan 21, changed phycon.te -> MAX2(phycon.te,gv.bin[nd]->tedust), PvH */
		HighEnergy = MAX2(HighEnergy,
		  MAX2(gv.bin[nd]->chrg[nz]->ThresInfInc,0.) + BoltzFac*MAX2(phycon.te,gv.bin[nd]->tedust));
	}
	HighEnergy = MIN2(HighEnergy,rfield.anu[rfield.nupper-1]);
	gv.bin[nd]->qnflux2 = ipoint(HighEnergy);
	gv.bin[nd]->qnflux = MAX2(rfield.nflux,gv.bin[nd]->qnflux2);

	ASSERT( gv.bin[nd]->qnflux <= rfield.nupper );

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


/* calculate charge state populations */
STATIC void GetFracPop(long nd,
		       long Zlo,
		       /*@in@*/ double rate_up[], /* rate_up[NCHU] */
		       /*@in@*/ double rate_dn[], /* rate_dn[NCHU] */
		       /*@out@*/ long *newZlo)
{
	int lgRedo;
	long i,
	  nz;
	double netloss[2],
	  pop[2][NCHU-1];
		

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

	ASSERT( nd >= 0 && nd < gv.nBin );
	ASSERT( Zlo >= gv.bin[nd]->LowestZg );

	/* solve level populations for levels 0..nChrg-2 (i == 0) and
	 * levels 1..nChrg-1 (i == 1), and determine net electron loss rate
	 * for each of those subsystems. Next we demand that
	 *          netloss[1] <= 0 <= netloss[0]
	 * and determine FracPop by linearly adding the subsystems such that
	 *     0 == frac*netloss[0] + (1-frac)*netloss[1]
	 * this assures that all charge state populations are positive */
	do
	{
		for( i=0; i < 2; i++ )
		{
			long j, k;
			double sum;

			sum = pop[i][0] = 1.;
			for( j=1; j < gv.bin[nd]->nChrg-1; j++ )
			{
				nz = i + j;
				if( rate_dn[nz] > 10.*rate_up[nz-1]/sqrt(DBL_MAX) )
				{
					pop[i][j] = pop[i][j-1]*rate_up[nz-1]/rate_dn[nz];
					sum += pop[i][j];
				}
				else
				{
					for( k=0; k < j; k++ )
					{
						pop[i][k] = 0.;
					}
					pop[i][j] = 1.;
					sum = pop[i][j];
				}
				/* guard against overflow */
				if( pop[i][j] > sqrt(DBL_MAX) )
				{
					for( k=0; k <= j; k++ )
					{
						pop[i][k] /= DBL_MAX/10.;
					}
					sum /= DBL_MAX/10.;
				}
			}
			netloss[i] = 0.;
			for( j=0; j < gv.bin[nd]->nChrg-1; j++ )
			{
				nz = i + j;
				pop[i][j] /= sum;
				netloss[i] += pop[i][j]*(rate_up[nz] - rate_dn[nz]);
			}
		}

		/* ascertain that the choice of Zlo was correct, this is to ensure positive
		 * level populations and continuous emission and recombination rates */
		if( netloss[0]*netloss[1] > 0. )
			*newZlo = ( netloss[1] > 0. ) ? Zlo + 1 : Zlo - 1;
		else
			*newZlo = Zlo;

		/* >>chng 04 feb 15, add protection for roundoff error when nChrg is much too large;
		 * netloss[0/1] can be almost zero, but may have arbitrary sign due to roundoff error;
		 * this can lead to a spurious lowering of Zlo below LowestZg, orion_pdr10.in,
		 * since newZlo cannot be lowered, lower nChrg instead and recalculate, PvH */
		/* >>chng 04 feb 15, also lower nChrg if population of highest charge state is marginal */
		if( gv.bin[nd]->nChrg > 2 &&
		    ( *newZlo < gv.bin[nd]->LowestZg ||
		    ( *newZlo == Zlo && pop[1][gv.bin[nd]->nChrg-2] < DBL_EPSILON ) ) )
		{
			gv.bin[nd]->nChrg--;
			lgRedo = TRUE;
		}
		else
		{
			lgRedo = FALSE;
		}

#		if 0
		printf( " fnzone %.2f nd %ld Zlo %ld newZlo %ld netloss %.4e %.4e nChrg %ld lgRedo %d\n",
			fnzone, nd, Zlo, *newZlo, netloss[0], netloss[1], gv.bin[nd]->nChrg, lgRedo );
#		endif
	}
	while( lgRedo );

	if( *newZlo < gv.bin[nd]->LowestZg )
	{
		fprintf( ioQQQ, " could not converge charge state populations for %s\n", gv.bin[nd]->chDstLab );
		ShowMe();
		puts( "[Stop in GetFracPop]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( *newZlo == Zlo )
	{
		double frac0, frac1;
#		ifndef NDEBUG
		double test1, test2, test3, x1, x2;
#		endif

		/* frac1 == 1-frac0, but calculating it like this avoids cancellation errors */
		frac0 = netloss[1]/(netloss[1]-netloss[0]);
		frac1 = -netloss[0]/(netloss[1]-netloss[0]);

		gv.bin[nd]->chrg[0]->FracPop = frac0*pop[0][0];
		gv.bin[nd]->chrg[gv.bin[nd]->nChrg-1]->FracPop = frac1*pop[1][gv.bin[nd]->nChrg-2];
		for( nz=1; nz < gv.bin[nd]->nChrg-1; nz++ )
		{
			gv.bin[nd]->chrg[nz]->FracPop = frac0*pop[0][nz] + frac1*pop[1][nz-1];
		}

#		ifndef NDEBUG
		test1 = test2 = test3 = 0.;
		for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
		{
			ASSERT( gv.bin[nd]->chrg[nz]->FracPop >= 0. );
			test1 += gv.bin[nd]->chrg[nz]->FracPop;
			test2 += gv.bin[nd]->chrg[nz]->FracPop*rate_up[nz];
			test3 += gv.bin[nd]->chrg[nz]->FracPop*(rate_up[nz]-rate_dn[nz]);
		}
		x1 = fabs(test1-1.);
		x2 = fabs(test3/test2);
		ASSERT( MAX2(x1,x2) < 10.*sqrt(gv.bin[nd]->nChrg)*DBL_EPSILON );
#		endif
	}

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


/* this routine updates all quantities that depend only on grain charge and radius
 *
 * NB NB - All global data in grain.c and grainvar.h that are charge dependent should
 *         be calculated here or in UpdatePot2 (except gv.bin[nd]->chrg[nz]->FracPop
 *         which is special).
 *
 * NB NB - the code assumes that the data calculated here remain valid throughout
 *         the model, i.e. do NOT depend on grain temperature, etc.
 */
STATIC void UpdatePot1(long nd,
		       long nz,
		       long Zg,
		       long ipStart)
{
	long i,
	  ipLo,
	  ipHi,
	  nfill;
	size_t alloc_size;
	pe_type pcase;
	double c1,
	  cnv_GR_pH,
	  d[2],
	  DustWorkFcn,
	  Elo,
	  PotSurf,
	  ThresInf,
	  ThresInfVal,
	  ThresSurfVal;

	float *yone = gv.bin[nd]->y1;
	float *anu = rfield.anu;
	float *yhat;
	double *cs_pdt, *fac1, *fac2;

	/* constants in the expression for the photodetachment cross section
	 * taken from 
	 * >>refer	grain	physics	Weingartner & Draine, ApJS, 2001, 134, 263 */
	const double INV_DELTA_E = EVRYD/3.;
	const double CS_PDT = 1.2e-17;

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

	/* >>chng 04 jan 23, replaced rfield.nflux with rfield.nupper so
	 *                   that data remains valid throughout the model, PvH */
	/* >>chng 04 jan 26, start using ipStart to solve the validity problem 
	 *      ipStart == 0 means that a full initialization needs to be done
	 *      ipStart > 0  means that rfield.nflux was updated and yhat, cs_pdt,
	 *                   fac1, and fac2 need to be reallocated, PvH */
	if( ipStart == 0 )
	{
		gv.bin[nd]->chrg[nz]->DustZ = Zg;

		/* invalidate eta and xi storage */
		memset( gv.bin[nd]->chrg[nz]->eta, 0, (LIMELM+2)*sizeof(double) );
		memset( gv.bin[nd]->chrg[nz]->xi, 0, (LIMELM+2)*sizeof(double) );

		GetPotValues(nd,Zg,&gv.bin[nd]->chrg[nz]->ThresInf,&gv.bin[nd]->chrg[nz]->ThresInfVal,
			     &gv.bin[nd]->chrg[nz]->ThresSurf,&gv.bin[nd]->chrg[nz]->ThresSurfVal,
			     &gv.bin[nd]->chrg[nz]->PotSurf,INCL_TUNNEL);

		/* >>chng 01 may 09, do not use tunneling corrections for incoming electrons */
		/* >>chng 01 nov 25, add gv.bin[nd]->chrg[nz]->ThresInfInc, PvH */
		GetPotValues(nd,Zg-1,&gv.bin[nd]->chrg[nz]->ThresInfInc,&d[0],&gv.bin[nd]->chrg[nz]->ThresSurfInc,
			     &d[1],&gv.bin[nd]->chrg[nz]->PotSurfInc,NO_TUNNEL);

		ipLo = 0;
		/* >>chng 01 nov 29, nflux -> nupper so that pointer is always valid, even above nflux */
		ipHi = rfield.nupper-1;
		/* find anu[ipLo]-0.5*widflx[ipLo] <= ThresInfVal < anu[ipHi]-0.5*widflx[ipHi] */
		while( ipHi-ipLo > 1 )
		{
			long ipMd = (ipLo+ipHi)/2;
			if( gv.anumin[ipMd] > (float)gv.bin[nd]->chrg[nz]->ThresInfVal )
				ipHi = ipMd;
			else
				ipLo = ipMd;
		}
		gv.bin[nd]->chrg[nz]->ipThresInfVal = ipLo;
	}
	else
	{
		ipLo = gv.bin[nd]->chrg[nz]->ipThresInfVal;
	}

	/* remember how far the yhat, cs_pdt, fac1, and fac2 arrays were filled in */
	gv.bin[nd]->chrg[nz]->nfill = rfield.nflux;
	nfill = gv.bin[nd]->chrg[nz]->nfill;

	/* >>chng 04 feb 07, only allocate arrays from ipLo to nfill to save memory, PvH */
	FREE_SAFE( gv.bin[nd]->chrg[nz]->yhat );

	if( nfill > ipLo )
	{
		alloc_size = (size_t)((unsigned)(nfill-ipLo)*sizeof(float));

		if( ( gv.bin[nd]->chrg[nz]->yhat = (float*)MALLOC(alloc_size) ) == NULL )
			BadMalloc();

		yhat = gv.bin[nd]->chrg[nz]->yhat;
		ThresSurfVal = gv.bin[nd]->chrg[nz]->ThresSurfVal;
		DustWorkFcn = gv.bin[nd]->DustWorkFcn;
		pcase = gv.which_pe[gv.bin[nd]->matType];
		Elo = -gv.bin[nd]->chrg[nz]->PotSurf;
		ThresInfVal = gv.bin[nd]->chrg[nz]->ThresInfVal;

		/* the cost of the switch and the if statements are negligible here */
		for( i=ipLo; i < nfill; i++ )
		{
			double xv,yzero,y2,Ehi;
			xv = MAX2((anu[i] - ThresSurfVal)/DustWorkFcn,0.);

			switch( pcase )
			{
			case PE_CAR:
				/* >>refer	grain	physics	Bakes & Tielens, 1994, ApJ, 427, 822 */
				xv = xv*POW2(xv)*POW2(xv);
				yzero = xv/((1./9.e-3) + (3.7e-2/9.e-3)*xv);
				break;
			case PE_SIL:
				/* >>refer	grain	physics	Weingartner & Draine, 2000 */
				yzero = xv/(2.+10.*xv);
				break;
			default:
				fprintf( ioQQQ, " UpdatePot1 detected unknown type for PE effect: %d\n" , pcase );
				puts( "[Stop in UpdatePot1]" );
				cdEXIT(1);
			}
			if( Zg > -1 )
			{
				Ehi = anu[i] - ThresInfVal;
				y2 = ( Ehi > 0. ) ? POW2(Ehi) * (Ehi - 3.*Elo) / POW3(Ehi - Elo) : 0.;
			}
			else 
			{
				y2 = 1.;
			}
			yhat[i-ipLo] = (float)(y2*MIN2(yzero*yone[i],1.));
		}
	}

	if( ipStart == 0 )
	{
		/* >>chng 01 jan 08, ThresInf[nz] and ThresInfVal[nz] may become zero in
		 * initial stages of grain potential search, PvH */
		/* >>chng 01 oct 10, use bisection search to find ipThresInf, ipThresInfVal. On C scale now */
		ipLo = 0;
		/* >>chng 01 nov 29, nflux -> nupper so that pointer is always valid, even above nflux */
		ipHi = rfield.nupper-1;
		/* find anu[ipLo]-0.5*widflx[ipLo] <= ThresInf < anu[ipHi]-0.5*widflx[ipHi] */
		while( ipHi-ipLo > 1 )
		{
			long ipMd = (ipLo+ipHi)/2;
			if( gv.anumin[ipMd] > (float)gv.bin[nd]->chrg[nz]->ThresInf )
				ipHi = ipMd;
			else
				ipLo = ipMd;
		}
		gv.bin[nd]->chrg[nz]->ipThresInf = ipLo;
	}
	else
	{
		ipLo = gv.bin[nd]->chrg[nz]->ipThresInf;
	}

	if( Zg > -1 )
		ipLo = nfill;

	/* >>chng 04 feb 07, only allocate arrays from ipLo to nfill to save memory, PvH */
	FREE_SAFE( gv.bin[nd]->chrg[nz]->cs_pdt );

	if( nfill > ipLo )
	{
		alloc_size = (size_t)((unsigned)(nfill-ipLo)*sizeof(double));

		if( ( gv.bin[nd]->chrg[nz]->cs_pdt = (double*)MALLOC(alloc_size) ) == NULL )
			BadMalloc();

		cs_pdt = gv.bin[nd]->chrg[nz]->cs_pdt;
		c1 = -CS_PDT*(double)Zg;
		ThresInf = gv.bin[nd]->chrg[nz]->ThresInf;
		cnv_GR_pH = gv.bin[nd]->cnv_GR_pH;

		for( i=ipLo; i < nfill; i++ )
		{
			double x = (anu[i] - ThresInf)*INV_DELTA_E;
			double cs = c1*x/POW2(1.+(1./3.)*POW2(x));
			/* this cross section must be at default depletion for consistency
			 * with dstab1, hence no factor dstAbund should be applied here */
			cs_pdt[i-ipLo] = MAX2(cs,0.)*cnv_GR_pH;
		}
	}

	/* >>chng 04 feb 07, added fac1, fac2 to optimize loop for photoelectric heating, PvH */
	FREE_SAFE( gv.bin[nd]->chrg[nz]->fac1 );
	FREE_SAFE( gv.bin[nd]->chrg[nz]->fac2 );

	ipLo = gv.bin[nd]->chrg[nz]->ipThresInf;
	if( nfill > ipLo )
	{
		alloc_size = (size_t)((unsigned)(nfill-ipLo)*sizeof(double));

		if( ( gv.bin[nd]->chrg[nz]->fac1 = (double*)MALLOC(alloc_size) ) == NULL )
			BadMalloc();
		if( ( gv.bin[nd]->chrg[nz]->fac2 = (double*)MALLOC(alloc_size) ) == NULL )
			BadMalloc();

		fac1 = gv.bin[nd]->chrg[nz]->fac1;
		fac2 = gv.bin[nd]->chrg[nz]->fac2;
		PotSurf = gv.bin[nd]->chrg[nz]->PotSurf;

		for( i=ipLo; i < nfill; i++ )
		{
			double cs1,cs2,cs_tot,cool1,cool2,ehat1,ehat2;

			/* >>chng 04 jan 25, created separate routine PE_init, PvH */
			PE_init(nd,nz,i,&cs1,&cs2,&cs_tot,&cool1,&cool2,&ehat1,&ehat2);

			fac1[i-ipLo] = MAX2(cs_tot*anu[i]-cs1*cool1-cs2*cool2,0.);
			fac2[i-ipLo] = cs1*(ehat1-PotSurf) + cs2*(ehat2-PotSurf);
		}
	}

	if( ipStart == 0 )
	{
		/* >>chng 00 jul 05, determine ionization stage Z0 the ion recombines to */
		/* >>chng 04 jan 20, use ALL_STAGES here so that result remains valid throughout the model */
		/* UpdateRecomZ0(nd,nz,NONZERO_STAGES); */
		UpdateRecomZ0(nd,nz,ALL_STAGES);
	}

	/* invalidate the remaining fields */
	gv.bin[nd]->chrg[nz]->FracPop = -DBL_MAX;

	gv.bin[nd]->chrg[nz]->RSum1 = -DBL_MAX;
	gv.bin[nd]->chrg[nz]->RSum2 = -DBL_MAX;
	gv.bin[nd]->chrg[nz]->ESum1a = -DBL_MAX;
	gv.bin[nd]->chrg[nz]->ESum1b = -DBL_MAX;
	gv.bin[nd]->chrg[nz]->ESum1c = -DBL_MAX;
	gv.bin[nd]->chrg[nz]->ESum2 = -DBL_MAX;

	gv.bin[nd]->chrg[nz]->tedust = 1.f;

	gv.bin[nd]->chrg[nz]->hcon1 = -DBL_MAX;
	gv.bin[nd]->chrg[nz]->hots1 = -DBL_MAX;
	gv.bin[nd]->chrg[nz]->bolflux1 = -DBL_MAX;
	gv.bin[nd]->chrg[nz]->pe1 = -DBL_MAX;

	gv.bin[nd]->chrg[nz]->BolFlux = -DBL_MAX;
	gv.bin[nd]->chrg[nz]->GrainHeat = -DBL_MAX;
	gv.bin[nd]->chrg[nz]->GrainHeatColl = -DBL_MAX;
	gv.bin[nd]->chrg[nz]->GasHeatPhotoEl = -DBL_MAX;
	gv.bin[nd]->chrg[nz]->GasHeatTherm = -DBL_MAX;
	gv.bin[nd]->chrg[nz]->GrainCoolTherm = -DBL_MAX;
	gv.bin[nd]->chrg[nz]->ChemEnIon = -DBL_MAX;
	gv.bin[nd]->chrg[nz]->ChemEnH2 = -DBL_MAX;

	gv.bin[nd]->chrg[nz]->HeatingRate2 = -DBL_MAX;

	/* sanity check */
	ASSERT( gv.bin[nd]->chrg[nz]->ipThresInf <= gv.bin[nd]->chrg[nz]->ipThresInfVal );

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


/* this routine updates all quantities that depend on grain charge, radius and temperature
 *
 * NB NB - All global data in grain.c and grainvar.h that are charge dependent should
 *         be calculated here or in UpdatePot1 (except gv.bin[nd]->chrg[nz]->FracPop
 *         which is special).
 *
 * NB NB - the code assumes that the data calculated here may vary throughout the model,
 *         e.g. because of a dependence on grain temperature
 */
STATIC void UpdatePot2(long nd,
		       long nz)
{
	double ThermExp;

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

	/* >>chng 00 jun 19, add in loss rate due to thermionic emission of electrons, PvH */
	ThermExp = gv.bin[nd]->chrg[nz]->ThresInf*TE1RYD/gv.bin[nd]->tedust;
	/* ThermExp is guaranteed to be >= 0. */
	gv.bin[nd]->chrg[nz]->ThermRate = THERMCONST*gv.bin[nd]->ThermEff*POW2(gv.bin[nd]->tedust)*exp(-ThermExp);
#	if defined( WD_TEST2 ) || defined( IGNORE_THERMIONIC )
	gv.bin[nd]->chrg[nz]->ThermRate = 0.;
#	endif

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


/* find highest ionization stage with non-zero population */
STATIC long HighestIonStage(void)
{
	long high,
	  ion,
	  nelem;

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

	high = 0;
	for( nelem=LIMELM-1; nelem >= 0; nelem-- )
	{
		if( dense.lgElmtOn[nelem] )
		{
			for( ion=nelem+1; ion >= 0; ion-- )
			{
				if( ion == high || dense.xIonDense[nelem][ion] > 0. )
					break;
			}
			high = MAX2(high,ion);
		}
		if( nelem <= high )
			break;
	}

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


STATIC void UpdateRecomZ0(long nd,
			  long nz,
			  int lgAllIonStages)
{
	long hi_ion,
	  i,
	  ion,
	  nelem,
	  Zg;
	double d[4],
	  phi_s_up[LIMELM+1],
	  phi_s_dn[2];

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

	Zg = gv.bin[nd]->chrg[nz]->DustZ;

	hi_ion = ( lgAllIonStages ) ? LIMELM : gv.HighestIon;

	phi_s_up[0] = gv.bin[nd]->chrg[nz]->ThresSurf;
	for( i=1; i <= LIMELM; i++ )
	{
		if( i <= hi_ion )
			GetPotValues(nd,Zg+i,&d[0],&d[1],&phi_s_up[i],&d[2],&d[3],INCL_TUNNEL);
		else
			phi_s_up[i] = -DBL_MAX;
	}
	phi_s_dn[0] = gv.bin[nd]->chrg[nz]->ThresSurfInc;
	GetPotValues(nd,Zg-2,&d[0],&d[1],&phi_s_dn[1],&d[2],&d[3],NO_TUNNEL);

	/* >>chng 01 may 09, use GrainIonColl which properly tracks step-by-step charge changes */
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		if( dense.lgElmtOn[nelem] )
		{
			for( ion=0; ion <= nelem+1; ion++ )
			{
				if( lgAllIonStages || dense.xIonDense[nelem][ion] > 0. )
				{
					GrainIonColl(nd,nz,nelem,ion,phi_s_up,phi_s_dn,
						     &gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion],
						     &gv.bin[nd]->chrg[nz]->RecomEn[nelem][ion],
						     &gv.bin[nd]->chrg[nz]->ChemEn[nelem][ion]);
				}
				else
				{
					gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion] = ion;
					gv.bin[nd]->chrg[nz]->RecomEn[nelem][ion] = 0.f;
					gv.bin[nd]->chrg[nz]->ChemEn[nelem][ion] = 0.f;
				}
			}
		}
	}

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

STATIC void GetPotValues(long nd,
			 long Zg,
			 /*@out@*/ double *ThresInf,
			 /*@out@*/ double *ThresInfVal,
			 /*@out@*/ double *ThresSurf,
			 /*@out@*/ double *ThresSurfVal,
			 /*@out@*/ double *PotSurf,
			 int lgUseTunnelCorr)
{
	double dstpot,
	  dZg = (double)Zg,
	  IP_v;

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

	/* >>chng 01 may 07, this routine now completely supports the hybrid grain charge model,
	 * the change for this routine is that now it is only fed integral charge states; calculation
	 * of IP has also been changed in accordance with Weingartner & Draine, 2001, PvH */

	/* this is average grain potential in Rydberg */
	dstpot = CHRG2POT(dZg);

	/* >>chng 01 mar 20, include O(a^-2) correction terms in ionization potential */
	/* these are correction terms for the ionization potential that are
	 * important for small grains. See Weingartner & Draine, 2000, Eq. 2 */
	IP_v = gv.bin[nd]->DustWorkFcn + dstpot - 0.5*ONE_ELEC + (dZg+2.)*AC0/gv.bin[nd]->AvRadius*ONE_ELEC;

	/* >>chng 01 mar 01, use new expresssion for ThresInfVal, ThresSurfVal following the discussion
	 * with Joe Weingartner. Also include the Schottky effect (see 
	 * >>refer	grain	physics	Spitzer, 1948, ApJ, 107, 6,
	 * >>refer	grain	physics	Draine & Sutin, 1987, ApJ, 320, 803), PvH */
	if( Zg <= -1 )
	{
		pot_type pcase = gv.which_pot[gv.bin[nd]->matType];
		double Emin,IP;

		IP = gv.bin[nd]->DustWorkFcn - gv.bin[nd]->BandGap + dstpot - 0.5*ONE_ELEC;
		switch( pcase )
		{
		case POT_CAR:
			IP -= AC1G/(gv.bin[nd]->AvRadius+AC2G)*ONE_ELEC;
			break;
		case POT_SIL:
			/* do nothing */
			break;
		default:
			fprintf( ioQQQ, " GetPotValues detected unknown type for ionization pot: %d\n" , pcase );
			puts( "[Stop in GetPotValues]" );
			cdEXIT(1);
		}

		/* prevent valence electron from becoming less bound than attached electron; this
		 * can happen for very negative, large graphitic grains and is not physical, PvH */
		IP_v = MAX2(IP,IP_v);

		if( Zg < -1 )
		{
			/* >>chng 01 apr 20, use improved expression for tunneling effect, PvH */
			double help = fabs(dZg+1);
			/* this is the barrier height solely due to the Schottky effect */
			Emin = -ThetaNu(help)*ONE_ELEC;
			if( lgUseTunnelCorr )
			{
				/* this is the barrier height corrected for tunneling effects */
				Emin *= 1. - 2.124e-4/(pow(gv.bin[nd]->AvRadius,0.45f)*pow(help,0.26));
			}
		}
		else
		{
			Emin = 0.;
		}

		*ThresInf = IP - Emin;
		*ThresInfVal = IP_v - Emin;
		*ThresSurf = *ThresInf;
		*ThresSurfVal = *ThresInfVal;
		*PotSurf = Emin;
	}
	else
	{
		*ThresInf = IP_v;
		*ThresInfVal = IP_v;
		*ThresSurf = *ThresInf - dstpot;
		*ThresSurfVal = *ThresInfVal - dstpot;
		*PotSurf = dstpot;
	}

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


/* given grain nd in charge state nz, and incoming ion (nelem,ion),
 * detemine outgoing ion (nelem,Z0) and chemical energy ChEn released
 * ChemEn is net contribution of ion recombination to grain heating */
STATIC void GrainIonColl(long int nd,
			 long int nz,
			 long int nelem,
			 long int ion,
			 const double phi_s_up[], /* phi_s_up[LIMELM+1] */
			 const double phi_s_dn[], /* phi_s_dn[2] */
			 /*@out@*/long *Z0,
			 /*@out@*/float *ChEn,
			 /*@out@*/float *ChemEn)
{
	long Zg;
	double d[4];
	double phi_s;
	
#	ifdef DEBUG_FUN
	fputs( "<+>GrainIonColl()\n", debug_fp );
#	endif

	long save = ion;
	if( ion > 0 && rfield.anu[Heavy.ipHeavy[nelem][ion-1]-1] > (float)phi_s_up[0] )
	{
		/* ion will get electron(s) */
		*ChEn = 0.f;
		*ChemEn = 0.f;
		Zg = gv.bin[nd]->chrg[nz]->DustZ;
		phi_s = phi_s_up[0];
		do 
		{
			*ChEn += rfield.anu[Heavy.ipHeavy[nelem][ion-1]-1] - (float)phi_s;
			*ChemEn += rfield.anu[Heavy.ipHeavy[nelem][ion-1]-1];
			/* this is a correction for the imperfections in the n-charge state model:
			 * since the transfer gets modeled as n single-electron transfers, instead of one
			 * n-electron transfer, a correction for the difference in binding energy is needed */
			*ChemEn -= (float)(phi_s - phi_s_up[0]);
			--ion;
			++Zg;
			phi_s = phi_s_up[save-ion];
		} while( ion > 0 && rfield.anu[Heavy.ipHeavy[nelem][ion-1]-1] > (float)phi_s );

		*Z0 = ion;
	}
	else if( ion <= nelem && gv.bin[nd]->chrg[nz]->DustZ > gv.bin[nd]->LowestZg &&
		 rfield.anu[Heavy.ipHeavy[nelem][ion]-1] < (float)phi_s_dn[0] )
	{
		/* grain will get electron(s) */
		*ChEn = 0.f;
		*ChemEn = 0.f;
		Zg = gv.bin[nd]->chrg[nz]->DustZ;
		phi_s = phi_s_dn[0];
		do 
		{
			*ChEn += (float)phi_s - rfield.anu[Heavy.ipHeavy[nelem][ion]-1];
			*ChemEn -= rfield.anu[Heavy.ipHeavy[nelem][ion]-1];
			/* this is a correction for the imperfections in the n-charge state model:
			 * since the transfer gets modeled as n single-electron transfers, instead of one
			 * n-electron transfer, a correction for the difference in binding energy is needed */
			*ChemEn += (float)(phi_s - phi_s_dn[0]);
			++ion;
			--Zg;

			if( ion-save < 2 )
				phi_s = phi_s_dn[ion-save];
			else
				GetPotValues(nd,Zg-1,&d[0],&d[1],&phi_s,&d[2],&d[3],NO_TUNNEL);

		} while( ion <= nelem && Zg > gv.bin[nd]->LowestZg &&
			 rfield.anu[Heavy.ipHeavy[nelem][ion]-1] < (float)phi_s );
		*Z0 = ion;
	}
	else
	{
		/* nothing happens */
		*ChEn = 0.f;
		*ChemEn = 0.f;
		*Z0 = ion;
	}
/*  	printf(" GrainIonColl: nelem %ld ion %ld -> %ld, ChEn %.6f\n",nelem,save,*Z0,*ChEn); */

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


/* initialize grain-ion charge transfer rates on grain species nd */
STATIC void GrainChrgTransferRates(long nd)
{
#	ifdef DEBUG_FUN
	fputs( "<+>GrainChrgTransferRates()\n", debug_fp );
#	endif

#	ifndef IGNORE_GRAIN_ION_COLLISIONS
	long nz;
	double fac0 = STICK_ION*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;

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

		if( fac1 == 0. )
			continue;

		for( ion=0; ion <= LIMELM; ion++ )
		{
			long nelem;
			double eta, fac2, xi;

			/* >>chng 00 jul 19, replace classical results with results including image potential
			 * to correct for polarization of the grain as charged particle approaches. */
			GrainScreen(ion,nd,nz,&eta,&xi);

			fac2 = eta*fac1;

			if( fac2 == 0. )
				continue;

			for( nelem=MAX2(0,ion-1); nelem < LIMELM; nelem++ )
			{
				if( dense.lgElmtOn[nelem] && ion != gptr->RecomZ0[nelem][ion] )
				{
					gv.GrainChTrRate[nelem][ion][gptr->RecomZ0[nelem][ion]] +=
						/* ionbal.lgGrainIonRecom is 1 by default, set to 0 with
						 * no grain neutralization command */
						(float)(fac2*DoppVel.AveVel[nelem]) * ionbal.lgGrainIonRecom;
				}
			}
		}
	}
#	endif

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


/* this routine updates all grain quantities that depend on radius,
 * except gv.dstab and gv.dstsc which are updated in GrainUpdateRadius2() */
STATIC void GrainUpdateRadius1(void)
{
	long nelem,
	  nd;

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

	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		gv.elmSumAbund[nelem] = 0.f;
	}

	/* grain abundance may be a function of depth */
	for( nd=0; nd < gv.nBin; nd++ )
	{
		/* refresh grain abundance if variable grain abundances enabled */
		if( gv.bin[nd]->dstAbund <= 0.f || gv.bin[nd]->lgDustVary )
		{
			gv.bin[nd]->dstAbund = (float)(gv.bin[nd]->dstfactor*GrnVryDpth(nd));
			ASSERT( gv.bin[nd]->dstAbund > 0.f );
		}

		/* grain unit conversion, <unit>/H (default depl) -> <unit>/cm^3 (actual depl) */
		gv.bin[nd]->cnv_H_pCM3 = dense.gas_phase[ipHYDROGEN]*gv.bin[nd]->dstAbund;
		gv.bin[nd]->cnv_CM3_pH = 1./gv.bin[nd]->cnv_H_pCM3;
		/* grain unit conversion, <unit>/cm^3 (actual depl) -> <unit>/grain */
		gv.bin[nd]->cnv_CM3_pGR = gv.bin[nd]->cnv_H_pGR/gv.bin[nd]->cnv_H_pCM3;
		gv.bin[nd]->cnv_GR_pCM3 = 1./gv.bin[nd]->cnv_CM3_pGR;

		/* >>chng 01 dec 05, calculate the number density of each element locked in grains,
		 * summed over all grain bins. this number uses the actual depletion of the grains
		 * and is already multiplied with hden, units cm^-3. */
		for( nelem=0; nelem < LIMELM; nelem++ )
		{
			gv.elmSumAbund[nelem] += gv.bin[nd]->elmAbund[nelem]*(float)gv.bin[nd]->cnv_H_pCM3;
		}
	}

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


/* this routine adds all the grain opacities in gv.dstab and gv.dstsc, this could not be
 * done in GrainUpdateRadius1 since charge and FracPop must be converged first */
STATIC void GrainUpdateRadius2(int lgAnyNegCharge)
{
	int lgChangeCS_PDT;
	long i,
	  nd,
	  nz;

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

	/* if either previous or this iteration has negative charges, cs_pdt is likely to change */
	lgChangeCS_PDT = gv.lgAnyNegCharge || lgAnyNegCharge;

	if( rfield.nflux > gv.nfill || ( gv.lgAnyDustVary && nzone != gv.nzone ) || lgChangeCS_PDT )
	{
		/* remember how far the dtsab and dstsc arrays were filled in */
		gv.nfill = rfield.nflux;
		gv.lgAnyNegCharge = lgAnyNegCharge;

		for( i=0; i < rfield.nupper; i++ )
		{
			gv.dstab[i] = 0.;
			gv.dstsc[i] = 0.;
		}

		/* >>chng 01 mar 26, from nupper to nflux */
		for( i=0; i < rfield.nflux; i++ )
		{
			for( nd=0; nd < gv.nBin; nd++ )
			{
				/* these are total absorption and scattering cross sections,
				 * the latter should contain the asymmetry factor (1-g) */
				/* this is effective area per proton, scaled by depletion
				 * dareff(nd) = darea(nd) * dstAbund(nd) */
				/* grain abundance may be a function of depth */
				/* >>chng 02 dec 30, separated scattering cross section and asymmetry factor, PvH */
				gv.dstab[i] += gv.bin[nd]->dstab1[i]*gv.bin[nd]->dstAbund;
				gv.dstsc[i] += gv.bin[nd]->pure_sc1[i]*gv.bin[nd]->asym[i]*gv.bin[nd]->dstAbund;

				/* >>chng 01 mar 24, added in cs for photodetachment from negative grains, PvH */
				/* >>chng 01 may 07, use two charge state approximation */
				/* >>chng 01 may 30, replace upper limit of loop gv.bin[nd]->nChrg -> nz_max */
				/* >>chng 04 jan 26, change back to gv.bin[nd]->nChrg, use "if" for clarity, PvH */
				for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
				{
					ChargeBin *gptr = gv.bin[nd]->chrg[nz];
					long ipLo = gptr->ipThresInf;

					/* >>chng 01 may 30, cs_pdt is only non-zero for negative charge states */
					if( gptr->DustZ <= -1 && i >= ipLo )
						gv.dstab[i] +=
							gptr->FracPop*gptr->cs_pdt[i-ipLo]*gv.bin[nd]->dstAbund;
					}
			}
			/* this must be positive, zero in case of uncontrolled underflow */
			ASSERT( gv.dstab[i] > 0. && gv.dstsc[i] > 0. );
		}
	}

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


/* GrainTemperature computes grains temperature, and gas cooling */
STATIC void GrainTemperature(long int nd,
			     /*@out@*/ float *dccool,
			     /*@out@*/ double *hcon,
			     /*@out@*/ double *hots,
			     /*@out@*/ double *hla)
{
	long int i,
	  ipLya,
	  nz;
	double csrheat[NCHU],
	  csthcool[NCHU],
	  EhatThermionic,
	  norm,
	  rate,
	  x,
	  y;
	float dcheat,
	  cscheat[NCHU];

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

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

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

	/* >>chng 01 may 07, this routine now completely supports the hybrid grain
	 * charge model, and the average charge state is not used anywhere anymore, PvH */

	/* direct heating by incident continuum (all energies) */
	*hcon = 0.;
	/* heating by diffuse ots fields */
	*hots = 0.;
	/* heating by Ly alpha alone, for output only, is already included in hots */
	*hla = 0.;

	ipLya = EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont - 1;

	/* integrate over ionizing continuum; energy goes to dust and gas
	 * GasHeatPhotoEl is what heats the gas */
	gv.bin[nd]->GasHeatPhotoEl = 0.;

	gv.bin[nd]->GrainCoolTherm = 0.;
	gv.bin[nd]->thermionic = 0.;

	dcheat = 0.f;
	*dccool = 0.f;

	gv.bin[nd]->BolFlux = 0.;

	/* >>chng 04 jan 25, moved initialization of phiTilde to qheat_init(), PvH */

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

		double hcon1 = 0.;
		double hots1 = 0.;
		double hla1 = 0.;
		double bolflux1 = 0.;
		double pe1 = 0.;

		/* >>chng 04 may 31, introduced lgReEvaluate2 to save time when iterating Tdust, PvH */
		int lgReEvaluate1 = gptr->hcon1 < 0.;
		int lgReEvaluate2 = gptr->hots1 < 0.;
		int lgReEvaluate = lgReEvaluate1 || lgReEvaluate2;

		/* integrate over incident continuum for non-ionizing energies */
		if( lgReEvaluate )
		{
			for( i=0; i < MIN2(gptr->ipThresInf,rfield.nflux); i++ )
			{
				double fac = gv.bin[nd]->dstab1[i]*rfield.anu[i];

				if( lgReEvaluate1 )
					hcon1 += rfield.flux[i]*fac;

				if( lgReEvaluate2 )
				{
					hots1 += rfield.SummedDif[i]*fac;
#					ifndef NDEBUG
					bolflux1 += rfield.SummedCon[i]*fac;
#					endif
				}
			}
		}

		/* >>chng 01 mar 02, use new expresssions for grain cooling and absorption
		 * cross sections following the discussion with Joe Weingartner, PvH */
		/* >>chng 04 feb 07, use fac1, fac2 to optimize this loop, PvH */
		if( lgReEvaluate )
		{
			for( i=gptr->ipThresInf; i < rfield.nflux; i++ )
			{
				/* this is heating by incident radiation field */
				if( lgReEvaluate1 )
					hcon1 += rfield.flux[i]*gptr->fac1[i-ipLo];

				if( lgReEvaluate2 )
				{
					/* this is heating by all diffuse fields:
					 * SummedDif has all continua and lines */
					hots1 += rfield.SummedDif[i]*gptr->fac1[i-ipLo];
					/*  GasHeatPhotoEl is rate grain photoionization heats the gas */
#ifdef WD_TEST2
					pe1 += rfield.flux[i]*gptr->fac2[i-ipLo];
#else
					pe1 += rfield.SummedCon[i]*gptr->fac2[i-ipLo];
#endif
#					ifndef NDEBUG
					bolflux1 += rfield.SummedCon[i]*gv.bin[nd]->dstab1[i]*rfield.anu[i];
					if( gptr->DustZ <= -1 )
						bolflux1 += rfield.SummedCon[i]*gptr->cs_pdt[i-ipLo]*rfield.anu[i];
#					endif
				}
			}
		}

		if( lgReEvaluate1 )
		{
			/* >>chng 04 feb 07, remember hcon1 for possible later use, PvH */
			gptr->hcon1 = hcon1;
		}
		else
		{
			hcon1 = gptr->hcon1;
		}

		if( lgReEvaluate2 )
		{
			gptr->hots1 = hots1;
			gptr->bolflux1 = bolflux1;
			gptr->pe1 = pe1;
		}
		else
		{
			hots1 = gptr->hots1;
			bolflux1 = gptr->bolflux1;
			pe1 = gptr->pe1;
		}

		/*  heating by Ly A on dust in this zone,
		 *  only used for printout; Ly-a is already in OTS fields */
		/* >>chng 00 apr 18, moved calculation of hla, by PvH */
		/* >>chng 04 feb 01, moved calculation of hla1 outside loop for optimization, PvH */
		if( ipLya < MIN2(gptr->ipThresInf,rfield.nflux) )
		{
			hla1 = rfield.otslin[ipLya]*gv.bin[nd]->dstab1[ipLya]*0.75;
		}
		else if( ipLya < rfield.nflux )
		{
			/* >>chng 00 apr 18, include photo-electric effect, by PvH */
			hla1 = rfield.otslin[ipLya]*gptr->fac1[ipLya-ipLo];
		}
		else
		{
			hla1 = 0.;
		}

		ASSERT( hcon1 >= 0. && hots1 >= 0. && hla1 >= 0. && bolflux1 >= 0. && pe1 >= 0. );

		*hcon += gptr->FracPop*hcon1;
		*hots += gptr->FracPop*hots1;
		*hla += gptr->FracPop*hla1;
		gv.bin[nd]->BolFlux += gptr->FracPop*bolflux1;
		if( gv.lgDHetOn )
			gv.bin[nd]->GasHeatPhotoEl += gptr->FracPop*pe1;

		/* >>chng 04 feb 08, remember radiative heating for each charge state, PvH */
		csrheat[nz] = (hcon1 + hots1)*EN1RYD*gv.bin[nd]->cnv_H_pCM3;

#		ifndef NDEBUG
		if( trace.lgTrace && trace.lgDustBug )
		{
			fprintf( ioQQQ, "    Zg %ld bolflux: %.4e\n", gptr->DustZ,
			  gptr->FracPop*bolflux1*EN1RYD*gv.bin[nd]->cnv_H_pCM3 );
		}
#		endif

		/* add in thermionic emissions (thermal evaporation of electrons), it gives a cooling
		 * term for the grain. thermionic emissions will not be treated separately in quantum
		 * heating since they are only important when grains are heated to near-sublimation 
		 * temperatures; under those conditions quantum heating effects will never be important.
		 * in order to maintain energy balance they will be added to the ion contribution though */
		/* ThermRate is normalized per cm^2 of grain surface area, scales with total grain area */
		rate = gptr->FracPop*gptr->ThermRate*gv.bin[nd]->IntArea*gv.bin[nd]->cnv_H_pCM3;
		/* >>chng 01 mar 02, PotSurf[nz] term was incorrectly taken into account, PvH */
		EhatThermionic = 2.*BOLTZMANN*gv.bin[nd]->tedust + MAX2(gptr->PotSurf*EN1RYD,0.);
		gv.bin[nd]->GrainCoolTherm += rate * (EhatThermionic + gptr->ThresSurf*EN1RYD);
		gv.bin[nd]->thermionic += rate * (EhatThermionic - gptr->PotSurf*EN1RYD);

		/* >>chng 04 feb 08, remember thermionic cooling for each charge state, PvH */
		csthcool[nz] = gptr->ThermRate*gv.bin[nd]->IntArea*gv.bin[nd]->cnv_H_pCM3*
			(EhatThermionic + gptr->ThresSurf*EN1RYD);
	}

	/* norm is used to convert all heating rates to erg/cm^3/s */
	norm = EN1RYD*gv.bin[nd]->cnv_H_pCM3;

	/* hcon is radiative heating by incident radiation field */
	*hcon *= norm;

	/* hots is total heating of the grain by diffuse fields */
	*hots *= norm;

	/* heating by Ly alpha alone, for output only, is already included in hots */
	*hla *= norm;

	gv.bin[nd]->BolFlux *= norm;

	/* heating by thermal collisions with gas does work
	 * DCHEAT is grain collisional heating by gas
	 * DCCOOL is gas cooling due to collisions with grains
	 * they are different since grain surface recombinations
	 * heat the grains, but do not cool the gas ! */
	/* >>chng 03 nov 06, moved call after renorm of BolFlux, so that GrainCollHeating can look at it, PvH */
	GrainCollHeating(nd,&dcheat,dccool,cscheat);

	/* GasHeatPhotoEl is what heats the gas */
	gv.bin[nd]->GasHeatPhotoEl *= norm;

	if( gv.lgLeidenBakesPAH_heat )
	{
		/* this is a dirty hack to get BT94 PE heating rate
		 * for PAH's included, for Lorentz Center 2004 PDR meeting, PvH */
		/*>>>refer	PAH	heating	Bakes, E.L.O., & Tielens, A.G.G.M. 1994, ApJ, 427, 822 */
		gv.bin[nd]->GasHeatPhotoEl += 1.e-24*hmi.UV_Cont_rel2_Habing_TH85_depth*
			dense.gas_phase[ipHYDROGEN]*(4.87e-2/(1.0+4e-3*pow(hmi.UV_Cont_rel2_Habing_TH85_depth*
			sqrt(phycon.te)/dense.eden,0.73)) + 3.65e-2*pow(phycon.te/1.e4,0.7)/
			(1.+2.e-4*(hmi.UV_Cont_rel2_Habing_TH85_depth*sqrt(phycon.te)/dense.eden)))/gv.nBin;
	}

	/* >>chng 01 nov 29, removed next statement, PvH */
	/*  dust often hotter than gas during initial TE search */
	/* if( nzone <= 2 ) */
	/* 	dcheat = MAX2(0.f,dcheat); */

	/*  find power absorbed by dust and resulting temperature
	 *
	 * hcon is heating from incident continuum (all energies)
	 * hots is heating from ots continua and lines
	 * dcheat is net grain collisional and chemical heating by
	 *    particle collisions and recombinations
	 * GrainCoolTherm is grain cooling by thermionic emissions
	 *
	 * GrainHeat is net heating of this grain type,
	 *    to be balanced by radiative cooling */
	gv.bin[nd]->GrainHeat = *hcon + *hots + dcheat - gv.bin[nd]->GrainCoolTherm;

	/* remember collisional heating for this grain species */
	gv.bin[nd]->GrainHeatColl = dcheat;

	/* >>chng 04 may 31, replace ASSERT of GrainHeat > 0. with if-statement and let
	 * GrainChargeTemp sort out the consquences of GrainHeat becoming negative, PvH */
	/* in case where the thermionic rates become very large,
	 * or collisional cooling dominates, this may become negative */
	if( gv.bin[nd]->GrainHeat > 0. )
	{
		int lgOutOfBounds;
		/*  now find temperature, GrainHeat is sum of total heating of grain
		 *  >>chng 97 jul 17, divide by abundance here */
		x = log(MAX2(DBL_MIN,gv.bin[nd]->GrainHeat*gv.bin[nd]->cnv_CM3_pH));
		/* >>chng 96 apr 27, as per Peter van Hoof comment */
		splint_safe(gv.bin[nd]->dstems,gv.dsttmp,gv.bin[nd]->dstslp,NDEMS,x,&y,&lgOutOfBounds);
		gv.bin[nd]->tedust = (float)exp(y);
	}
	else
	{
		gv.bin[nd]->GrainHeat = -1.;
		gv.bin[nd]->tedust = -1.;
	}

	if( thermal.ConstGrainTemp > 0. )
	{
		int lgOutOfBounds;
		/* use temperature set with constant grain temperature command */
		gv.bin[nd]->tedust = thermal.ConstGrainTemp;
		/* >>chng 04 jun 01, make sure GrainHeat is consistent with value of tedust, PvH */
		x = log(gv.bin[nd]->tedust);
		splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,x,&y,&lgOutOfBounds);
		gv.bin[nd]->GrainHeat = exp(y)*gv.bin[nd]->cnv_H_pCM3;
	}

	/*  save for later possible printout */
	gv.bin[nd]->TeGrainMax = (float)MAX2(gv.bin[nd]->TeGrainMax,gv.bin[nd]->tedust);

#if 0
	/* >>chng 04 feb 08, calculate grain temperature for each charge state, PvH */
	for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
	{
		int lgOutOfBounds;

		gv.bin[nd]->chrg[nz]->GrainHeat = csrheat[nz] + cscheat[nz] - csthcool[nz];

		ASSERT( gv.bin[nd]->chrg[nz]->GrainHeat > 0. );

		/*  now find temperature, GrainHeat is sum of total heating of this charge state */
		x = log(MAX2(DBL_MIN,gv.bin[nd]->chrg[nz]->GrainHeat*gv.bin[nd]->cnv_CM3_pH));
		splint_safe(gv.bin[nd]->dstems,gv.dsttmp,gv.bin[nd]->dstslp,NDEMS,x,&y,&lgOutOfBounds);
		gv.bin[nd]->chrg[nz]->tedust = (float)exp(y);
	}

#	ifndef NDEBUG
	{
		double z = 0.;
		for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
		{
			z += gv.bin[nd]->chrg[nz]->FracPop*gv.bin[nd]->chrg[nz]->GrainHeat;
		}
		ASSERT( fabs(gv.bin[nd]->GrainHeat/z - 1.) <= 10.*sqrt(gv.bin[nd]->nChrg)*FLT_EPSILON );
	}
#	endif
#endif

	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "  >GrainTemperature finds %s Tdst %.5e hcon %.4e ",
			 gv.bin[nd]->chDstLab, gv.bin[nd]->tedust, *hcon);
		fprintf( ioQQQ, "hots %.4e dcheat %.4e GrainCoolTherm %.4e\n", 
			 *hots, dcheat, gv.bin[nd]->GrainCoolTherm );
#if 0
		fprintf( ioQQQ, "     Charge State Temperature" );
		for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
		{
			fprintf( ioQQQ, " Zg %ld Tgr %.5e",
				 gv.bin[nd]->chrg[nz]->DustZ, gv.bin[nd]->chrg[nz]->tedust );
		}
		fprintf( ioQQQ, "\n" );
#endif
	}

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


/* helper routine for initializing quantities related to the photo-electric effect */
void PE_init(long nd,
	     long nz,
	     long i,
	     /*@out@*/ double *cs1,
	     /*@out@*/ double *cs2,
	     /*@out@*/ double *cs_tot,
	     /*@out@*/ double *cool1,
	     /*@out@*/ double *cool2,
	     /*@out@*/ double *ehat1,
	     /*@out@*/ double *ehat2)
{
	ChargeBin *gptr = gv.bin[nd]->chrg[nz];
	long ipLo1 = gptr->ipThresInfVal;
	long ipLo2 = gptr->ipThresInf;

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

	/* sanity checks */
	ASSERT( nd >= 0 && nd < gv.nBin );
	ASSERT( nz >= 0 && nz < gv.bin[nd]->nChrg );
	ASSERT( i >= 0 && i < rfield.nflux );

	/* effective cross sections for photo-ejection, 1 = valence band, 2 = conduction band */
	*cs1 = ( i >= ipLo1 ) ? gv.bin[nd]->dstab1[i]*gptr->yhat[i-ipLo1] : 0.;
	*cs2 = ( gptr->DustZ <= -1 && i >= ipLo2 ) ? gptr->cs_pdt[i-ipLo2] : 0.;
	*cs_tot = gv.bin[nd]->dstab1[i] + *cs2;

	/* >>chng 00 jul 17, use description of Weingartner & Draine, 2000 */
	/* ehat1a,1b,2 is the average energy of the esacping electron at the grain surface */
	if( gptr->DustZ > -1 ) 
	{
		/* contribution from valence band */
		double Elo = -gptr->PotSurf;
		double Ehi = rfield.anu[i] - gptr->ThresInfVal;
		*ehat1 = ( Ehi > 0. ) ?
			0.5*Ehi*(Ehi-2.*Elo)/(Ehi-3.*Elo) + gptr->PotSurf : 0.;
		*ehat2 = 0.;
	}
	else 
	{
		/* contribution from valence band */
		*ehat1 = 0.5*MAX2(rfield.anu[i] - gptr->ThresSurfVal,0.);
		/* contribution from conduction band */
		*ehat2 = MAX2(rfield.anu[i] - gptr->ThresInf,0.);
	}

	/* cool1,2 is the amount by which photo-ejection cools the grain */
	*cool1 = gptr->ThresSurfVal + *ehat1;
	/* this term in cool1 describes the heating by an electron that de-excites
	 * from the conduction band into the valence band. it is assumed that all
	 * energy is absorbed by the grain, i.e. no photons are emitted */
	/* >>chng 01 nov 27, changed de-excitation energy to conserve energy, PvH */
	if( gptr->DustZ <= -1 )
		*cool1 += gptr->ThresSurf - gptr->ThresSurfVal;

	*cool2 = gptr->ThresSurf + *ehat2;

	/* sanity checks */
	ASSERT( *ehat1 > 0. || ( *ehat1 == 0. && *cs1 == 0. ) );
	ASSERT( *ehat2 > 0. || ( *ehat2 == 0. && *cs2 == 0. ) );
	ASSERT( *cool1 >= 0. && *cool2 >= 0. );

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


#if 1
/* GrainCollHeating compute grains collisional heating cooling */
STATIC void GrainCollHeating(long int nd,
			     /*@out@*/ float *dcheat,
			     /*@out@*/ float *dccool,
			     /*@out@*/ float cscheat[]) /* cscheat[NCHU] */
{
	long int ion,
	  nelem,
	  nz;
	H2_type ipH2;
	double Accommodation,
	  CollisionRateElectr,      /* rate electrons strike grains */
	  CollisionRateMol,         /* rate molecules strike grains */
	  CollisionRateIon,         /* rate ions strike grains */
	  CoolTot,
	  CoolBounce,
	  CoolEmitted,
	  CoolElectrons,
	  CoolMolecules,
	  CoolPotential,
	  CoolPotentialGas,
	  eta,
	  HeatTot,
	  HeatBounce,
	  HeatCollisions,
	  HeatElectrons,
	  HeatIons,
	  HeatMolecules,
	  HeatRecombination, /* sum of abundances of ions times velocity times ionization potential times eta */
	  HeatChem,
	  HeatCor,
	  Stick,
	  ve,
	  WeightMol,
	  xi;

	/* energy deposited into grain by formation of a single H2 molecule, in eV,
	 * >>refer  Takahashi J., Uehara H., 2001, ApJ, 561, 843 */
	const double H2_FORMATION_GRAIN_HEATING[H2_TOP] = { 0.20, 0.4, 1.72 };

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


	/* >>chng 01 may 07, this routine now completely supports the hybrid grain
	 * charge model, and the average charge state is not used anywhere anymore, PvH */

	/* this subroutine evaluates the gas heating-cooling rate
	 * (erg cm^-3 s^-1) due to grain gas collisions.
	 * the net effect can be positive or negative,
	 * depending on whether the grains or gas are hotter
	 * the physics is described in 
	 * >>refer	grain	physics	Baldwin, Ferland, Martin et al., 1991, ApJ 374, 580 */

	HeatTot = 0.;
	CoolTot = 0.;

	HeatIons = 0.;

	gv.bin[nd]->ChemEn = 0.;

	/* loop over the charge states */
	for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
	{
		ChargeBin *gptr = gv.bin[nd]->chrg[nz];

		/* HEAT1 will be rate collisions heat the grain
		 * COOL1 will be rate collisions cool the gas kinetics */
		double Heat1 = 0.;
		double Cool1 = 0.;
		double ChemEn1 = 0.;

		/* ============================================================================= */
		/* heating/cooling due to neutrals and positive ions */

		/* loop over all stages of ionization */
		for( ion=0; ion <= LIMELM; ion++ )
		{
			/* this is heating of grains due to recombination energy of species,
			 * and assumes that every ion is fully neutralized upon striking the grain surface.
			 * all radiation produced in the recombination process is absorbed within the grain
			 *
			 * ion=0 are neutrals, ion=1 are single ions, etc
			 * each population is weighted by the AVERAGE velocity
			 * */
			CollisionRateIon = 0.;
			CoolPotential = 0.;
			CoolPotentialGas = 0.;
			HeatRecombination = 0.;
			HeatChem = 0.;

			/* >>chng 00 jul 19, replace classical results with results including image potential
			 * to correct for polarization of the grain as charged particle approaches. */
			GrainScreen(ion,nd,nz,&eta,&xi);

			for( nelem=MAX2(0,ion-1); nelem < LIMELM; nelem++ )
			{
				if( dense.lgElmtOn[nelem] && dense.xIonDense[nelem][ion] > 0. )
				{
					double CollisionRateOne;

					/* >>chng 00 apr 05, use correct accomodation coefficient, by PvH
					 * the coefficient is defined at the end of appendix A.10 of BFM
					 * assume ion sticking prob is unity */
#if defined( IGNORE_GRAIN_ION_COLLISIONS )
					Stick = 0.;
#elif defined( WD_TEST2 )
					Stick = ( ion == gptr->RecomZ0[nelem][ion] ) ?
						0. : STICK_ION;
#else
					Stick = ( ion == gptr->RecomZ0[nelem][ion] ) ?
						gv.bin[nd]->AccomCoef[nelem] : STICK_ION;
#endif
					/* this is rate with which charged ion strikes grain */
					/* >>chng 00 may 02, this had left 2./SQRTPI off */
					/* >>chng 00 may 05, use average speed instead of 2./SQRTPI*Doppler, PvH */
					CollisionRateOne = Stick*dense.xIonDense[nelem][ion]*DoppVel.AveVel[nelem];
					CollisionRateIon += CollisionRateOne;
					/* >>chng 01 nov 26, use PotSurfInc when appropriate:
					 * the values for the surface potential used here make it
					 * consistent with the rest of the code and preserve energy.
					 * NOTE: For incoming particles one should use PotSurfInc with
					 * Schottky effect for positive ion, for outgoing particles
					 * one should use PotSurf for Zg+ion-Z_0-1 (-1 because PotSurf
					 * assumes electron going out), these corrections are small
					 * and will be neglected for now, PvH */
					if( ion >= gptr->RecomZ0[nelem][ion] )
					{
						CoolPotential += CollisionRateOne * (double)ion *
							gptr->PotSurf;
						CoolPotentialGas += CollisionRateOne *
							(double)gptr->RecomZ0[nelem][ion] *
							gptr->PotSurf;
					}
					else
					{
						CoolPotential += CollisionRateOne * (double)ion *
							gptr->PotSurfInc;
						CoolPotentialGas += CollisionRateOne *
							(double)gptr->RecomZ0[nelem][ion] *
							gptr->PotSurfInc;
					}
					/* this is sum of all energy liberated as ion recombines to Z0 in grain */
					/* >>chng 00 jul 05, subtract energy needed to get 
					 * electron out of grain potential well, PvH */
					/* >>chng 01 may 09, chemical energy now calculated in GrainIonColl, PvH */
					HeatRecombination += CollisionRateOne *
						gptr->RecomEn[nelem][ion];
					HeatChem += CollisionRateOne * gptr->ChemEn[nelem][ion];
				}
			}

			/* >>chng 00 may 01, boltzmann factor had multiplied all of factor instead
			 * of only first and last term.  pvh */

			/* equation 29 from Balwin et al 91 */
			/* this is direct collision rate, 2kT * xi, first term in eq 29 */
			HeatCollisions = CollisionRateIon * 2.*BOLTZMANN*phycon.te*xi;
			/* this is change in energy due to charge acceleration within grain's potential 
			 * this is exactly balanced by deceleration of incoming electrons and accelaration
			 * of outgoing photo-electrons and thermionic emissions; all these terms should
			 * add up to zero (total charge of grain should remain constant) */
			CoolPotential *= eta*EN1RYD;
			CoolPotentialGas *= eta*EN1RYD;
			/* this is recombination energy released within grain */
			HeatRecombination *= eta*EN1RYD;
			HeatChem *= eta*EN1RYD;
			/* energy carried away by neutrals after recombination, so a cooling term */
			CoolEmitted = CollisionRateIon * 2.*BOLTZMANN*gv.bin[nd]->tedust*eta;

			/* total GraC 0 in the emission line output */
			Heat1 += HeatCollisions - CoolPotential + HeatRecombination - CoolEmitted;

			/* rate kinetic energy lost from gas - gas cooling - eq 32 in BFM */
			/* this GrGC 0 in the main output */
			/* >>chng 00 may 05, reversed sign of gas cooling contribution */
			Cool1 += HeatCollisions - CoolEmitted - CoolPotentialGas;

			ChemEn1 += HeatChem;
		}

		/* remember grain heating by ion collisions for quantum heating treatment */
		HeatIons += gptr->FracPop*Heat1;

		if( trace.lgTrace && trace.lgDustBug )
		{
			fprintf( ioQQQ, "    Zg %ld ions heat/cool: %.4e %.4e\n", gptr->DustZ,
			  gptr->FracPop*Heat1*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
			  gptr->FracPop*Cool1*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3 );
		}

		/* ============================================================================= */
		/* heating/cooling due to electrons */

		ion = -1;
		Stick = ( gptr->DustZ <= -1 ) ? gv.bin[nd]->StickElecNeg : gv.bin[nd]->StickElecPos;
		/* VE is mean (not RMS) electron velocity */
		/*ve = TePowers.sqrte*6.2124e5;*/
		ve = sqrt(8.*BOLTZMANN/PI/ELECTRON_MASS*phycon.te);

		/* electron arrival rate - eqn 29 again */
		CollisionRateElectr = Stick*dense.eden*ve;

		/* >>chng 00 jul 19, replace classical results with results including image potential
		 * to correct for polarization of the grain as charged particle approaches. */
		GrainScreen(ion,nd,nz,&eta,&xi);

		if( gptr->DustZ > gv.bin[nd]->LowestZg )
		{
			HeatCollisions = CollisionRateElectr*2.*BOLTZMANN*phycon.te*xi;
			/* this is change in energy due to charge acceleration within grain's potential 
			 * this term (perhaps) adds up to zero when summed over all charged particles */
			CoolPotential = CollisionRateElectr * (double)ion*gptr->PotSurfInc*eta*EN1RYD;
			/* >>chng 00 jul 05, this is term for energy released due to recombination, PvH */
			HeatRecombination = CollisionRateElectr * gptr->ThresSurfInc*eta*EN1RYD;
			HeatBounce = 0.;
			CoolBounce = 0.;
		}
		else
		{
			HeatCollisions = 0.;
			CoolPotential = 0.;
			HeatRecombination = 0.;
			/* >>chng 00 jul 05, add in terms for electrons that bounce off grain, PvH */
			/* >>chng 01 mar 09, remove these terms, their contribution is negligible, and replace
			 * them with similar terms that describe electrons that are captured by grains at Z_min,
			 * these electrons are not in a bound state and the grain will quickly autoionize, PvH */
			HeatBounce = CollisionRateElectr * 2.*BOLTZMANN*phycon.te*xi;
			/* >>chng 01 mar 14, replace (2kT_g - phi_g) term with -EA; for autoionizing states EA is
			 * usually higher than phi_g, so more energy is released back into the electron gas, PvH */ 
			CoolBounce = CollisionRateElectr *
				(-gptr->ThresSurfInc-gptr->PotSurfInc)*EN1RYD*eta;
			CoolBounce = MAX2(CoolBounce,0.);
		}
			
		/* >>chng 00 may 02, CoolPotential had not been included */
		/* >>chng 00 jul 05, HeatRecombination had not been included */
		HeatElectrons = HeatCollisions-CoolPotential+HeatRecombination+HeatBounce-CoolBounce;
		Heat1 += HeatElectrons;

		CoolElectrons = HeatCollisions+HeatBounce-CoolBounce;
		Cool1 += CoolElectrons;

		if( trace.lgTrace && trace.lgDustBug )
		{
			fprintf( ioQQQ, "    Zg %ld electrons heat/cool: %.4e %.4e\n", gptr->DustZ,
			  gptr->FracPop*HeatElectrons*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
			  gptr->FracPop*CoolElectrons*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3 );
		}

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

		/* calculate net heating rate in erg/H/s at standard depl
		 * include contributions for recombining electrons, autoionizing electrons
		 * and subtract thermionic emissions here since it is inverse process
		 *
		 * NB - in extreme conditions this rate may become 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. */
		gptr->HeatingRate2 = HeatElectrons*gv.bin[nd]->IntArea/4. -
			gv.bin[nd]->GrainCoolTherm*gv.bin[nd]->cnv_CM3_pH;

		/* >>chng 04 jan 25, moved inclusion into phitilde to qheat_init(), PvH */

		/* heating/cooling above is in erg/s/cm^2 -> multiply with projected grain area per cm^3 */
		/* GraC 0 is integral of dcheat, the total collisional heating of the grain */
		HeatTot += gptr->FracPop*Heat1;

		/* GrGC 0 total cooling of gas integrated */
		CoolTot += gptr->FracPop*Cool1;

		gv.bin[nd]->ChemEn += gptr->FracPop*ChemEn1;

		cscheat[nz] = (float)Heat1;
	}

	/* ============================================================================= */
	/* heating/cooling due to molecules */

	/* these rates do not depend on charge, hence they are outside of nz loop */

	/* sticking prob for H2 onto grain,
	 * estimated from accomodation coefficient defined at end of A.10 in BFM */
	WeightMol = 2.*dense.AtomicWeight[ipHYDROGEN];
	Accommodation = 2.*gv.bin[nd]->atomWeight*WeightMol/POW2(gv.bin[nd]->atomWeight+WeightMol);
	/* molecular hydrogen onto grains */
#ifndef IGNORE_GRAIN_ION_COLLISIONS
	/*CollisionRateMol = Accommodation*hmi.Hmolec[ipMH2g]* */
	CollisionRateMol = Accommodation*hmi.H2_total*
		sqrt(8.*BOLTZMANN/PI/ATOMIC_MASS_UNIT/WeightMol*phycon.te);
	/* >>chng 03 feb 12, added grain heating by H2 formation on the surface, PvH 
	 * >>refer	grain	H2 heat	Takahashi & Uehara, ApJ, 561, 843 */
	ipH2 = gv.which_H2distr[gv.bin[nd]->matType];
	/* this is rate in erg/cm^3/s */
	/* >>chng 04 may 26, changed dense.gas_phase[ipHYDROGEN] -> dense.xIonDense[ipHYDROGEN][0], PvH */
	gv.bin[nd]->ChemEnH2 = gv.bin[nd]->rate_h2_form_grains_used*dense.xIonDense[ipHYDROGEN][0]*
		H2_FORMATION_GRAIN_HEATING[ipH2]*EN1EV;
	/* convert to rate per cm^2 of projected grain surface area used here */
	gv.bin[nd]->ChemEnH2 /=	gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
#else
	CollisionRateMol = 0.;
	gv.bin[nd]->ChemEnH2 = 0.;
#endif

	/* now add in CO */
	WeightMol = dense.AtomicWeight[ipCARBON] + dense.AtomicWeight[ipOXYGEN];
	Accommodation = 2.*gv.bin[nd]->atomWeight*WeightMol/POW2(gv.bin[nd]->atomWeight+WeightMol);
#ifndef IGNORE_GRAIN_ION_COLLISIONS
	CollisionRateMol += Accommodation*co.hevmol[ipCO]*
		sqrt(8.*BOLTZMANN/PI/ATOMIC_MASS_UNIT/WeightMol*phycon.te);
#else
	CollisionRateMol = 0.;
#endif

	/* xi and eta are unity for neutrals and so ignored */
	HeatCollisions = CollisionRateMol * 2.*BOLTZMANN*phycon.te;
	CoolEmitted = CollisionRateMol * 2.*BOLTZMANN*gv.bin[nd]->tedust;

	HeatMolecules = HeatCollisions - CoolEmitted + gv.bin[nd]->ChemEnH2;
	HeatTot += HeatMolecules;

	/* >>chng 00 may 05, reversed sign of gas cooling contribution */
	CoolMolecules = HeatCollisions - CoolEmitted;
	CoolTot += CoolMolecules;

	gv.bin[nd]->RateUp = 0.;
	gv.bin[nd]->RateDn = 0.;
	HeatCor = 0.;
	for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
	{
		double d[4],Auger;
		double rate_dn = GrainElecRecomb1(nd,nz,&d[0],&d[1]);
		double rate_up = GrainElecEmis1(nd,nz,&d[0],&Auger,&d[1],&d[2],&d[3]);

		gv.bin[nd]->RateUp += gv.bin[nd]->chrg[nz]->FracPop*rate_up;
		gv.bin[nd]->RateDn += gv.bin[nd]->chrg[nz]->FracPop*rate_dn;

		/* >>chng 01 dec 19, subtract the Auger rate, it should not be used in heating correction */
		/* Auger/Compton-recoil processes are disabled for now */
		/*TODO	2   a self-consistent treatment for the heating by Auger electrons should be used */
		/*TODO	2   a self-consistent treatment for the heating by Compton recoil electrons should be used */
		HeatCor += (gv.bin[nd]->chrg[nz]->FracPop*(rate_up-Auger)*gv.bin[nd]->chrg[nz]->ThresSurf -
			    gv.bin[nd]->chrg[nz]->FracPop*rate_dn*gv.bin[nd]->chrg[nz]->ThresSurfInc +
			    gv.bin[nd]->chrg[nz]->FracPop*(rate_up-Auger)*gv.bin[nd]->chrg[nz]->PotSurf -
			    gv.bin[nd]->chrg[nz]->FracPop*rate_dn*gv.bin[nd]->chrg[nz]->PotSurfInc)*EN1RYD;
	}
	/* >>chng 01 nov 24, correct for imperfections in the n-charge state model,
	 * these corrections should add up to zero, but are actually small but non-zero, PvH */
	HeatTot += HeatCor;

	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "    molecules heat/cool: %.4e %.4e heatcor: %.4e\n",
			 HeatMolecules*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
			 CoolMolecules*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
			 HeatCor*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3 );
	}

	*dcheat = (float)(HeatTot*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3);
	*dccool = ( gv.lgDColOn ) ? (float)(CoolTot*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3) : 0.f;

	gv.bin[nd]->ChemEn *= gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
	gv.bin[nd]->ChemEnH2 *= gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;

	/* >>chng 04 feb 08, remember collisional heating for each charge state, PvH */
	for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
	{
		cscheat[nz] += (float)(HeatMolecules + HeatCor);
		cscheat[nz] *= (float)(gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3);
	}

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

	/* calculate heating rate in erg/H/s at standard depl
	 * include 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 01 nov 26, add in HeatCor as well, otherwise energy imbalance will result, PvH */
	gv.bin[nd]->HeatingRate1 = (HeatMolecules+HeatIons+HeatCor)*gv.bin[nd]->IntArea/4.;

	/* >>chng 04 jan 25, moved inclusion into phiTilde to qheat_init(), PvH */

#	ifdef DEBUG_FUN
	fputs( " <->GrainCollHeating()\n", debug_fp );
#	endif
	return;
}
#else
/* GrainCollHeating compute grains collisional heating cooling */
STATIC void GrainCollHeating(long int nd,
			     /*@out@*/ float *dcheat,
			     /*@out@*/ float *dccool,
			     /*@out@*/ float cscheat[]) /* cscheat[NCHU] */
{
	long int nz;
	double ChemEn,
	  ChemEn1,
	  ChemEnH2,
	  CoolTot,
	  CoolTot1,
	  HeatingRate2,
	  HeatIons,
	  HeatIons1,
	  HeatTot,
	  HeatTot1;

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

	HeatTot = 0.;
	CoolTot = 0.;
	HeatIons = 0.;
	ChemEn = 0.;

	for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
	{
		GrainCollHeating1(nd,nz,&HeatTot1,&CoolTot1,&HeatIons1,&ChemEn1,&HeatingRate2);

		HeatTot += gv.bin[nd]->chrg[nz]->FracPop*HeatTot1;
		CoolTot += gv.bin[nd]->chrg[nz]->FracPop*CoolTot1;
		HeatIons += gv.bin[nd]->chrg[nz]->FracPop*HeatIons1;
		ChemEn += gv.bin[nd]->chrg[nz]->FracPop*ChemEn1;

		/* >>chng 04 feb 08, remember collisional heating for each charge state, PvH */
		cscheat[nz] = (float)HeatTot1;
		gv.bin[nd]->chrg[nz]->HeatingRate2 = HeatingRate2;
	}

	GrainCollHeating2(nd,&HeatTot1,&CoolTot1,&ChemEnH2);

	HeatTot += HeatTot1;
	CoolTot += CoolTot1;

	*dcheat = (float)(HeatTot*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3);
	*dccool = ( gv.lgDColOn ) ? (float)(CoolTot*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3) : 0.f;

	for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
	{
		cscheat[nz] += (float)HeatTot1;
		cscheat[nz] *= (float)(gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3);
	}

	gv.bin[nd]->ChemEn = ChemEn * gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
	gv.bin[nd]->ChemEnH2 = ChemEnH2 * gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;

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

	/* calculate heating rate in erg/H/s at standard depl
	 * include 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 01 nov 26, add in HeatCor as well, otherwise energy imbalance will result, PvH */
	gv.bin[nd]->HeatingRate1 = (HeatTot1+HeatIons)*gv.bin[nd]->IntArea/4.;

	/* >>chng 04 jan 25, moved inclusion into phiTilde to qheat_init(), PvH */

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


/* GrainCollHeating compute grains collisional heating cooling dur to ions/electrons */
STATIC void GrainCollHeating1(long nd,
			      long nz,
			      /*@out@*/ double *HeatTot,
			      /*@out@*/ double *CoolTot,
			      /*@out@*/ double *HeatIons,
			      /*@out@*/ double *ChemEn,
			      /*@out@*/ double *HeatingRate2)
{
	long int ion,
	  nelem;
	double CollisionRateElectr, /* rate electrons strike grains */
	  CollisionRateIon,         /* rate ions strike grains */
	  CoolBounce,
	  CoolEmitted,
	  CoolElectrons,
	  CoolPotential,
	  CoolPotentialGas,
	  eta,
	  HeatBounce,
	  HeatCollisions,
	  HeatElectrons,
	  HeatRecombination, /* sum of abundances of ions times velocity times ionization potential times eta */
	  HeatChem,
	  Stick,
	  ve,
	  xi;

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


	/* >>chng 01 may 07, this routine now completely supports the hybrid grain
	 * charge model, and the average charge state is not used anywhere anymore, PvH */

	/* this subroutine evaluates the gas heating-cooling rate
	 * (erg cm^-3 s^-1) due to grain gas collisions.
	 * the net effect can be positive or negative,
	 * depending on whether the grains or gas are hotter
	 * the physics is described in 
	 * >>refer	grain	physics	Baldwin, Ferland, Martin et al., 1991, ApJ 374, 580 */

	ChargeBin *gptr = gv.bin[nd]->chrg[nz];

	/* HEAT1 will be rate collisions heat the grain
	 * COOL1 will be rate collisions cool the gas kinetics */
	double Heat1 = 0.;
	double Cool1 = 0.;
	double ChemEn1 = 0.;

	/* ============================================================================= */
	/* heating/cooling due to neutrals and positive ions */

	/* loop over all stages of ionization */
	for( ion=0; ion <= LIMELM; ion++ )
	{
		/* this is heating of grains due to recombination energy of species,
		 * and assumes that every ion is fully neutralized upon striking the grain surface.
		 * all radiation produced in the recombination process is absorbed within the grain
		 *
		 * ion=0 are neutrals, ion=1 are single ions, etc
		 * each population is weighted by the AVERAGE velocity
		 * */
		CollisionRateIon = 0.;
		CoolPotential = 0.;
		CoolPotentialGas = 0.;
		HeatRecombination = 0.;
		HeatChem = 0.;

		/* >>chng 00 jul 19, replace classical results with results including image potential
		 * to correct for polarization of the grain as charged particle approaches. */
		GrainScreen(ion,nd,nz,&eta,&xi);

		for( nelem=MAX2(0,ion-1); nelem < LIMELM; nelem++ )
		{
			if( dense.lgElmtOn[nelem] && dense.xIonDense[nelem][ion] > 0. )
			{
				double CollisionRateOne;

				/* >>chng 00 apr 05, use correct accomodation coefficient, by PvH
				 * the coefficient is defined at the end of appendix A.10 of BFM
				 * assume ion sticking prob is unity */
#if defined( IGNORE_GRAIN_ION_COLLISIONS )
				Stick = 0.;
#elif defined( WD_TEST2 )
				Stick = ( ion == gptr->RecomZ0[nelem][ion] ) ? 0. : STICK_ION;
#else
				Stick = ( ion == gptr->RecomZ0[nelem][ion] ) ?
					gv.bin[nd]->AccomCoef[nelem] : STICK_ION;
#endif
				/* this is rate with which charged ion strikes grain */
				/* >>chng 00 may 02, this had left 2./SQRTPI off */
				/* >>chng 00 may 05, use average speed instead of 2./SQRTPI*Doppler, PvH */
				CollisionRateOne = Stick*dense.xIonDense[nelem][ion]*DoppVel.AveVel[nelem];
				CollisionRateIon += CollisionRateOne;
				/* >>chng 01 nov 26, use PotSurfInc when appropriate:
				 * the values for the surface potential used here make it
				 * consistent with the rest of the code and preserve energy.
				 * NOTE: For incoming particles one should use PotSurfInc with
				 * Schottky effect for positive ion, for outgoing particles
				 * one should use PotSurf for Zg+ion-Z_0-1 (-1 because PotSurf
				 * assumes electron going out), these corrections are small
				 * and will be neglected for now, PvH */
				if( ion >= gptr->RecomZ0[nelem][ion] )
				{
					CoolPotential += CollisionRateOne * (double)ion * gptr->PotSurf;
					CoolPotentialGas += CollisionRateOne *
						(double)gptr->RecomZ0[nelem][ion] * gptr->PotSurf;
				}
				else
				{
					CoolPotential += CollisionRateOne * (double)ion * gptr->PotSurfInc;
					CoolPotentialGas += CollisionRateOne *
						(double)gptr->RecomZ0[nelem][ion] * gptr->PotSurfInc;
				}
				/* this is sum of all energy liberated as ion recombines to Z0 in grain */
				/* >>chng 00 jul 05, subtract energy needed to get 
				 * electron out of grain potential well, PvH */
				/* >>chng 01 may 09, chemical energy now calculated in GrainIonColl, PvH */
				HeatRecombination += CollisionRateOne * gptr->RecomEn[nelem][ion];
				HeatChem += CollisionRateOne * gptr->ChemEn[nelem][ion];
			}
		}

		/* >>chng 00 may 01, boltzmann factor had multiplied all of factor instead
		 * of only first and last term.  pvh */

		/* equation 29 from Balwin et al 91 */
		/* this is direct collision rate, 2kT * xi, first term in eq 29 */
		HeatCollisions = CollisionRateIon * 2.*BOLTZMANN*phycon.te*xi;
		/* this is change in energy due to charge acceleration within grain's potential 
		 * this is exactly balanced by deceleration of incoming electrons and accelaration
		 * of outgoing photo-electrons and thermionic emissions; all these terms should
		 * add up to zero (total charge of grain should remain constant) */
		CoolPotential *= eta*EN1RYD;
		CoolPotentialGas *= eta*EN1RYD;
		/* this is recombination energy released within grain */
		HeatRecombination *= eta*EN1RYD;
		HeatChem *= eta*EN1RYD;
		/* energy carried away by neutrals after recombination, so a cooling term */
		CoolEmitted = CollisionRateIon * 2.*BOLTZMANN*gv.bin[nd]->tedust*eta;

		/* total GraC 0 in the emission line output */
		Heat1 += HeatCollisions - CoolPotential + HeatRecombination - CoolEmitted;

		/* rate kinetic energy lost from gas - gas cooling - eq 32 in BFM */
		/* this GrGC 0 in the main output */
		/* >>chng 00 may 05, reversed sign of gas cooling contribution */
		Cool1 += HeatCollisions - CoolEmitted - CoolPotentialGas;

		ChemEn1 += HeatChem;
	}

	/* remember grain heating by ion collisions for quantum heating treatment */
	*HeatIons = Heat1;

	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "    Zg %ld ions heat/cool: %.4e %.4e\n", gptr->DustZ,
			 gptr->FracPop*Heat1*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
			 gptr->FracPop*Cool1*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3 );
	}

	/* ============================================================================= */
	/* heating/cooling due to electrons */

	ion = -1;
	Stick = ( gptr->DustZ <= -1 ) ? gv.bin[nd]->StickElecNeg : gv.bin[nd]->StickElecPos;
	/* VE is mean (not RMS) electron velocity */
	/*ve = TePowers.sqrte*6.2124e5;*/
	ve = sqrt(8.*BOLTZMANN/PI/ELECTRON_MASS*phycon.te);

	/* electron arrival rate - eqn 29 again */
	CollisionRateElectr = Stick*dense.eden*ve;

	/* >>chng 00 jul 19, replace classical results with results including image potential
	 * to correct for polarization of the grain as charged particle approaches. */
	GrainScreen(ion,nd,nz,&eta,&xi);

	if( gptr->DustZ > gv.bin[nd]->LowestZg )
	{
		HeatCollisions = CollisionRateElectr*2.*BOLTZMANN*phycon.te*xi;
		/* this is change in energy due to charge acceleration within grain's potential 
		 * this term (perhaps) adds up to zero when summed over all charged particles */
		CoolPotential = CollisionRateElectr * (double)ion*gptr->PotSurfInc*eta*EN1RYD;
		/* >>chng 00 jul 05, this is term for energy released due to recombination, PvH */
		HeatRecombination = CollisionRateElectr * gptr->ThresSurfInc*eta*EN1RYD;
		HeatBounce = 0.;
		CoolBounce = 0.;
	}
	else
	{
		HeatCollisions = 0.;
		CoolPotential = 0.;
		HeatRecombination = 0.;
		/* >>chng 00 jul 05, add in terms for electrons that bounce off grain, PvH */
		/* >>chng 01 mar 09, remove these terms, their contribution is negligible, and replace
		 * them with similar terms that describe electrons that are captured by grains at Z_min,
		 * these electrons are not in a bound state and the grain will quickly autoionize, PvH */
		HeatBounce = CollisionRateElectr * 2.*BOLTZMANN*phycon.te*xi;
		/* >>chng 01 mar 14, replace (2kT_g - phi_g) term with -EA; for autoionizing states EA is
		 * usually higher than phi_g, so more energy is released back into the electron gas, PvH */ 
		CoolBounce = CollisionRateElectr *
			(-gptr->ThresSurfInc-gptr->PotSurfInc)*EN1RYD*eta;
		CoolBounce = MAX2(CoolBounce,0.);
	}

	/* >>chng 00 may 02, CoolPotential had not been included */
	/* >>chng 00 jul 05, HeatRecombination had not been included */
	HeatElectrons = HeatCollisions-CoolPotential+HeatRecombination+HeatBounce-CoolBounce;
	Heat1 += HeatElectrons;

	CoolElectrons = HeatCollisions+HeatBounce-CoolBounce;
	Cool1 += CoolElectrons;

	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "    Zg %ld electrons heat/cool: %.4e %.4e\n", gptr->DustZ,
			 gptr->FracPop*HeatElectrons*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
			 gptr->FracPop*CoolElectrons*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3 );
	}

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

	/* calculate net heating rate in erg/H/s at standard depl
	 * include contributions for recombining electrons, autoionizing electrons
	 * and subtract thermionic emissions here since it is inverse process
	 *
	 * NB - in extreme conditions this rate may become 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. */
	*HeatingRate2 = HeatElectrons*gv.bin[nd]->IntArea/4. - gv.bin[nd]->GrainCoolTherm*gv.bin[nd]->cnv_CM3_pH;

	/* >>chng 04 jan 25, moved inclusion into phitilde to qheat_init(), PvH */

	/* heating/cooling is in erg/s/cm^2 -> needs to be multiplied with projected grain area per cm^3 */
	/* GraC 0 is integral of dcheat, the total collisional heating of the grain */
	*HeatTot = Heat1;

	/* GrGC 0 total cooling of gas integrated */
	*CoolTot = Cool1;

	*ChemEn = ChemEn1;

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


/* GrainCollHeating2 compute grains collisional heating/cooling due to molecules */
STATIC void GrainCollHeating2(long int nd,
			      /*@out@*/ double *HeatTot,
			      /*@out@*/ double *CoolTot,
			      /*@out@*/ double *ChemEnH2)
{
	long nz;
	H2_type ipH2;
	double Accommodation,
	  CollisionRateMol,         /* rate molecules strike grains */
	  CoolEmitted,
	  CoolMolecules,
	  HeatCollisions,
	  HeatMolecules,
	  HeatCor,
	  WeightMol;

	/* energy deposited into grain by formation of a single H2 molecule, in eV,
	 * >>refer  Takahashi J., Uehara H., 2001, ApJ, 561, 843 */
	const double H2_FORMATION_GRAIN_HEATING[H2_TOP] = { 0.20, 0.4, 1.72 };

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

	/* ============================================================================= */
	/* heating/cooling due to molecules */

	/* these rates do not depend on charge, hence they are outside of nz loop */

	/* sticking prob for H2 onto grain,
	 * estimated from accomodation coefficient defined at end of A.10 in BFM */
	WeightMol = 2.*dense.AtomicWeight[ipHYDROGEN];
	Accommodation = 2.*gv.bin[nd]->atomWeight*WeightMol/POW2(gv.bin[nd]->atomWeight+WeightMol);
	/* molecular hydrogen onto grains */
#ifndef IGNORE_GRAIN_ION_COLLISIONS
	/*CollisionRateMol = Accommodation*hmi.Hmolec[ipMH2g]* */
	CollisionRateMol = Accommodation*hmi.H2_total*
		sqrt(8.*BOLTZMANN/PI/ATOMIC_MASS_UNIT/WeightMol*phycon.te);
	/* >>chng 03 feb 12, added grain heating by H2 formation on the surface, PvH 
	 * >>refer	grain	H2 heat	Takahashi & Uehara, ApJ, 561, 843 */
	ipH2 = gv.which_H2distr[gv.bin[nd]->matType];
	/* this is rate in erg/cm^3/s */
	/* >>chng 04 may 26, changed dense.gas_phase[ipHYDROGEN] -> dense.xIonDense[ipHYDROGEN][0], PvH */
	*ChemEnH2 = gv.bin[nd]->rate_h2_form_grains_used*dense.xIonDense[ipHYDROGEN][0]*
		H2_FORMATION_GRAIN_HEATING[ipH2]*EN1EV;
	/* convert to rate per cm^2 of projected grain surface area used here */
	*ChemEnH2 /= gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
#else
	CollisionRateMol = 0.;
	*ChemEnH2 = 0.;
#endif

	/* now add in CO */
	WeightMol = dense.AtomicWeight[ipCARBON] + dense.AtomicWeight[ipOXYGEN];
	Accommodation = 2.*gv.bin[nd]->atomWeight*WeightMol/POW2(gv.bin[nd]->atomWeight+WeightMol);
#ifndef IGNORE_GRAIN_ION_COLLISIONS
	CollisionRateMol += Accommodation*co.hevmol[ipCO]*
		sqrt(8.*BOLTZMANN/PI/ATOMIC_MASS_UNIT/WeightMol*phycon.te);
#else
	CollisionRateMol = 0.;
#endif

	/* xi and eta are unity for neutrals and so ignored */
	HeatCollisions = CollisionRateMol * 2.*BOLTZMANN*phycon.te;
	CoolEmitted = CollisionRateMol * 2.*BOLTZMANN*gv.bin[nd]->tedust;

	HeatMolecules = HeatCollisions - CoolEmitted + *ChemEnH2;
	*HeatTot = HeatMolecules;

	/* >>chng 00 may 05, reversed sign of gas cooling contribution */
	CoolMolecules = HeatCollisions - CoolEmitted;
	*CoolTot = CoolMolecules;

	gv.bin[nd]->RateUp = 0.;
	gv.bin[nd]->RateDn = 0.;
	HeatCor = 0.;
	for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
	{
		double d[4],Auger;
		double rate_dn = GrainElecRecomb1(nd,nz,&d[0],&d[1]);
		double rate_up = GrainElecEmis1(nd,nz,&d[0],&Auger,&d[1],&d[2],&d[3]);

		gv.bin[nd]->RateUp += gv.bin[nd]->chrg[nz]->FracPop*rate_up;
		gv.bin[nd]->RateDn += gv.bin[nd]->chrg[nz]->FracPop*rate_dn;

		/* >>chng 01 dec 19, subtract the Auger rate, it should not be used in heating correction */
		/* Auger/Compton-recoil processes are disabled for now */
		/*TODO	2   a self-consistent treatment for the heating by Auger electrons should be used */
		/*TODO	2   a self-consistent treatment for the heating by Compton recoil electrons should be used */
		HeatCor += (gv.bin[nd]->chrg[nz]->FracPop*(rate_up-Auger)*gv.bin[nd]->chrg[nz]->ThresSurf -
			    gv.bin[nd]->chrg[nz]->FracPop*rate_dn*gv.bin[nd]->chrg[nz]->ThresSurfInc +
			    gv.bin[nd]->chrg[nz]->FracPop*(rate_up-Auger)*gv.bin[nd]->chrg[nz]->PotSurf -
			    gv.bin[nd]->chrg[nz]->FracPop*rate_dn*gv.bin[nd]->chrg[nz]->PotSurfInc)*EN1RYD;
	}
	/* >>chng 01 nov 24, correct for imperfections in the n-charge state model,
	 * these corrections should add up to zero, but are actually small but non-zero, PvH */
	*HeatTot += HeatCor;

	if( trace.lgTrace && trace.lgDustBug )
	{
		fprintf( ioQQQ, "    molecules heat/cool: %.4e %.4e heatcor: %.4e\n",
			 HeatMolecules*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
			 CoolMolecules*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
			 HeatCor*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3 );
	}

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


/* GrainDrift computes grains drift velocity */
void GrainDrift(void)
{
	long int i, 
	  loop, 
	  nd;
	float *help;
	double alam, 
	  corr, 
	  dmomen, 
	  fac, 
	  fdrag, 
	  g0, 
	  g2, 
	  phi2lm, 
	  psi, 
	  rdust, 
	  si, 
	  vdold, 
	  volmom;

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

	/* >>chng 04 jan 31, use help array to store intermediate results, PvH */
	if( ( help = (float*)MALLOC((size_t)((unsigned)rfield.nflux*sizeof(float))) ) == NULL ) 
		BadMalloc();
	for( i=0; i < rfield.nflux; i++ )
	{
		help[i] = (rfield.flux[i]+rfield.ConInterOut[i]+rfield.outlin[i]+rfield.outlin_noplot[i])*
			rfield.anu[i];
	}

	for( nd=0; nd < gv.nBin; nd++ )
	{
		/* find momentum absorbed by grain */
		dmomen = 0.;
		for( i=0; i < rfield.nflux; i++ )
		{
			/* >>chng 02 dec 30, separated scattering cross section and asymmetry factor, PvH */
			dmomen += help[i]*(gv.bin[nd]->dstab1[i] + gv.bin[nd]->pure_sc1[i]*gv.bin[nd]->asym[i]);
		}
		ASSERT( dmomen >= 0. );
		dmomen *= EN1RYD*4./gv.bin[nd]->IntArea;

		/* now find force on grain, and drift velocity */
		fac = 2*BOLTZMANN*phycon.te;

		/* now PSI defined by 
		 * >>refer	grain	physics	Draine and Salpeter 79 Ap.J. 231, 77 (1979) */
		psi = gv.bin[nd]->dstpot*TE1RYD/phycon.te;
		if( psi > 0. )
		{
			rdust = 1.e-6;
			alam = log(20.702/rdust/psi*phycon.sqrte/dense.SqrtEden);
		}
		else
		{
			alam = 0.;
		}

		phi2lm = POW2(psi)*alam;
		corr = 2.;
		/* >>chng 04 jan 31, increased loop limit 10 -> 50, precision -> 0.001, PvH */
		for( loop = 0; loop < 50 && fabs(corr-1.) > 0.001; loop++ )
		{
			vdold = gv.bin[nd]->DustDftVel;

			/* interactions with protons */
			si = gv.bin[nd]->DustDftVel/phycon.sqrte*7.755e-5;
			g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
			g2 = si/(1.329 + POW3(si));

			/* drag force due to protons, both linear and square in velocity
			 * equation 4 from D+S Ap.J. 231, p77. */
			fdrag = fac*dense.xIonDense[ipHYDROGEN][1]*(g0 + phi2lm*g2);

			/* drag force due to interactions with electrons */
			si = gv.bin[nd]->DustDftVel/phycon.sqrte*1.816e-6;
			g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
			g2 = si/(1.329 + POW3(si));
			fdrag += fac*dense.eden*(g0 + phi2lm*g2);

			/* drag force due to collisions with hydrogen and helium atoms */
			si = gv.bin[nd]->DustDftVel/phycon.sqrte*7.755e-5;
			g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
			fdrag += fac*(dense.xIonDense[ipHYDROGEN][0] + 1.1*dense.xIonDense[ipHELIUM][0])*g0;

			/* drag force due to interactions with helium ions */
			si = gv.bin[nd]->DustDftVel/phycon.sqrte*1.551e-4;
			g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
			g2 = si/(1.329 + POW3(si));
			fdrag += fac*dense.xIonDense[ipHELIUM][1]*(g0 + phi2lm*g2);

			/* this term does not work
			 *  2      HEIII*(G0+4.*PSI**2*(ALAM-0.693)*G2) )
			 * this is total momentum absorbed by dust per unit vol */
			volmom = dmomen/SPEEDLIGHT;

			if( fdrag > 0. )
			{
				corr = sqrt(volmom/fdrag);
				gv.bin[nd]->DustDftVel = (float)(vdold*corr);
			}
			else
			{
				corr = 1.;
				negdrg.lgNegGrnDrg = TRUE;
				gv.bin[nd]->DustDftVel = 0.;
			}

			if( trace.lgTrace && trace.lgDustBug )
			{
				fprintf( ioQQQ, "     %2ld new drift velocity:%10.2e momentum absorbed:%10.2e\n", 
				  loop, gv.bin[nd]->DustDftVel, volmom );
			}
		}
	}

	free( help );

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

/*GrnVryDpth set grains abundance as a function of depth into cloud*/
STATIC double GrnVryDpth(

/* nd is the number of the grain bin. The values are listed in the Cloudy output,
 * under "Average Grain Properties", and can easily be obtained by doing a trial
 * run without varying the grain abundance and setting stop zone to 1 */

	long int nd)
{
	double GrnVryDpth_v;

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

	/* set grains abundance as a function of depth into cloud
	 * NB most quantities are undefined for first calls to this sub */
	/* nd is the index of the grain bin.  This routine must return
	 * a scale factor for the grain abundance at this position. */

	if( gv.bin[nd]->lgDustVary )
	{
		/* the scale factor is the hydrogen atomic fraction, small when gas is ionized or molecular
		 * and unity when atomic.  This function is observed for PAHs across the Orion Bar, the
		 * PAHs are strong near the ionization front and weak in the ionized and molecular gas */
		/* >>chng 04 sep 28, propto atomic graction */
		GrnVryDpth_v = dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN];

		ASSERT( GrnVryDpth_v > 0. && GrnVryDpth_v <= 1.000001 );
	}
	else
	{
		GrnVryDpth_v = 1.;
	}


#	ifdef DEBUG_FUN
	fputs( " <->GrnVryDpth()\n", debug_fp );
#	endif
	return GrnVryDpth_v;
}
/*lint +e662 creation of out of bounds pointer */
