/* 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 */
#include <ctype.h>
#include "cddefines.h"
#include "elementnames.h"
#include "physconst.h"
#include "path.h"
#include "dense.h"
#include "called.h"
#include "version.h"
#include "grainvar.h"
#include "rfield.h"
#include "atmdat.h"
#include "math_complex.h"
#include "grains.h"

/*=======================================================*
 *
 * Mie code for spherical grains.
 *
 * Calculates <pi*a^2*Q_abs>, <pi*a^2*Q_sct>, and (1-<g>)
 * for arbitrary grain species and size distribution.
 *
 * This code is derived from the program cmieuvx.f
 *
 * Written by: P.G. Martin (CITA), based on the code described in
 * >>refer Hansen, J. E., Travis, L. D. 1974, Space Sci. Rev., 16, 527
 *
 * Adapted for Cloudy by Peter A.M. van Hoof (University of Kentucky,
 *   Canadian Institute for Theoretical Astrophysics,
 *   Queen's University of Belfast)
 *
 *=======================================================*/


/* these are the magic numbers for the .rfi, .szd, .opc, and .mix files
 * the first digit is file type, the rest is date (YYMMDD) */
const long MAGIC_RFI = 1030103L;
const long MAGIC_SZD = 2010403L;
const long MAGIC_OPC = 3030307L;
const long MAGIC_MIX = 4030103L;

#define NINT(X) ((long)((X) < 0. ? (X)-0.5 : (X)+0.5))
#define SIGN3(X) ((X) < 0. ? -1L : ((X) == 0. ? 0L : 1L))
#define POW4(X) (POW2(X)*POW2(X))

/* >>chng 02 may 28, by Ryan, moved struct complex to cddefines.h to make it available to entire code. */

/* these are the absolute smallest and largest grain sizes we will
 * consider (in micron). the lower limit gives a grain with on the
 * order of one atom in it, so it is physically motivated. the upper
 * limit comes from the series expansions used in the mie theory,
 * they will have increasingly more problems converging for larger
 * grains, so this limit is numerically motivated */
const double SMALLEST_GRAIN = 0.0001*(1.-10.*DBL_EPSILON);
const double LARGEST_GRAIN = 10.*(1.+10.*DBL_EPSILON);

/* maximum no. of parameters for grain size distribution */
#define NSD 7

/* these are the indices into the parameter array a[NSD],
 * NB NB -- the numbers defined below should range from 0 to NSD-1 */
const int ipSize  = 0; /* single size */
const int ipBLo   = 0; /* lower bound */
const int ipBHi   = 1; /* upper bound */
const int ipExp   = 2; /* exponent for powerlaw */
const int ipBeta  = 3; /* beta parameter for powerlaw */
const int ipSLo   = 4; /* scale size for lower exp. cutoff */
const int ipSHi   = 5; /* scale size for upper exp. cutoff */
const int ipAlpha = 6; /* alpha parameter for exp. cutoff */
const int ipGCen  = 2; /* center of gaussian distribution */
const int ipGSig  = 3; /* 1-sigma width of gaussian distribution */

/* these are the types of refractive index files we recognize */
typedef enum {
	RFI_TABLE, OPC_TABLE, OPC_GREY, OPC_PAH1
} rfi_type;

/* these are the types of EMT's we recognize */
typedef enum {
	FARAFONOV00, STOGNIENKO95, BRUGGEMAN35
} emt_type;

/* these are all the size distribution cases we support */
typedef enum {
	SD_ILLEGAL, SD_SINGLE_SIZE, SD_POWERLAW, SD_EXP_CUTOFF1, SD_EXP_CUTOFF2,
	SD_EXP_CUTOFF3, SD_LOG_NORMAL, SD_LIN_NORMAL, SD_TABLE
} sd_type;

typedef struct {
	double a[NSD];    /* parameters for size distribution */
	double lim[2];    /* holds lower and upper size limit for entire distribution */
	double clim[2];   /* holds lower and upper size limit for current bin */
	double *xx;       /* xx[nn]: abcissas for Gauss quadrature on [-1,1] */
	double *aa;       /* aa[nn]: weights for Gauss quadrature on [-1,1] */
	double *rr;       /* rr[nn]: abcissas for Gauss quadrature */
	double *ww;       /* ww[nn]: weights for Gauss quadrature */
	double unity;     /* normalization for integrals over size distribution */
	double unity_bin; /* normalization for integrals over size distribution */
	double cSize;     /* the size currently being used for the calculations */
	double radius;    /* average grain radius for current bin <a> */
	double area;      /* average grain surface area for current bin <4pi*a^2> */
	double vol;       /* average grain volume for current bin <4pi/3*a^3> */
	double *ln_a;     /* ln(a)[npts]: log of grain radii for user-supplied size distr */
	double *ln_a4dNda;/* ln(a^4*dN/da)[npts]: log of user-supplied size distr */
	sd_type sdCase;   /* SD_SINGLE_SIZE, SD_POWERLAW, ... */
	long int magic;   /* magic number */
	long int cPart;   /* current partition no. for size distribution */
	long int nPart;   /* total no. of partitions for size distribution */
	long int nmul;    /* multiplier for obtaining no. of abscissas in gaussian quadrature */
	long int nn;      /* no. of abscissas used in gaussian quadrature (per partition) */
	long int npts;    /* no. of points in user-supplied size distr */
	int lgLogScale;   /* use logarithmic mesh for integration over size ? */
} sd_data;

/* maximum no. of principal axes for crystalline grains */
#define NAX 3
#define NDAT 4

typedef struct {
	double *wavlen[NAX];    /* wavelength grid for rfi for all axes (micron) */
	complex *n[NAX];        /* refractive index n for all axes */
	double *nr1[NAX];       /* re(n)-1 for all axes */
	double *opcAnu;         /* energies for data points in OPC_TABLE file */
	double *opcData[NDAT];  /* data values from OPC_TABLE file */
	double wt[NAX];         /* relative weight of each axis */
	double abun;            /* abundance of grain molecule rel. to hydrogen for max depletion */
	double depl;            /* depletion efficiency */
	double elmAbun[LIMELM]; /* abundances of constituent elements rel. to hydrogen */
	double mol_weight;      /* molecular weight of grain molecule (amu) */
	double atom_weight;     /* molecular weight of grain molecule per atom (amu) */
	double rho;             /* specific weight (g/cm^3) */
	double norm;            /* number of protons in plasma per average grain */
	double work;            /* work function (Ryd) */
	double bandgap;         /* gap between valence and conduction band (Ryd) */
	double therm_eff;       /* efficiency of thermionic emission, between 0 and 1 */
	double subl_temp;       /* sublimation temperature (K) */
	long int magic;         /* magic number */
	long int cAxis;         /* number of axis currently being used */
	long int nAxes;         /* no. of principal axes for this grain */
	long int ndata[NAX];    /* no. of wavelength points for each axis */
	long int nOpcCols;      /* no. of data columns in OPC_TABLE file */
	long int nOpcData;      /* no. of data points in OPC_TABLE file */
	mat_type matType;       /* material type, determines enthalpy function, etc. */
	rfi_type rfiType;       /* type of data in rfi file: rfi table, grey grain, pah, etc. */
} grain_data;

/* used in mie_read_rfi, mie_read_opc, mie_write_opc */
/* >>chng 01 oct 23, to DATA_PATH_LENGTH_2 the largest possible path */
/*#define FILENAME_PATH_LENGTH_2 140*/

#define WORDLEN 5

/* maximum size for grain type labels */
#define LABELSUB1 3
#define LABELSUB2 5
#define LABELSIZE (LABELSUB1+LABELSUB2+4)

/* this is the number of data points used to set up table of optical constants for a mixed grain */
const long MIX_TABLE_SIZE = 2000L;

static void mie_auxiliary(/*@partial@*/sd_data*,/*@in@*/const char*);
static void mie_integrate(/*@partial@*/sd_data*,double,double,/*@out@*/double*,int);
static void mie_cs_size_distr(double,/*@in@*/sd_data*,/*@in@*/grain_data*,
			      void(*)(double,/*@in@*/sd_data*,/*@in@*/grain_data*,/*@out@*/double*,
				      /*@out@*/double*,/*@out@*/double*,/*@out@*/int*),
			      /*@out@*/double*,/*@out@*/double*,/*@out@*/double*,/*@out@*/int*);
static void mie_cs(double,/*@in@*/sd_data*,/*@in@*/grain_data*,/*@out@*/double*,/*@out@*/double*,
		   /*@out@*/double*,/*@out@*/int*);
static void pah1_fun(double,/*@in@*/sd_data*,/*@in@*/grain_data*,/*@out@*/double*,/*@out@*/double*,
		     /*@out@*/double*,/*@out@*/int*);
static void tbl_fun(double,/*@in@*/sd_data*,/*@in@*/grain_data*,/*@out@*/double*,/*@out@*/double*,
		    /*@out@*/double*,/*@out@*/int*);
static double size_distr(double,/*@in@*/sd_data*);
static double search_limit(double,double,double,sd_data);
static void mie_calc_ial(/*@in@*/grain_data*,long,/*@out@*/double[],/*@in@*/const char*,/*@in@*/int*);
static void mie_repair(/*@in@*/const char*,long,int,int,/*@in@*/float[],double[],/*@in@*/int[],
		       int,/*@in@*/int*);
static double mie_find_slope(/*@in@*/const float[],/*@in@*/const double[],/*@in@*/const int[],
			     long,long,int,int,/*@in@*/int*);
static void mie_read_rfi(/*@in@*/const char*,/*@out@*/grain_data*);
static void mie_read_mix(/*@in@*/const char*,/*@out@*/grain_data*);
static void init_eps(double,long,/*@in@*/grain_data[],/*@out@*/complex[]);
static complex cnewton(void(*)(complex,double[],complex[],long,complex*,double*,double*),
		       double[],complex[],long,complex,double);
static void Stognienko(complex,double[],complex[],long,complex*,double*,double*);
static void Bruggeman(complex,double[],complex[],long,complex*,double*,double*);
static void mie_read_szd(/*@in@*/const char*,/*@out@*/sd_data*);
static void mie_read_long(/*@in@*/const char*,/*@in@*/const char[],/*@out@*/long int*,int,long int);
static void mie_read_float(/*@in@*/const char*,/*@in@*/const char[],/*@out@*/float*,int,long int);
static void mie_read_double(/*@in@*/const char*,/*@in@*/const char[],/*@out@*/double*,int,long int);
static void mie_read_form(/*@in@*/const char*,/*@out@*/double[],/*@out@*/double*,/*@out@*/double*);
static void mie_write_form(/*@in@*/const double[],/*@out@*/char[]);
static void mie_read_word(/*@in@*/const char[],/*@out@*/char[],long,int);
static void mie_next_data(/*@in@*/const char*,/*@in@*/FILE*,/*@out@*/char*,/*@in@*/long int*);
static void mie_next_line(/*@in@*/const char*,/*@in@*/FILE*,/*@out@*/char*,/*@in@*/long int*);

/*=======================================================*/
/* the following five routines are the core of the Mie code supplied by Peter Martin */

static void sinpar(double,double,double,/*@out@*/double*,/*@out@*/double*,/*@out@*/double*,
		   /*@out@*/double*,/*@out@*/double*,/*@out@*/long*);
static void anomal(double,/*@out@*/double*,/*@out@*/double*,/*@out@*/double*,/*@out@*/double*,double,double);
static void bigk(complex,/*@out@*/complex*);
static void ritodf(double,double,/*@out@*/double*,/*@out@*/double*);
static void dftori(/*@out@*/double*,/*@out@*/double*,double,double);

/* >>chng 02 may 28, by Ryan, moved complex functions to math.h */

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

void mie_write_opc(/*@in@*/ const char *rfi_file,
		   /*@in@*/ const char *szd_file,
		   long int nbin)
{
	int Error = 0,
		/* >>chng 01 aug 01, MALLOC the space, no more ncell */
	  *ErrorIndex/*[NC_ELL]*/,
	  lgErr,
	  lgErrorOccurred,
	  lgWarning;
	long int i,
	  j,
	  nelem,
	  p;
	double **acs_abs,
	  **acs_sct,
	  **a1g,
	  *inv_att_len,
	  cosb,
	  cs_abs,
	  cs_sct,
	  volfrac,
	  volnorm,
	  wavlen;
	char chGrainLabel[LABELSIZE+1],
	  ext[3],
	  chFile[FILENAME_PATH_LENGTH_2],
	  chFile2[FILENAME_PATH_LENGTH_2],
	  hlp1[LABELSUB1+2],
	  hlp2[LABELSUB2+2],
	  *str,
	  string[FILENAME_PATH_LENGTH_2];
	sd_data sd;
	grain_data gd,
	  gd2;
	time_t timer;
	FILE *fdes;
	

	/* no. of logarithmic intervals in table printout of size distribution function */
	const long NPTS_TABLE = 100L;

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

	/* >>chng 01 aug 22, MALLOC this space */
	if( (ErrorIndex = (int*)MALLOC(sizeof(int)*(unsigned)rfield.nupper) ) == NULL )
		BadMalloc();

	gd.nAxes = 0;
	for( j=0; j < NAX; j++ ) 
	{
		gd.wavlen[j] = NULL;
		gd.n[j] = NULL;
		gd.nr1[j] = NULL;
	}
	gd.nOpcCols = 0;
	gd.opcAnu = NULL;
	for( j=0; j < NDAT; j++ )
	{
		gd.opcData[j] = NULL;
	}
	gd2.nAxes = 0;
	for( j=0; j < NAX; j++ ) 
	{
		gd2.wavlen[j] = NULL;
		gd2.n[j] = NULL;
		gd2.nr1[j] = NULL;
	}
	gd2.nOpcCols = 0;
	gd2.opcAnu = NULL;
	for( j=0; j < NDAT; j++ )
	{
		gd2.opcData[j] = NULL;
	}
	sd.ln_a = NULL;
	sd.ln_a4dNda = NULL;

	mie_read_szd( szd_file , &sd );

	sd.nPart = ( sd.sdCase == SD_SINGLE_SIZE ) ? 1 : nbin;
	if( sd.nPart <= 0 || sd.nPart >= 100 ) 
	{
		fprintf( ioQQQ, " Illegal number of size distribution bins: %ld\n",sd.nPart );
		fprintf( ioQQQ, " The number should be between 1 and 99.\n" );
		puts( "[Stop in mie_write_opc]" );
		cdEXIT(EXIT_FAILURE);
	}
	sd.lgLogScale = TRUE;

	str = strstr(rfi_file,".");
	if( str != NULL && lgMatch( str, ".rfi" ) )
	{
		mie_read_rfi( rfi_file , &gd );
	}
	else if( str != NULL && lgMatch( str, ".mix" ) )
	{
		mie_read_mix( rfi_file , &gd );
	}
	else
	{
		fprintf( ioQQQ, " Refractive index file name %s has wrong extention\n", rfi_file );
		fprintf( ioQQQ, " It should have extention .rfi or .mix.\n" );
		puts( "[Stop in mie_write_opc]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( gd.rfiType == OPC_TABLE && sd.nPart > 1 )
	{
		fprintf( ioQQQ, " Illegal number of size distribution bins: %ld\n",sd.nPart );
		fprintf( ioQQQ, " The number should always be 1 for OPC_TABLE files.\n" );
		puts( "[Stop in mie_write_opc]" );
		cdEXIT(EXIT_FAILURE);
	}
	if( gd.rho <= 0. )
	{
		fprintf( ioQQQ, " Illegal value for the specific density: %.4e\n",gd.rho );
		puts( "[Stop in mie_write_opc]" );
		cdEXIT(EXIT_FAILURE);
	}
	if( gd.mol_weight <= 0. )
	{
		fprintf( ioQQQ, " Illegal value for the molecular weight: %.4e\n",gd.mol_weight );
		puts( "[Stop in mie_write_opc]" );
		cdEXIT(EXIT_FAILURE);
	}

	lgWarning = FALSE;

	/* generate output file name from input file names */
	strcpy(chFile,rfi_file);
	str = strstr(chFile,".");

	if( str != NULL )
		*str = '\0';

	strcpy(chFile2,szd_file);
	str = strstr(chFile2,".");

	if( str != NULL )
		*str = '\0';

	if( sd.sdCase != SD_SINGLE_SIZE ) 
	{
		sprintf(ext,"%02ld",nbin);
		strcat(strcat(strcat(strcat(strcat(chFile,"_"),chFile2),"_"),ext),".opc");
	}
	else 
	{
		strcat(strcat(strcat(chFile,"_"),chFile2),".opc");
	}

	fprintf( ioQQQ, "\n Starting mie_write_opc, output will be written to %s\n\n",chFile );

	mie_auxiliary(&sd,"init");

	/* number of protons in plasma per average grain volume */
	gd.norm = sd.vol*gd.rho/(ATOMIC_MASS_UNIT*gd.mol_weight*gd.abun*gd.depl);
	volnorm = sd.vol;
	volfrac = 1.;

	acs_abs = (double **)MALLOC(sizeof(double *)*(unsigned)sd.nPart);
	acs_sct = (double **)MALLOC(sizeof(double *)*(unsigned)sd.nPart);
	a1g = (double **)MALLOC(sizeof(double *)*(unsigned)sd.nPart);
	inv_att_len = (double *)MALLOC(sizeof(double)*(unsigned)rfield.nupper);

	if( acs_abs == NULL || acs_sct == NULL || a1g == NULL || inv_att_len == NULL ) 
	{
		fprintf( ioQQQ, " Could not MALLOC arrays\n" );
		puts( "[Stop in mie_write_opc]" );
		cdEXIT(EXIT_FAILURE);
	}
	for( p=0; p < sd.nPart; p++ ) 
	{
		acs_abs[p] = (double *)MALLOC(sizeof(double)*(unsigned)rfield.nupper);
		acs_sct[p] = (double *)MALLOC(sizeof(double)*(unsigned)rfield.nupper);
		a1g[p] = (double *)MALLOC(sizeof(double)*(unsigned)rfield.nupper);
		if( acs_abs[p] == NULL || acs_sct[p] == NULL || a1g[p] == NULL ) 
		{
			fprintf( ioQQQ, " Could not MALLOC arrays\n" );
			puts( "[Stop in mie_write_opc]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	if( (fdes = fopen( chFile, "w" )) == NULL ) 
	{
		fprintf( ioQQQ, " Could not open %s for writing\n", chFile );
		puts( "[Stop in mie_write_opc]" );
		cdEXIT(EXIT_FAILURE);
	}
	lgErr = FALSE;

	(void)time(&timer);
	lgErr = lgErr || ( fprintf(fdes,"# this file was created by Cloudy %s (%s) on %s",
				   version.chVersion,version.chDate,ctime(&timer)) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n#\n") < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%12ld # magic number opacity file\n",MAGIC_OPC) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%12ld # magic number rfi/mix file\n",gd.magic) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%12ld # magic number szd file\n",sd.magic) < 0 );

	/* generate grain label for Cloudy output
	 * adjust LABELSIZE in mie.h when the format defined below is changed ! */
	strncpy(hlp1,chFile,(size_t)(LABELSUB1+1));
	hlp1[LABELSUB1+1] = '\0';
	str = strstr(hlp1,"-");

	if( str != NULL )
		*str = '\0';

	strncpy(hlp2,chFile2,(size_t)(LABELSUB2+1));
	hlp2[LABELSUB2+1] = '\0';
	str = strstr(hlp2,"-");

	if( str != NULL )
		*str = '\0';

	strcpy(chGrainLabel," ");
	if( sd.nPart > 1 ) 
	{
		hlp1[LABELSUB1] = '\0';
		hlp2[LABELSUB2] = '\0';
		strcat(strcat(strcat(strcat(chGrainLabel,hlp1),"-"),hlp2),"xx");
		lgErr = lgErr || ( fprintf(fdes,"%-12.12s # grain type label, xx will be replaced by bin no.\n",
					   chGrainLabel) < 0 );
	}
	else 
	{
		strcat(strcat(strcat(chGrainLabel,hlp1),"-"),hlp2);
		lgErr = lgErr || ( fprintf(fdes,"%-12.12s # grain type label\n", chGrainLabel) < 0 );
	}

	lgErr = lgErr || ( fprintf(fdes,"%.6e # specific weight (g/cm^3)\n",gd.rho) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%.6e # molecular weight of grain molecule (amu)\n",gd.mol_weight) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%.6e # average molecular weight per atom (amu)\n", gd.atom_weight) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%.6e # abundance of grain molecule at max depletion\n",gd.abun) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%.6e # depletion efficiency\n",gd.depl) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%.6e # average grain radius <a^3>/<a^2>, full size distr (cm)\n",
			   3.*sd.vol/sd.area) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%.6e # average grain surface area <4pi*a^2>, full size distr (cm^2)\n",
			   sd.area) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%.6e # average grain volume <4/3pi*a^3>, full size distr (cm^3)\n",
			   sd.vol) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain radius Int(a) per H, full size distr (cm/H)\n",
			   sd.radius/gd.norm) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain area Int(4pi*a^2) per H, full size distr (cm^2/H)\n",
			   sd.area/gd.norm) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain vol Int(4/3pi*a^3) per H, full size distr (cm^3/H)\n",
			   sd.vol/gd.norm) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%.6e # work function (Ryd)\n",gd.work) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%.6e # gap between valence and conduction band (Ryd)\n",gd.bandgap) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%.6e # efficiency of thermionic emission\n",gd.therm_eff) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%.6e # sublimation temperature (K)\n",gd.subl_temp) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%12d # material type, 1=carbonaceous, 2=silicate\n",gd.matType) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"#\n# abundances of constituent elements rel. to hydrogen\n#\n") < 0 );

	for( nelem=0; nelem < LIMELM; nelem++ ) 
	{
		lgErr = lgErr || ( fprintf(fdes,"%.6e # %s\n",gd.elmAbun[nelem],
					   elementnames.chElementSym[nelem]) < 0 );
	}

	if( sd.sdCase != SD_SINGLE_SIZE )
	{
		lgErr = lgErr || ( fprintf(fdes,"#\n# entire size distribution, amin=%.5f amax=%.5f micron\n",
					   sd.lim[ipBLo],sd.lim[ipBHi]) < 0 );
		lgErr = lgErr || ( fprintf(fdes,"#\n%.6e # ratio a_max/a_min in each size bin\n",
					   pow(sd.lim[ipBHi]/sd.lim[ipBLo],1./(double)sd.nPart) ) < 0 );
		lgErr = lgErr || ( fprintf(fdes,"#\n# size distribution function\n#\n") < 0 );
		lgErr = lgErr || ( fprintf(fdes,"%12ld # number of table entries\n#\n",NPTS_TABLE+1) < 0 );
		lgErr = lgErr || ( fprintf(fdes,"# ============================\n") < 0 );
		lgErr = lgErr || ( fprintf(fdes,"# size (micr) a^4*dN/da (cm^3/H)\n#\n") < 0 );
		for( i=0; i <= NPTS_TABLE; i++ )
		{
			double radius, a4dNda;
			radius = sd.lim[ipBLo]*exp((double)i/(double)NPTS_TABLE*log(sd.lim[ipBHi]/sd.lim[ipBLo]));
			radius = MAX2(MIN2(radius,sd.lim[ipBHi]),sd.lim[ipBLo]);
			a4dNda = POW4(radius)*size_distr(radius,&sd)/gd.norm*1.e-12/sd.unity;
			lgErr = lgErr || ( fprintf(fdes,"%.6e %.6e\n",radius,a4dNda) < 0 );
		}
	}
	else
	{
		lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
		lgErr = lgErr || ( fprintf(fdes,"%.6e # a_max/a_min = 1 for single sized grain\n", 1. ) < 0 );
		lgErr = lgErr || ( fprintf(fdes,"%12ld # no size distribution table\n",0L) < 0 );
	}

	lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%12ld # rfield.nupper\n",rfield.nupper) < 0 );
	lgErr = lgErr || ( fprintf(fdes,"%12ld # number of size distr. bins\n#\n",sd.nPart) < 0 );

	for( p=0; p < sd.nPart; p++ ) 
	{
		sd.cPart = p;

		mie_auxiliary(&sd,"step");

		if( sd.nPart > 1 ) 
		{
			/* >>chng 01 mar 20, creating mie_integrate introduced a change in the normalization
			 * of sd.radius, sd.area, and sd.vol; they now give average quantities for this bin.
			 * gd.norm converts average quanties to integrated quantities per H assuming the
			 * number of grains for the entire size distribution, hence multiplication by frac is
			 * needed to convert to the number of grains in this particular size bin, PvH */
			double frac = sd.unity_bin/sd.unity;
			volfrac = sd.vol*frac/volnorm;
			fprintf( ioQQQ, " Starting size bin %ld, amin=%.5f amax=%.5f micron\n",
				 p+1,sd.clim[ipBLo],sd.clim[ipBHi] );
			lgErr = lgErr || ( fprintf(fdes,"# size bin %ld, amin=%.5f amax=%.5f micron\n",
						   p+1,sd.clim[ipBLo],sd.clim[ipBHi]) < 0 );
			lgErr = lgErr || ( fprintf(fdes,"%.6e # average grain ",3.*sd.vol/sd.area) < 0 );
			lgErr = lgErr || ( fprintf(fdes,"radius <a^3>/<a^2>, this bin (cm)\n") < 0 );
			lgErr = lgErr || ( fprintf(fdes,"%.6e # average ",sd.area) < 0 );
			lgErr = lgErr || ( fprintf(fdes,"grain area <4pi*a^2>, this bin (cm^2)\n") < 0 );
			lgErr = lgErr || ( fprintf(fdes,"%.6e # average ",sd.vol) < 0 );
			lgErr = lgErr || ( fprintf(fdes,"grain volume <4/3pi*a^3>, this bin (cm^3)\n") < 0 );
			lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain ",sd.radius*frac/gd.norm) < 0 );
			lgErr = lgErr || ( fprintf(fdes,"radius Int(a) per H, this bin (cm/H)\n") < 0 );
			lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain area ",sd.area*frac/gd.norm) < 0 );
			lgErr = lgErr || ( fprintf(fdes,"Int(4pi*a^2) per H, this bin (cm^2/H)\n") < 0 );
			lgErr = lgErr || ( fprintf(fdes,"%.6e # total grain volume ",sd.vol*frac/gd.norm) < 0 );
			lgErr = lgErr || ( fprintf(fdes,"Int(4/3pi*a^3) per H, this bin (cm^3/H)\n#\n") < 0 );
		}

		lgErrorOccurred = FALSE;

		/* >> chng 02 sep 18, created infrastructure for new special files like PAH's, PvH */
		for( i=0; i < rfield.nupper; i++ ) 
		{
			wavlen = WAVNRYD/rfield.anu[i]*1.e4;

			switch( gd.rfiType )
			{
			case RFI_TABLE:
				ErrorIndex[i] = 0;
				acs_abs[p][i] = 0.;
				acs_sct[p][i] = 0.;
				a1g[p][i] = 0.;

				for( gd.cAxis=0; gd.cAxis < gd.nAxes; gd.cAxis++ ) 
				{
					mie_cs_size_distr(wavlen,&sd,&gd,mie_cs,&cs_abs,&cs_sct,&cosb,&Error);
					ErrorIndex[i] = MAX2(ErrorIndex[i],Error);
					acs_abs[p][i] += cs_abs*gd.wt[gd.cAxis];
					acs_sct[p][i] += cs_sct*gd.wt[gd.cAxis];
					a1g[p][i] += cs_sct*(1.-cosb)*gd.wt[gd.cAxis];
				}

				if( ErrorIndex[i] > 0 ) 
				{
					ErrorIndex[i] = MIN2(ErrorIndex[i],2);
					lgErrorOccurred = TRUE;
				}

				switch( ErrorIndex[i] ) 
				{
				/*lint -e616 */
				case 2:
					acs_abs[p][i] = 0.;
					acs_sct[p][i] = 0.;
					/*lint -fallthrough */
					/* controls is supposed to flow to the next case */
				case 1:
					a1g[p][i] = 0.;
					break;
				/*lint +e616 */
				case 0:
					a1g[p][i] /= acs_sct[p][i];
					break;
				default:
					fprintf( ioQQQ, " Insane value for ErrorIndex: %d\n", ErrorIndex[i] );
					ShowMe();
					puts( "[Stop in mie_write_opc]" );
					cdEXIT(EXIT_FAILURE);
				}

				/* sanity checks */
				if( ErrorIndex[i] < 2 )
					ASSERT( acs_abs[p][i] > 0. && acs_sct[p][i] > 0. );
				if( ErrorIndex[i] < 1 )
					ASSERT( a1g[p][i] > 0. );
				break;
			case OPC_TABLE:
				/* >>chng 02 oct 22, added OPC_TABLE case, PvH */
				gd.cAxis = 0;
				mie_cs_size_distr(wavlen,&sd,&gd,tbl_fun,&cs_abs,&cs_sct,&cosb,&Error);
				ErrorIndex[i] = MIN2(Error,2);
				lgErrorOccurred = lgErrorOccurred || ( Error > 0 );
				acs_abs[p][i] = cs_abs*gd.norm;
				acs_sct[p][i] = cs_sct*gd.norm;
				a1g[p][i] = 1.-cosb;
				break;
			case OPC_GREY:
				ErrorIndex[i] = 0;
				acs_abs[p][i] = 1.3121e-23*volfrac*gd.norm;
				acs_sct[p][i] = 2.6242e-23*volfrac*gd.norm;
				a1g[p][i] = 1.;
				break;
			case OPC_PAH1:
				/* >>chng 02 oct 17, added OPC_PAH1 case, PvH */
				gd.cAxis = 0;
				mie_cs_size_distr(wavlen,&sd,&gd,pah1_fun,&cs_abs,&cs_sct,&cosb,&Error);
				ErrorIndex[i] = MIN2(Error,2);
				lgErrorOccurred = lgErrorOccurred || ( Error > 0 );
				acs_abs[p][i] = cs_abs;
				acs_sct[p][i] = cs_sct;
				a1g[p][i] = 1.-cosb;
				break;
			default:
				fprintf( ioQQQ, " Insanity in mie_write_opc\n" );
				ShowMe();
				puts( "[Stop in mie_write_opc]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		/* extrapolate/interpolate for missing data */
		if( lgErrorOccurred ) 
		{
			strcpy(string,"absorption cs");
			mie_repair(string,rfield.nupper,2,0,rfield.anu,acs_abs[p],ErrorIndex,FALSE,&lgWarning);
			strcpy(string,"scattering cs");
			mie_repair(string,rfield.nupper,2,1,rfield.anu,acs_sct[p],ErrorIndex,FALSE,&lgWarning);
			strcpy(string,"asymmetry parameter");
			mie_repair(string,rfield.nupper,1,1,rfield.anu,a1g[p],ErrorIndex,TRUE,&lgWarning);
		}

		for( i=0; i < rfield.nupper; i++ ) 
		{
			acs_abs[p][i] /= gd.norm;
			/* >>chng 02 dec 30, do not multiply with (1-g) and write this factor out
			 * separately; this is useful for calculating extinction properties of grains, PvH */
			acs_sct[p][i] /= gd.norm;
		}
#if 0
		{
			FILE* ftmp;
			char name[20] = "asymxx";
			sprintf(&name[4],"%2.2li",p+1);
			ftmp = fopen(name,"w");
			for( i=0; i < rfield.nupper; i++ ) 
			{
				fprintf( ftmp, "%.6e %.6e\n", rfield.anu[i],a1g[p][i]);
			}
			fclose(ftmp);
		}
#endif

		mie_auxiliary(&sd,"cleanup");
	}

	lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
	lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n") < 0 );
	lgErr = lgErr || ( fprintf(fdes,"# anu (Ryd) abs_cs_01 (cm^2/H) abs_cs_02.....\n#\n") < 0 );

	for( i=0; i < rfield.nupper; i++ ) 
	{
		lgErr = lgErr || ( fprintf(fdes,"%.6e ",rfield.anu[i]) < 0 );
		for( p=0; p < sd.nPart; p++ ) 
		{
			lgErr = lgErr || ( fprintf(fdes,"%.6e ",acs_abs[p][i]) < 0 );
		}
		lgErr = lgErr || ( fprintf(fdes,"\n") < 0 );
	}

	lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
	lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n") < 0 );
	lgErr = lgErr || ( fprintf(fdes,"# anu (Ryd) sct_cs_01 (cm^2/H) sct_cs_02.....\n#\n") < 0 );

	for( i=0; i < rfield.nupper; i++ ) 
	{
		lgErr = lgErr || ( fprintf(fdes,"%.6e ",rfield.anu[i]) < 0 );
		for( p=0; p < sd.nPart; p++ ) 
		{
			lgErr = lgErr || ( fprintf(fdes,"%.6e ",acs_sct[p][i]) < 0 );
		}
		lgErr = lgErr || ( fprintf(fdes,"\n") < 0 );
	}

	lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
	lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n") < 0 );
	lgErr = lgErr || ( fprintf(fdes,"# anu (Ryd) (1-g)_bin_01 (1-g)_bin_02.....\n#\n") < 0 );

	for( i=0; i < rfield.nupper; i++ ) 
	{
		lgErr = lgErr || ( fprintf(fdes,"%.6e ",rfield.anu[i]) < 0 );
		for( p=0; p < sd.nPart; p++ ) 
		{
			lgErr = lgErr || ( fprintf(fdes,"%.6e ",a1g[p][i]) < 0 );
		}
		lgErr = lgErr || ( fprintf(fdes,"\n") < 0 );
	}

	fprintf( ioQQQ, " Starting calculation of inverse attenuation length\n" );
	strcpy(string,"inverse attenuation length");
	if( gd.rfiType != RFI_TABLE ) 
	{
		/* >>chng 02 sep 18, added graphite case for special files like PAH's, PvH */
		ial_type icase = gv.which_ial[gd.matType];
		switch( icase )
		{
		case IAL_CAR:
			mie_read_rfi("graphite.rfi",&gd2);
			mie_calc_ial(&gd2,rfield.nupper,inv_att_len,string,&lgWarning);
			break;
		case IAL_SIL:
			mie_read_rfi("silicate.rfi",&gd2);
			mie_calc_ial(&gd2,rfield.nupper,inv_att_len,string,&lgWarning);
			break;
		default:
			fprintf( ioQQQ, " mie_write_opc detected unknown ial type: %d\n" , icase );
			puts( "[Stop in mie_write_opc]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
	else 
	{
		mie_calc_ial(&gd,rfield.nupper,inv_att_len,string,&lgWarning);
	}

	lgErr = lgErr || ( fprintf(fdes,"#\n") < 0 );
	lgErr = lgErr || ( fprintf(fdes,"# ===========================================\n") < 0 );
	lgErr = lgErr || ( fprintf(fdes,"# anu (Ryd) inverse attenuation length (cm^-1)\n#\n") < 0 );

	for( i=0; i < rfield.nupper; i++ ) 
	{
		lgErr = lgErr || ( fprintf(fdes,"%.6e %.6e\n",rfield.anu[i],inv_att_len[i]) < 0 );
	}

	fclose(fdes);

	if( lgErr ) 
	{
		fprintf( ioQQQ, "\n Error writing file: %s\n", chFile );
		if( remove(chFile) == 0 )
		{
			fprintf( ioQQQ, " The file has been removed\n" );
			puts( "[Stop in mie_write_opc]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
	else 
	{
		fprintf( ioQQQ, "\n Opacity file %s written succesfully\n\n", chFile );
		if( lgWarning )
		{
			fprintf( ioQQQ, "\n !!! Warnings were detected !!!\n\n" );
		}
	}

	for( p=0; p < sd.nPart; p++ ) 
	{
		free(a1g[p]);
		free(acs_sct[p]);
		free(acs_abs[p]);
	}

	free(inv_att_len);
	free(a1g);
	free(acs_sct);
	free(acs_abs);

	/* these arrays were allocated in mie_read_szd */
	if( sd.ln_a != NULL )
		free(sd.ln_a);
	if( sd.ln_a4dNda != NULL )
		free(sd.ln_a4dNda);

	/* these arrays were allocated in mie_read_rfi */
	for( j=0; j < gd2.nOpcCols; j++ ) 
	{
		if( gd2.opcData[j] != NULL )
			free(gd2.opcData[j]);
	}
	if( gd2.opcAnu != NULL )
		free(gd2.opcAnu);
	for( j=0; j < gd2.nAxes; j++ ) 
	{
		if( gd2.nr1[j] != NULL )
			free(gd2.nr1[j]);
		if( gd2.n[j] != NULL )
			free(gd2.n[j]);
		if( gd2.wavlen[j] != NULL )
			free(gd2.wavlen[j]);
	}

	for( j=0; j < gd.nOpcCols; j++ ) 
	{
		if( gd.opcData[j] != NULL )
			free(gd.opcData[j]);
	}
	if( gd.opcAnu != NULL )
		free(gd.opcAnu);
	for( j=0; j < gd.nAxes; j++ ) 
	{
		if( gd.nr1[j] != NULL )
			free(gd.nr1[j]);
		if( gd.n[j] != NULL )
			free(gd.n[j]);
		if( gd.wavlen[j] != NULL )
			free(gd.wavlen[j]);
	}

	free( ErrorIndex );

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


static void mie_auxiliary(/*@partial@*/ sd_data *sd,
			  /*@in@*/ const char *auxCase)
{
	double amin,
	  amax,
	  delta,
	  oldvol,
	  step;

	/* desired relative accuracy of integration over size distribution */
	const double TOLER = 1.e-3;

#	ifdef DEBUG_FUN
	fputs( "<+>mie_auxiliary()\n", debug_fp );
#	endif
	if( strcmp(auxCase,"init") == 0 )
	{
		sd->xx = NULL;
		sd->aa = NULL;
		sd->rr = NULL;
		sd->ww = NULL;

		/* this is the initial estimate for the multiplier needed to get the
		 * number of abscissas in the gaussian quadrature, the correct value
		 * will be iterated below */
		sd->nmul = 1;

		/* calculate average grain surface area and volume over size distribution */
		switch( sd->sdCase ) 
		{
		case SD_SINGLE_SIZE:
			sd->radius = sd->a[ipSize]*1.e-4;
			sd->area = 4.*PI*POW2(sd->a[ipSize])*1.e-8;
			sd->vol = 4./3.*PI*POW3(sd->a[ipSize])*1.e-12;
			break;
		case SD_POWERLAW:
		case SD_EXP_CUTOFF1:
		case SD_EXP_CUTOFF2:
		case SD_EXP_CUTOFF3:
		case SD_LOG_NORMAL:
		case SD_LIN_NORMAL:
		case SD_TABLE:
			/* set up Gaussian quadrature for entire size range,
			 * first estimate no. of abscissas needed */
			amin = sd->lgLogScale ? log(sd->lim[ipBLo]) : sd->lim[ipBLo];
			amax = sd->lgLogScale ? log(sd->lim[ipBHi]) : sd->lim[ipBHi];

			sd->clim[ipBLo] = sd->lim[ipBLo];
			sd->clim[ipBHi] = sd->lim[ipBHi];

			/* iterate nmul until the integrated volume has converged sufficiently */
			oldvol= 0.;
			do 
			{
				sd->nmul *= 2;
				mie_integrate(sd,amin,amax,&sd->unity,TRUE);
				delta = fabs(sd->vol-oldvol)/sd->vol;
				oldvol = sd->vol;
			} while( sd->nmul <= 1024 && delta > TOLER );

			if( delta > TOLER )
			{
				fprintf( ioQQQ, " could not converge integration of size distribution\n" );
				puts( "[Stop in mie_auxiliary]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* we can safely reduce nmul by a factor of 2 and
			 * still reach a relative accuracy of TOLER */
			sd->nmul /= 2;
			mie_integrate(sd,amin,amax,&sd->unity,TRUE);
			break;
		case SD_ILLEGAL:
		default:
			fprintf( ioQQQ, " insane case for grain size distribution: %d\n" , sd->sdCase );
			ShowMe();
			puts( "[Stop in mie_auxiliary]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
	else if( strcmp(auxCase,"step") == 0 ) 
	{
		/* calculate average grain surface area and volume over size bin */
		switch( sd->sdCase ) 
		{
		case SD_SINGLE_SIZE:
			break;
		case SD_POWERLAW:
		case SD_EXP_CUTOFF1:
		case SD_EXP_CUTOFF2:
		case SD_EXP_CUTOFF3:
		case SD_LOG_NORMAL:
		case SD_LIN_NORMAL:
		case SD_TABLE:
			amin = sd->lgLogScale ? log(sd->lim[ipBLo]) : sd->lim[ipBLo];
			amax = sd->lgLogScale ? log(sd->lim[ipBHi]) : sd->lim[ipBHi];
			step = (amax - amin)/(double)sd->nPart;
			amin = amin + (double)sd->cPart*step;
			amax = MIN2(amax,amin + step);

			sd->clim[ipBLo] = sd->lgLogScale ? exp(amin) : amin;
			sd->clim[ipBHi] = sd->lgLogScale ? exp(amax) : amax;

			mie_integrate(sd,amin,amax,&sd->unity_bin,FALSE);

			break;
		case SD_ILLEGAL:	
		default:
			fprintf( ioQQQ, " insane case for grain size distribution: %d\n" , sd->sdCase );
			ShowMe();
			puts( "[Stop in mie_auxiliary]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
	else if( strcmp(auxCase,"cleanup") == 0 ) 
	{
		if( sd->xx != NULL )
			free(sd->xx);
		sd->xx = NULL;
		if( sd->aa != NULL )
			free(sd->aa);
		sd->aa = NULL;
		if( sd->rr != NULL )
			free(sd->rr);
		sd->rr = NULL;
		if( sd->ww != NULL )
			free(sd->ww);
		sd->ww = NULL;
	}
	else 
	{
		fprintf( ioQQQ, " mie_auxiliary called with insane argument: %s\n", auxCase );
		ShowMe();
		puts( "[Stop in mie_auxiliary]" );
		cdEXIT(EXIT_FAILURE);
	}

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

static void mie_integrate(/*@partial@*/ sd_data *sd,
			  double amin,
			  double amax,
			  /*@out@*/ double *normalization,
			  int lgFreeMem)
{
	long int j;
	double unity;

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

	/* set up Gaussian quadrature for size range,
	 * first estimate no. of abscissas needed */
	sd->nn = sd->nmul*((long)(2.*log(sd->clim[ipBHi]/sd->clim[ipBLo])) + 1);
	sd->nn = MIN2(MAX2(sd->nn,2*sd->nmul),4096);
	sd->xx = (double *)MALLOC(sizeof(double)*(unsigned)sd->nn);
	sd->aa = (double *)MALLOC(sizeof(double)*(unsigned)sd->nn);
	sd->rr = (double *)MALLOC(sizeof(double)*(unsigned)sd->nn);
	sd->ww = (double *)MALLOC(sizeof(double)*(unsigned)sd->nn);
	if( sd->xx == NULL || sd->aa == NULL || sd->rr == NULL || sd->ww == NULL ) 
	{
		fprintf( ioQQQ, " Could not MALLOC arrays\n" );
		puts( "[Stop in mie_integrate]" );
		cdEXIT(EXIT_FAILURE);
	}
	gauss_legendre(sd->nn,sd->xx,sd->aa);
	gauss_init(sd->nn,amin,amax,sd->xx,sd->aa,sd->rr,sd->ww);

	/* now integrate surface area and volume */
	unity = 0.;
	sd->radius = 0.;
	sd->area = 0.;
	sd->vol = 0.;

	for( j=0; j < sd->nn; j++ ) 
	{
		double weight;

		/* use extra factor of size in weights when we use logarithmic mesh */
		if( sd->lgLogScale ) 
		{
			sd->rr[j] = exp(sd->rr[j]);
			sd->ww[j] *= sd->rr[j];
		}
		weight = sd->ww[j]*size_distr(sd->rr[j],sd);
		unity += weight;
		sd->radius += weight*sd->rr[j];
		sd->area += weight*POW2(sd->rr[j]);
		sd->vol += weight*POW3(sd->rr[j]);
	}
	*normalization = unity;
	sd->radius *= 1.e-4/unity;
	sd->area *= 4.*PI*1.e-8/unity;
	sd->vol *= 4./3.*PI*1.e-12/unity;

	if( lgFreeMem )
	{
		if( sd->xx != NULL )
			free(sd->xx);
		sd->xx = NULL;
		if( sd->aa != NULL )
			free(sd->aa);
		sd->aa = NULL;
		if( sd->rr != NULL )
			free(sd->rr);
		sd->rr = NULL;
		if( sd->ww != NULL )
			free(sd->ww);
		sd->ww = NULL;
	}

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

/* read in the *.opc file with opacities and other relevant information */
void mie_read_opc(/*@in@*/const char *chFile,
		  GrainPar gp)
{
	int res,
	  lgDefaultQHeat;
	long int dl,
	  help,
	  i,
	  nelem,
	  j,
	  magic,
	  nbin,
	  nd,
	  nd2,
	  nup;
	float RefAbund[LIMELM],
	  VolTotal;
	double anu;
	double RadiusRatio;
	char chLine[FILENAME_PATH_LENGTH_2],
	  *str,
	  string[FILENAME_PATH_LENGTH_2];
	FILE *io2;

	/* if a_max/a_min in a single size bin is less than
	 * RATIO_MAX quantum heating will be turned on by default */
	const double RATIO_MAX = pow(100.,1./3.);

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

	/* true is special path set, false if file in same dir */
	if( lgDataPathSet )
	{
		/* path is set, generate full path name with file */
		strcpy( string , chDataPath );
		strcat( string,  chFile );
	}
	else
	{
		/* path not set, look here */
		strcpy( string , chFile );
	}

	if( (io2 = fopen(string,"r")) == NULL )
	{
		/* check local dir since not on path above */
		strcpy( string , chFile );
		if( (io2 = fopen(string,"r")) == NULL )
		{
			fprintf( ioQQQ, " Could not open %s for reading\n",string);
			puts( "[Stop in mie_read_opc]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* include the name of the file we are reading in the Cloudy output */
	strcpy( chLine, "                                        " );
	if( strlen(chFile) <= 40 )
	{
		strncpy( &chLine[0], chFile, strlen(chFile) );
	}
	else
	{
		strncpy( &chLine[0], chFile, 37 );
		strncpy( &chLine[37], "...", 3 );
	}
	if( called.lgTalk )
		fprintf( ioQQQ, "                       * >>>> mie_read_opc reading file -- %40s <<<< *\n", chLine );

	/* >>chng 02 jan 30, check if file has already been read before, PvH */
	for( i=0; i < gv.ReadPtr; ++i )
	{
		if( strcmp(chFile,gv.ReadRecord[i]) == 0 )
		{
			fprintf( ioQQQ, " File %s has already been read before, was this intended ?\n", chFile );
			break;
		}
	}
	/* remember the name of the file we are reading now */
	if( gv.ReadPtr < MAX_READ_RECORDS )
	{
		strcpy(gv.ReadRecord[gv.ReadPtr],chFile);
		++gv.ReadPtr;
	}

	/* allocate memory for first bin */
	nd = NewGrainBin();

	dl = 0; /* line counter for input file */

	/* first read magic numbers */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_long(chFile,chLine,&magic,TRUE,dl);
	if( magic != MAGIC_OPC ) 
	{
		fprintf( ioQQQ, " Opacity file %s has obsolete magic number\n",chFile );
		fprintf( ioQQQ, " I found magic number %ld, but expected %ld on line #%ld\n",
			 magic,MAGIC_OPC,dl );
		fprintf( ioQQQ, " Please recompile this file\n" );
		puts( "[Stop in mie_read_opc]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* the following two magic numbers are for information only */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_long(chFile,chLine,&magic,TRUE,dl);

	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_long(chFile,chLine,&magic,TRUE,dl);

	/* the grain scale factor is set equal to the abundances scale factor
	 * that might have appeared on the grains command.  Later, in grains.c, 
	 * it will be further multiplied by gv.GrainMetal, the scale factor that
	 * appears on the metals & grains command.  That command may, or may not,
	 * have been parsed yet, so can't do it at this stage. */
	gv.bin[nd]->dstfactor = (float)gp.dep;

	/* grain type label */
	mie_next_data(chFile,io2,chLine,&dl);
	strncpy(gv.bin[nd]->chDstLab,chLine,(size_t)LABELSIZE);
	gv.bin[nd]->chDstLab[LABELSIZE] = '\0';

	/* specific weight (g/cm^3) */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_float(chFile,chLine,&gv.bin[nd]->dustp[0],TRUE,dl);

	/* molecular weight of grain molecule (amu) */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_float(chFile,chLine,&gv.bin[nd]->dustp[1],TRUE,dl);

	/* average molecular weight per atom (amu) */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_float(chFile,chLine,&gv.bin[nd]->atomWeight,TRUE,dl);

	/* abundance of grain molecule for max depletion */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_float(chFile,chLine,&gv.bin[nd]->dustp[2],TRUE,dl);

	/* depletion efficiency */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_float(chFile,chLine,&gv.bin[nd]->dustp[3],TRUE,dl);

	/* average grain radius <a^3>/<a^2> for entire size distr (cm) */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_float(chFile,chLine,&gv.bin[nd]->AvRadius,TRUE,dl);

	/* average grain area <4pi*a^2> for entire size distr (cm^2) */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_float(chFile,chLine,&gv.bin[nd]->AvArea,TRUE,dl);

	/* average grain volume <4/3pi*a^3> for entire size distr (cm^3) */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_float(chFile,chLine,&gv.bin[nd]->AvVol,TRUE,dl);

	/* total grain radius Int(a) per H for entire size distr (cm/H) */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_float(chFile,chLine,&gv.bin[nd]->IntRadius,TRUE,dl);

	/* total grain area Int(4pi*a^2) per H for entire size distr (cm^2/H) */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_float(chFile,chLine,&gv.bin[nd]->IntArea,TRUE,dl);

	/* total grain vol Int(4/3pi*a^3) per H for entire size distr (cm^3/H) */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_float(chFile,chLine,&gv.bin[nd]->IntVol,TRUE,dl);

	/* work function, in Rydberg */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_float(chFile,chLine,&gv.bin[nd]->DustWorkFcn,TRUE,dl);

	/* bandgap, in Rydberg */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_float(chFile,chLine,&gv.bin[nd]->BandGap,FALSE,dl);

	/* efficiency of thermionic emissions, between 0 and 1 */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_float(chFile,chLine,&gv.bin[nd]->ThermEff,TRUE,dl);

	/* sublimation temperature in K */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_float(chFile,chLine,&gv.bin[nd]->Tsublimat,TRUE,dl);

	/* material type, determines enthalpy function, etc... */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_long(chFile,chLine,&help,TRUE,dl);
	gv.bin[nd]->matType = (mat_type)help;
	
	for( nelem=0; nelem < LIMELM; nelem++ ) 
	{
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_float(chFile,chLine,&RefAbund[nelem],FALSE,dl);

		/* this coefficient is defined at the end of appendix A.10 of BFM */
		gv.bin[nd]->AccomCoef[nelem] = 2.*gv.bin[nd]->atomWeight*dense.AtomicWeight[nelem]/
			POW2(gv.bin[nd]->atomWeight+dense.AtomicWeight[nelem]);
	}

	/* ratio a_max/a_min for grains in a single size bin */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_double(chFile,chLine,&RadiusRatio,TRUE,dl);

	gv.bin[nd]->lgDustVary = gp.lgAbunVsDepth;
	if( gp.lgAbunVsDepth )
		gv.lgAnyDustVary = TRUE;
	lgDefaultQHeat = ( RadiusRatio < RATIO_MAX && !gp.lgGreyGrain );
	gv.bin[nd]->lgQHeat = ( gp.lgForbidQHeating ) ? FALSE : ( gp.lgRequestQHeating || lgDefaultQHeat );
	gv.bin[nd]->cnv_H_pGR = gv.bin[nd]->AvVol/gv.bin[nd]->IntVol;
	gv.bin[nd]->cnv_GR_pH = 1./gv.bin[nd]->cnv_H_pGR;

	/* this is capacity per grain, in Farad per grain */
	gv.bin[nd]->Capacity = PI4*ELECTRIC_CONST*gv.bin[nd]->IntRadius/100.*gv.bin[nd]->cnv_H_pGR;

	/* skip the table of the size distribution function (if present) */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_long(chFile,chLine,&nup,FALSE,dl);
	for( i=0; i < nup; i++ )
		mie_next_data(chFile,io2,chLine,&dl);

	/* nup is number of frequency bins stored in file, this should match rfield.nupper */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_long(chFile,chLine,&nup,TRUE,dl);

	gv.bin[nd]->NFPCheck = nup;

	/* no. of size distribution bins */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_long(chFile,chLine,&nbin,TRUE,dl);

	if( nbin == 1 )
	{
		/* >>chng 01 sep 12, allocate/free [rfield.nupper] arrays dynamically */
		ASSERT( gv.bin[nd]->dstab1 == NULL ); /* prevent memory leaks */
		if( ( gv.bin[nd]->dstab1 = (double*)MALLOC(sizeof(double)*(unsigned)nup) ) == NULL ) 
			BadMalloc();
		ASSERT( gv.bin[nd]->pure_sc1 == NULL ); /* prevent memory leaks */
		if( ( gv.bin[nd]->pure_sc1 = (double*)MALLOC(sizeof(double)*(unsigned)nup) ) == NULL ) 
			BadMalloc();
		ASSERT( gv.bin[nd]->asym == NULL ); /* prevent memory leaks */
		if( ( gv.bin[nd]->asym = (double*)MALLOC(sizeof(double)*(unsigned)nup) ) == NULL ) 
			BadMalloc();
		ASSERT( gv.bin[nd]->inv_att_len == NULL ); /* prevent memory leaks */
		if( ( gv.bin[nd]->inv_att_len = (float*)MALLOC(sizeof(float)*(unsigned)nup) ) == NULL ) 
			BadMalloc();

		gv.bin[nd]->dustp[4] = 1.;
		for( nelem=0; nelem < LIMELM; nelem++ )
		{
			gv.bin[nd]->elmAbund[nelem] = RefAbund[nelem];
		}
	}
	else if( nbin > 1 )
	{
		/* remember this number since it will be overwritten below */
		VolTotal = gv.bin[nd]->IntVol;

		for( i=0; i < nbin; i++ ) 
		{
			/* allocate memory for remaining bins */
			nd2 = ( i == 0 ) ? nd : NewGrainBin();

			/* >>chng 01 sep 12, allocate/free [rfield.nupper] arrays dynamically */
			ASSERT( gv.bin[nd2]->dstab1 == NULL ); /* prevent memory leaks */
			if( ( gv.bin[nd2]->dstab1 = (double*)MALLOC(sizeof(double)*(unsigned)nup) ) == NULL ) 
				BadMalloc();
			ASSERT( gv.bin[nd2]->pure_sc1 == NULL ); /* prevent memory leaks */
			if( ( gv.bin[nd2]->pure_sc1 = (double*)MALLOC(sizeof(double)*(unsigned)nup) ) == NULL ) 
				BadMalloc();
			ASSERT( gv.bin[nd2]->asym == NULL ); /* prevent memory leaks */
			if( ( gv.bin[nd2]->asym = (double*)MALLOC(sizeof(double)*(unsigned)nup) ) == NULL ) 
				BadMalloc();
			ASSERT( gv.bin[nd2]->inv_att_len == NULL ); /* prevent memory leaks */
			if( ( gv.bin[nd2]->inv_att_len = (float*)MALLOC(sizeof(float)*(unsigned)nup) ) == NULL ) 
				BadMalloc();

			/* average grain radius <a^3>/<a^2> for this bin (cm) */
			mie_next_data(chFile,io2,chLine,&dl);
			mie_read_float(chFile,chLine,&gv.bin[nd2]->AvRadius,TRUE,dl);

			/* average grain area in this bin (cm^2) */
			mie_next_data(chFile,io2,chLine,&dl);
			mie_read_float(chFile,chLine,&gv.bin[nd2]->AvArea,TRUE,dl);

			/* average grain volume in this bin (cm^3) */
			mie_next_data(chFile,io2,chLine,&dl);
			mie_read_float(chFile,chLine,&gv.bin[nd2]->AvVol,TRUE,dl);

			/* total grain radius Int(a) per H for this bin (cm/H) */
			mie_next_data(chFile,io2,chLine,&dl);
			mie_read_float(chFile,chLine,&gv.bin[nd2]->IntRadius,TRUE,dl);

			/* total grain area Int(4pi*a^2) per H for this bin (cm^2/H) */
			mie_next_data(chFile,io2,chLine,&dl);
			mie_read_float(chFile,chLine,&gv.bin[nd2]->IntArea,TRUE,dl);

			/* total grain vol Int(4/3pi*a^3) per H for this bin (cm^3/H) */
			mie_next_data(chFile,io2,chLine,&dl);
			mie_read_float(chFile,chLine,&gv.bin[nd2]->IntVol,TRUE,dl);

			gv.bin[nd2]->cnv_H_pGR = gv.bin[nd2]->AvVol/gv.bin[nd2]->IntVol;
			gv.bin[nd2]->cnv_GR_pH = 1./gv.bin[nd2]->cnv_H_pGR;

			

			/* this is capacity per grain, in Farad per grain */
			gv.bin[nd2]->Capacity =
				PI4*ELECTRIC_CONST*gv.bin[nd2]->IntRadius/100.*gv.bin[nd2]->cnv_H_pGR;

			/* dustp[4] gives the fraction of the grain abundance that is
			 * contained in a particular bin. for unresolved distributions it is
			 * by definition 1, for resolved distributions it is smaller than 1. */
			gv.bin[nd2]->dustp[4] = gv.bin[nd2]->IntVol/VolTotal;
			for( nelem=0; nelem < LIMELM; nelem++ )
			{
				gv.bin[nd2]->elmAbund[nelem] = RefAbund[nelem]*gv.bin[nd2]->dustp[4];
			}

			if( i > 0 ) 
			{
				gv.bin[nd2]->dstfactor = gv.bin[nd]->dstfactor;
				strcpy(gv.bin[nd2]->chDstLab,gv.bin[nd]->chDstLab);
				gv.bin[nd2]->atomWeight = gv.bin[nd]->atomWeight;
				for( j=0; j < 4; j++ ) 
				{
					gv.bin[nd2]->dustp[j] = gv.bin[nd]->dustp[j];
				}
				gv.bin[nd2]->DustWorkFcn = gv.bin[nd]->DustWorkFcn;
				gv.bin[nd2]->BandGap = gv.bin[nd]->BandGap;
				gv.bin[nd2]->ThermEff = gv.bin[nd]->ThermEff;
				gv.bin[nd2]->Tsublimat = gv.bin[nd]->Tsublimat;
				gv.bin[nd2]->matType = gv.bin[nd]->matType;
				gv.bin[nd2]->lgDustVary = gv.bin[nd]->lgDustVary;
				gv.bin[nd2]->lgQHeat = gv.bin[nd]->lgQHeat;
				gv.bin[nd2]->NFPCheck = gv.bin[nd]->NFPCheck;
				for( nelem=0; nelem < LIMELM; nelem++ )
				{
					gv.bin[nd2]->AccomCoef[nelem] = gv.bin[nd]->AccomCoef[nelem];
				}
			}
		}
		for( i=0; i < nbin; i++ ) 
		{
			nd2 = nd + i;
			/* modify grain labels */
			str = strstr(gv.bin[nd2]->chDstLab,"xx");
			if( str != NULL )
				sprintf(str,"%02ld",i+1);
		}
	}

	/* skip the next 5 lines */
	for( i=0; i < 5; i++ )
		mie_next_line(chFile,io2,chLine,&dl);

	/* now read absorption opacities */
	for( i=0; i < nup; i++ ) 
	{
		/* read in energy scale and then opacities */
		if( (res = fscanf(io2,"%le",&anu)) != 1 ) 
		{
			fprintf( ioQQQ, " Read failed on %s\n",chFile );
			if( res == EOF )
				fprintf( ioQQQ, " EOF reached prematurely\n" );
			puts( "[Stop in mie_read_opc]" );
			cdEXIT(EXIT_FAILURE);
		}
		for( j=0; j < nbin; j++ ) 
		{
			nd2 = nd + j;
			if( (res = fscanf(io2,"%le",&gv.bin[nd2]->dstab1[i])) != 1 ) 
			{
				fprintf( ioQQQ, " Read failed on %s\n",chFile );
				if( res == EOF )
					fprintf( ioQQQ, " EOF reached prematurely\n" );
				puts( "[Stop in mie_read_opc]" );
				cdEXIT(EXIT_FAILURE);
			}
			ASSERT( gv.bin[nd2]->dstab1[i] > 0. );
		}
	}

	/* skip to end-of-line and then skip next 4 lines */
	for( i=0; i < 5; i++ )
		mie_next_line(chFile,io2,chLine,&dl);

	/* now read scattering opacities */
	for( i=0; i < nup; i++ ) 
	{
		if( (res = fscanf(io2,"%le",&anu)) != 1 ) 
		{
			fprintf( ioQQQ, " Read failed on %s\n",chFile );
			if( res == EOF )
				fprintf( ioQQQ, " EOF reached prematurely\n" );
			puts( "[Stop in mie_read_opc]" );
			cdEXIT(EXIT_FAILURE);
		}
		for( j=0; j < nbin; j++ ) 
		{
			nd2 = nd + j;
			if( (res = fscanf(io2,"%le",&gv.bin[nd2]->pure_sc1[i])) != 1 ) 
			{
				fprintf( ioQQQ, " Read failed on %s\n",chFile );
				if( res == EOF )
					fprintf( ioQQQ, " EOF reached prematurely\n" );
				puts( "[Stop in mie_read_opc]" );
				cdEXIT(EXIT_FAILURE);
			}
			ASSERT( gv.bin[nd2]->pure_sc1[i] > 0. );
		}
	}

	
		  	
	/* skip to end-of-line and then skip next 4 lines */
	for( i=0; i < 5; i++ )
		mie_next_line(chFile,io2,chLine,&dl);

	/* now read asymmetry factor */
	for( i=0; i < nup; i++ ) 
	{
		if( (res = fscanf(io2,"%le",&anu)) != 1 ) 
		{
			fprintf( ioQQQ, " Read failed on %s\n",chFile );
			if( res == EOF )
				fprintf( ioQQQ, " EOF reached prematurely\n" );
			puts( "[Stop in mie_read_opc]" );
			cdEXIT(EXIT_FAILURE);
		}
		for( j=0; j < nbin; j++ ) 
		{
			nd2 = nd + j;
			if( (res = fscanf(io2,"%le",&gv.bin[nd2]->asym[i])) != 1 ) 
			{
				fprintf( ioQQQ, " Read failed on %s\n",chFile );
				if( res == EOF )
					fprintf( ioQQQ, " EOF reached prematurely\n" );
				puts( "[Stop in mie_read_opc]" );
				cdEXIT(EXIT_FAILURE);
			}
			ASSERT( gv.bin[nd2]->asym[i] > 0. );
		}
	}

	/* skip to end-of-line and then skip next 4 lines */
	for( i=0; i < 5; i++ )
		mie_next_line(chFile,io2,chLine,&dl);

	/* now read inverse attenuation length */
	for( i=0; i < nup; i++ ) 
	{
		if( (res = fscanf(io2,"%le %e",&anu,&gv.bin[nd]->inv_att_len[i])) != 2 ) 
		{
			fprintf( ioQQQ, " Read failed on %s\n",chFile );
			if( res == EOF )
				fprintf( ioQQQ, " EOF reached prematurely\n" );
			puts( "[Stop in mie_read_opc]" );
			cdEXIT(EXIT_FAILURE);
		}
		ASSERT( gv.bin[nd]->inv_att_len[i] > 0. );
		/* save energy for later tests of energy grid */
		gv.bin[nd]->EnergyCheck = (float)anu;

		for( j=1; j < nbin; j++ ) 
		{
			nd2 = nd + j;
			gv.bin[nd2]->inv_att_len[i] = gv.bin[nd]->inv_att_len[i];
			gv.bin[nd2]->EnergyCheck = gv.bin[nd]->EnergyCheck;
		}
	}

	fclose(io2);

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


/* calculate average absorption, scattering cross section (i.e. pi a^2 Q) and
 * average asymmetry parameter g for an arbitrary grain size distribution */
static void mie_cs_size_distr(double wavlen, /* micron */
			      /*@in@*/ sd_data *sd,
			      /*@in@*/ grain_data *gd,
			      void(*cs_fun)(double,/*@in@*/sd_data*,/*@in@*/grain_data*,
					    /*@out@*/double*,/*@out@*/double*,
					    /*@out@*/double*,/*@out@*/int*),
			      /*@out@*/ double *cs_abs, /* cm^2, average */
			      /*@out@*/ double *cs_sct, /* cm^2, average */
			      /*@out@*/ double *cosb,
			      /*@out@*/ int *error)
{
	int lgADLused;
	long int i;
	double absval,
	  g,
	  sctval,
	  weight;

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

	/* sanity checks */
	ASSERT( wavlen > 0. );
	ASSERT( gd->cAxis >= 0 && gd->cAxis < gd->nAxes && gd->cAxis < NAX );

	switch( sd->sdCase ) 
	{
	case SD_SINGLE_SIZE:
		/* do single sized grain */
		ASSERT( sd->a[ipSize] > 0. );
		sd->cSize = sd->a[ipSize];
		(*cs_fun)(wavlen,sd,gd,cs_abs,cs_sct,cosb,error);
		break;
	case SD_POWERLAW:
		/* simple powerlaw distribution */
	case SD_EXP_CUTOFF1:
	case SD_EXP_CUTOFF2:
	case SD_EXP_CUTOFF3:
		/* powerlaw distribution with exponential cutoff */
	case SD_LOG_NORMAL:
		/* gaussian distribution in ln(a) */
	case SD_LIN_NORMAL:
		/* gaussian distribution in a */
	case SD_TABLE:
		/* user supplied table of a^4*dN/da */
		ASSERT( sd->lim[ipBLo] > 0. && sd->lim[ipBHi] > 0. && sd->lim[ipBHi] > sd->lim[ipBLo] );
		lgADLused = FALSE;
		*cs_abs = 0.;
		*cs_sct = 0.;
		*cosb = 0.;
		for( i=0; i < sd->nn; i++ ) 
		{
			sd->cSize = sd->rr[i];
			(*cs_fun)(wavlen,sd,gd,&absval,&sctval,&g,error);
			if( *error >= 2 ) 
			{
				/* mie_cs failed to converge -> integration is invalid */
				*cs_abs = -1.;
				*cs_sct = -1.;
				*cosb = -2.;

#				ifdef DEBUG_FUN
				fputs( " <->mie_cs_size_distr()\n", debug_fp );
#				endif
				return;
			}
			else if( *error == 1 ) 
			{
				/* anomalous diffraction limit used -> g is not valid */
				lgADLused = TRUE;
			}
			weight = sd->ww[i]*size_distr(sd->rr[i],sd);
			*cs_abs += weight*absval;
			*cs_sct += weight*sctval;
			*cosb += weight*sctval*g;
		}
		if( lgADLused ) 
		{
			*error = 1;
			*cosb = -2.;
		}
		else 
		{
			*error = 0;
			*cosb /= *cs_sct;
		}
		*cs_abs /= sd->unity;
		*cs_sct /= sd->unity;
		break;
	case SD_ILLEGAL:
	default:
		fprintf( ioQQQ, " insane case for grain size distribution: %d\n" , sd->sdCase );
		ShowMe();
		puts( "[Stop in mie_cs_size_distr]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* sanity checks */
	if( *error < 2 )
		ASSERT( *cs_abs > 0. && *cs_sct > 0. );
	if( *error < 1 )
		ASSERT( fabs(*cosb) <= 1.+10.*DBL_EPSILON );

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


/* calculate absorption, scattering cross section (i.e. pi a^2 Q) and
 * asymmetry parameter g (=cosb) for a single sized grain defined by gd */
static void mie_cs(double wavlen,  /* micron */
		   /*@in@*/ sd_data *sd,
		   /*@in@*/ grain_data *gd,
		   /*@out@*/ double *cs_abs, /* cm^2 */
		   /*@out@*/ double *cs_sct, /* cm^2 */
		   /*@out@*/ double *cosb,
		   /*@out@*/ int *error)
{
	int lgOutOfBounds;
	long int iflag,
	  ind;
	double area,
	  aqabs,
	  aqext,
	  aqphas,
	  beta,
	  ctbrqs,
	  delta,
	  frac,
	  nim,
	  nre,
	  nr1,
	  qback,
	  qext,
	  qphase,
	  qscatt,
	  x,
	  xistar;

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

	/* sanity checks, should already have been checked further upstream */
	ASSERT( wavlen > 0. );
	ASSERT( sd->cSize > 0. );
	ASSERT( gd->wavlen[gd->cAxis] != NULL && gd->ndata[gd->cAxis] > 1 );

	/* first interpolate optical constants */
	/* >>chng 02 oct 22, moved calculation of optical constants from mie_cs_size_distr to mie_cs,
	 * this increases overhead, but makes the code in mie_cs_size_distr more transparent, PvH */
	find_arr(wavlen,gd->wavlen[gd->cAxis],gd->ndata[gd->cAxis],&ind,&lgOutOfBounds);

	if( lgOutOfBounds ) 
	{
		*error = 3;
		*cs_abs = -1.;
		*cs_sct = -1.;
		*cosb = -2.;

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

	frac = (wavlen-gd->wavlen[gd->cAxis][ind])/(gd->wavlen[gd->cAxis][ind+1]-gd->wavlen[gd->cAxis][ind]);
	ASSERT( frac > 0.-10.*DBL_EPSILON && frac < 1.+10.*DBL_EPSILON );
	nre = (1.-frac)*gd->n[gd->cAxis][ind].re + frac*gd->n[gd->cAxis][ind+1].re;
	ASSERT( nre > 0. );
	nim = (1.-frac)*gd->n[gd->cAxis][ind].im + frac*gd->n[gd->cAxis][ind+1].im;
	ASSERT( nim > 0. );
	nr1 = (1.-frac)*gd->nr1[gd->cAxis][ind] + frac*gd->nr1[gd->cAxis][ind+1];
	ASSERT( fabs(nre-1.-nr1) < 10.*MAX2(nre,1.)*DBL_EPSILON );

	/* size in micron, area in cm^2 */
	area = PI*POW2(sd->cSize)*1.e-8;

	x = sd->cSize/wavlen*2.*PI;

	/* note that in the following, only nre, nim are used in sinpar
	 * and also nr1 in anomalous diffraction limit */

	sinpar(nre,nim,x,&qext,&qphase,&qscatt,&ctbrqs,&qback,&iflag);

	/* iflag=0 normal exit, 1 failure to converge, 2 not even tried
	 * for exit 1,2, see whether anomalous diffraction is available */

	if( iflag == 0 ) 
	{
		*error = 0;
		*cs_abs = area*(qext - qscatt);
		*cs_sct = area*qscatt;
		*cosb = ctbrqs/qscatt;
	}
	else 
	{
		/* anomalous diffraction -- x >> 1 and |m-1| << 1 but any phase shift */
		if( x >= 100. && sqrt(nr1*nr1+nim*nim) <= 0.001 ) 
		{
			delta = -nr1;
			beta = nim;

			anomal(x,&aqext,&aqabs,&aqphas,&xistar,delta,beta);

			/* cosb is invalid */
			*error = 1;
			*cs_abs = area*aqabs;
			*cs_sct = area*(aqext - aqabs);
			*cosb = -2.;
		}
		/* nothing works */
		else 
		{
			*error = 2;
			*cs_abs = -1.;
			*cs_sct = -1.;
			*cosb = -2.;
		}
	}
	if( *error < 2 ) 
	{
		if( *cs_abs <= 0. || *cs_sct <= 0. ) 
		{
			fprintf( ioQQQ, " illegal opacity found: wavl=%.4e micron," , wavlen );
			fprintf( ioQQQ, " abs_cs=%.2e, sct_cs=%.2e\n" , *cs_abs , *cs_sct );
			fprintf( ioQQQ, " please check refractive index file...\n" );
			puts( "[Stop in mie_cs]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
	if( *error < 1 ) 
	{
		if( fabs(*cosb) > 1.+10.*DBL_EPSILON ) 
		{
			fprintf( ioQQQ, " illegal asymmetry parameter found: wavl=%.4e micron," , wavlen );
			fprintf( ioQQQ, " cosb=%.2e\n" , *cosb );
			fprintf( ioQQQ, " please check refractive index file...\n" );
			puts( "[Stop in mie_cs]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
		

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


/* this routine calculates the absorption cross sections of PAH molecules, it is based on:
 * >>refer Desert, F.-X., Boulanger, F., Puget, J. L. 1990, A&A, 237, 215
 *
 * the original version of this routine was written by Kevin Volk (University Of Calgary) */

static const double pah1_strength[7]={1.4e-21,1.8e-21,1.2e-20,6.0e-21,4.0e-20,1.9e-20,1.9e-20};
static const double pah1_wlBand[7]={3.3, 6.18, 7.7, 8.6, 11.3, 12.0, 13.3};
static const double pah1_width[7]={0.024, 0.102, 0.24, 0.168, 0.086, 0.174, 0.174};

static void pah1_fun(double wavl,    /* in micron */
		     /*@in@*/ sd_data *sd,
		     /*@in@*/ grain_data *gd,
		     /*@out@*/ double *cs_abs,
		     /*@out@*/ double *cs_sct,
		     /*@out@*/ double *cosb,
		     /*@out@*/ int *error)
{
	long int j;
	double cval1,
	  pah1_fun_v,
	  term,
	  term1,
	  term2,
	  term3,
	  x;

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

	const double p1 = 4.0e-18;
	const double p2 = 1.1e-18;
	const double p3 = 3.3e-21;
	const double p4 = 6.0e-21;
	const double p5 = 2.4e-21;
	const double wl1a = 5.0;
	const double wl1b = 7.0;
	const double wl1c = 9.0;
	const double wl1d = 9.5;
	const double wl2a = 11.0;
	const double delwl2 = 0.3;
	/* this is the rise interval for the second plateau */
	const double wl2b = wl2a + delwl2;
	const double wl2c = 15.0;
	const double wl3a = 3.2;
	const double wl3b = 3.57;
	const double wl3m = (wl3a+wl3b)/2.;
	const double wl3sig = 0.1476;
	const double x1 = 4.0;
	const double x2 = 5.9;
	const double x2a = RYD_INF/1.e4;
	const double x3 = 0.1;

	/* grain volume */
	double vol = 4.*PI/3.*POW3(sd->cSize)*1.e-12;
	/* number of carbon atoms in PAH molecule */
	double xnc = floor(vol*gd->rho/(ATOMIC_MASS_UNIT*dense.AtomicWeight[ipCARBON]));
	/* number of hydrogen atoms in PAH molecule */
	/* >>chng 02 oct 18, use integral number of hydrogen atoms instead of fractional number */
	double xnh = floor(sqrt(6.*xnc));
	/* this is the hydrogen over carbon ratio in the PAH molecule */
	double xnhoc = xnh/xnc;
	/* ftoc3p3 is the feature to continuum ratio at 3.3 micron */
	double ftoc3p3 = 100.;

	double csVal1 = 0.;
	double csVal2 = 0.;

/* 	printf(" no. of atoms: C %.0f H %.0f\n",xnc,xnh); */

	/*TODO	2	impose size limits on PAH's, changeover to graphite ?? */
	/*TODO	2	detect if PAH's are present in ionized regions ?? */

	x = 1./wavl;

	if( x >= x2a )
	{
		/* >>chng 02 oct 18, use atomic cross sections for energies larger than 1 Ryd */
		double anu_ev = x/x2a*EVRYD;

		term1 = atmdat_phfit(ipHYDROGEN+1,1,1,anu_ev,FALSE);
		term2 = 0.;
		for( j=1; j <= 3; ++j )
			term2 += atmdat_phfit(ipCARBON+1,6,j,anu_ev,FALSE);

		csVal1 = (xnh*term1 + xnc*term2)*1.e-18;
	}
	
	if( x <= 2.*x2a )
	{
		cval1 = log(sqrt(xnc)*ftoc3p3/1.2328)/12.2;

		term = POW2(MIN2(x,x1))*(3.*x1 - 2.*MIN2(x,x1))/POW3(x1);

		term1 = (0.1*x + 0.41)*POW2(MAX2(x-x2,0.));

		/* The following is an exponential cut-off in the continuum for 
		 * wavelengths larger than 2500 Angstroms; it is exponential in 
		 * wavelength so it varies as 1/x here.  This replaces the 
		 * sharp cut-off at 8000 Angstroms in the original paper.
		 *
		 * This choice of continuum shape is also arbitrary.  The continuum
		 * is never observed at these wavelengths.  For the "standard" ratio 
		 * value at 3.3 microns the continuum level in the optical is not that 
		 * much smaller than it was in the original paper.  If one wants to 
		 * exactly reproduce the original optical opacity, one can change 
		 * the x1 value to a value of 1.125.  Then there will be a discontinuity
		 * in the cross-section at 8000 Angstroms.
		 *
		 * My judgement was that the flat cross-section in the optical used by 
		 * Desert, Boulander, and Puget (1990) is just a rough value that is not 
		 * based upon much in the way of direct observations, and so I could 
		 * change the cross-section at wavelengths above 2500 Angstroms.  It is
		 * likely that one should really build in the ERE somehow, judging from 
		 * the spectrum of the Red Rectangle, and there is no trace of this in
		 * the original paper.  The main concern in adding this exponential 
		 * drop-off in the continuum was to have a finite infrared continuum 
		 * between the features. */
		term2 = exp(cval1*(1. - (x1/MIN2(x,x1))));

		term3 = p3*exp(-POW2(x/x3))*MIN2(x,x3)/x3;

		csVal2 = xnc*((p1*term + p2*term1)*term2 + term3);
	}

	if( x2a <= x && x <= 2.*x2a )
	{
		/* create gradual change from Desert et al to atomic cross sections */
		double frac = POW2(2.-x/x2a);
		pah1_fun_v = exp((1.-frac)*log(csVal1) + frac*log(csVal2));
	}
	else
	{
		/* one of these will be zero */
		pah1_fun_v = csVal1 + csVal2;
	}

	/* now add in the three plateau features. the first two are based upon
	 * >>refer Schutte, Tielens, and Allamandola (1993) ApJ, 415, 397. */
	if( wl1a <= wavl && wl1d >= wavl )
	{
		if( wavl < wl1b )
			term = p4*(wavl - wl1a)/(wl1b - wl1a);
		else
			term = ( wavl > wl1c ) ? p4*(wl1d - wavl)/(wl1d - wl1c) : p4;
		pah1_fun_v += term*xnc;
	}
	if( wl2a <= wavl && wl2c >= wavl )
	{
		term = ( wavl < wl2b ) ? p5*((wavl - wl2a)/delwl2) : p5*sqrt(1.-POW2((wavl-wl2a)/(wl2c-wl2a)));
		pah1_fun_v += term*xnc;
	}
	if( wl3m-10.*wl3sig <= wavl && wavl <= wl3m+10.*wl3sig )
	{
		/* >>chng 02 nov 08, replace top hat distribution with gaussian, PvH */
/* 		term = 1.1*pah1_strength[0]/(wl3b - wl3a); */
		term = 1.1*pah1_strength[0]*exp(-0.5*POW2((wavl-wl3m)/wl3sig));
		pah1_fun_v += term*xnh;
	}

	/* add in the various discrete features in the infrared: 3.3, 6.2, 7.6, etc. */
	for( j=0; j < 7; j++ )
	{
		term1 = (wavl - pah1_wlBand[j])/pah1_width[j];
		term = 0.;
		if( j == 2 )
		{
			/* This assumes linear interpolation between the points, which are
			 * located at -1, -0.5, +1.5, and +3 times the width, or a fine spacing that
			 * well samples the profile. Otherwise there will be an error in the total
			 * feature strength of order 50%  */
			if( term1 >= -1. && term1 < -0.5 )
			{
				term = pah1_strength[j]/(3.*pah1_width[j]);
				term *= 2. + 2.*term1;
			}
			if( term1 >= -0.5 && term1 <= 1.5 )
				term = pah1_strength[j]/(3.*pah1_width[j]);
			if( term1 > 1.5 && term1 <= 3. )
			{
				term = pah1_strength[j]/(3.*pah1_width[j]);
				term *= 2. - term1*2./3.;
			}
		}
		else
		{
			/* This assumes linear interpolation between the points, which are
			 * located at -2, -1, +1, and +2 times the width, or a fine spacing that
			 * well samples the profile. Otherwise there will be an error in the total
			 * feature strength of order 50%  */
			if( term1 >= -2. && term1 < -1. )
			{
				term = pah1_strength[j]/(3.*pah1_width[j]);
				term *= 2. + term1;
			}
			if( term1 >= -1. && term1 <= 1. )
				term = pah1_strength[j]/(3.*pah1_width[j]);
			if( term1 > 1. && term1 <= 2. )
			{
				term = pah1_strength[j]/(3.*pah1_width[j]);
				term *= 2. - term1;
			}
		}
		if( j == 0 || j > 2 )
			term *= xnhoc;
		pah1_fun_v += term*xnc;
	}

	*cs_abs = pah1_fun_v;
	/* the next two numbers are completely arbitrary, but the code requires them... */
	/* >>chng 02 oct 18, cs_sct was 1.e-30, but this is very high for X-ray photons, PvH */
	*cs_sct = 0.1*pah1_fun_v;
	*cosb = 0.;
	*error = 0;

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


static void tbl_fun(double wavl,    /* in micron */
		    /*@in@*/ sd_data *sd, /* NOT USED -- MUST BE KEPT FOR COMPATIBILITY */
		    /*@in@*/ grain_data *gd,
		    /*@out@*/ double *cs_abs,
		    /*@out@*/ double *cs_sct,
		    /*@out@*/ double *cosb,
		    /*@out@*/ int *error)
{
	int lgOutOfBounds;
	long int ind;

#	ifdef DEBUG_FUN
	fputs( "<+>tbl_fun()\n", debug_fp );
#	endif
	double anu = WAVNRYD/wavl*1.e4;

	/* >>chng 02 nov 17, add this test to prevent warning that this var not used */
	if( sd == NULL )
		TotalInsanity();

	/*TODO	2	include code for interpolating inv_att_len somewhere!! */
	/*TODO	2	why is charge of pah2_ism_01.opc so different?? */

	find_arr(anu,gd->opcAnu,gd->nOpcData,&ind,&lgOutOfBounds);
	if( !lgOutOfBounds ) 
	{
		double a1g;
		double frac = log(anu/gd->opcAnu[ind])/log(gd->opcAnu[ind+1]/gd->opcAnu[ind]);

		*cs_abs = exp((1.-frac)*log(gd->opcData[0][ind])+frac*log(gd->opcData[0][ind+1]));
		ASSERT( *cs_abs > 0. );
		if( gd->nOpcCols > 1 )
			*cs_sct = exp((1.-frac)*log(gd->opcData[1][ind])+frac*log(gd->opcData[1][ind+1]));
		else
			*cs_sct = 0.1*(*cs_abs);
		ASSERT( *cs_sct > 0. );
		if( gd->nOpcCols > 2 )
			a1g = exp((1.-frac)*log(gd->opcData[2][ind])+frac*log(gd->opcData[2][ind+1]));
		else
			a1g = 1.;
		ASSERT( a1g > 0. );
		*cosb = 1. - a1g;
		*error = 0;
	}
	else
	{
		*cs_abs = -1.;
		*cs_sct = -1.;
		*cosb = -2.;
		*error = 3;
	}

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


static double size_distr(double size,
			 /*@in@*/ sd_data *sd)
{
	int lgOutOfBounds;
	long ind;
	double frac,
	  res,
	  x;

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

	if( size >= sd->lim[ipBLo] && size <= sd->lim[ipBHi] )
		switch( sd->sdCase ) 
		{
		case SD_SINGLE_SIZE:
			res = 1.; /* should really not be used in this case */
			break;
		case SD_POWERLAW:
			/* simple powerlaw */
		case SD_EXP_CUTOFF1:
		case SD_EXP_CUTOFF2:
		case SD_EXP_CUTOFF3:
			/* powerlaw with exponential cutoff, inspired by Greenberg (1978)
			 * Cosmic Dust, ed. J.A.M. McDonnell, Wiley, p. 187 */
			res = pow(size,sd->a[ipExp]);
			if( sd->a[ipBeta] < 0. )
				res /= (1. - sd->a[ipBeta]*size);
			else if( sd->a[ipBeta] > 0. )
				res *= (1. + sd->a[ipBeta]*size);
			if( size < sd->a[ipBLo] && sd->a[ipSLo] > 0. )
				res *= exp(-powi((sd->a[ipBLo]-size)/sd->a[ipSLo],NINT(sd->a[ipAlpha])));
			if( size > sd->a[ipBHi] && sd->a[ipSHi] > 0. )
				res *= exp(-powi((size-sd->a[ipBHi])/sd->a[ipSHi],NINT(sd->a[ipAlpha])));
			break;
		case SD_LOG_NORMAL:
			x = log(size/sd->a[ipGCen])/sd->a[ipGSig];
			res = exp(-0.5*POW2(x))/size;
			break;
		case SD_LIN_NORMAL:
			x = (size-sd->a[ipGCen])/sd->a[ipGSig];
			res = exp(-0.5*POW2(x))/size;
			break;
		case SD_TABLE:
			find_arr(log(size),sd->ln_a,sd->npts,&ind,&lgOutOfBounds);
			if( lgOutOfBounds )
			{
				fprintf( ioQQQ, " size distribution table has insufficient range\n" );
				fprintf( ioQQQ, " requested size: %.5f table range %.5f - %.5f\n",
					 size, exp(sd->ln_a[0]), exp(sd->ln_a[sd->npts-1]) );
				puts( "[Stop in size_distr]" );
				cdEXIT(EXIT_FAILURE);
			}				
			frac = (log(size)-sd->ln_a[ind])/(sd->ln_a[ind+1]-sd->ln_a[ind]);
			ASSERT( frac > 0.-10.*DBL_EPSILON && frac < 1.+10.*DBL_EPSILON );
			res = (1.-frac)*sd->ln_a4dNda[ind] + frac*sd->ln_a4dNda[ind+1];
			/* convert from a^4*dN/da to dN/da */
			res = exp(res)/POW4(size);
			break;
		case SD_ILLEGAL:
		default:
			fprintf( ioQQQ, " insane case for grain size distribution: %d\n" , sd->sdCase );
			ShowMe();
			puts( "[Stop in size_distr]" );
			cdEXIT(EXIT_FAILURE);
		}
	else
		res = 0.;

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


/* search for upper/lower limit lim of size distribution such that
 * lim^4 * dN/da(lim) < rel_cutoff * ref^4 * dN/da(ref)
 * the initial estimate of lim = ref + step
 * step may be both positive (upper limit) or negative (lower limit) */ 
static double search_limit(double ref,
			   double step,
			   double rel_cutoff,
			   sd_data sd)
{
	long i;
	double f1,
	  f2,
	  fmid,
	  renorm,
	  x1,
	  x2 = DBL_MAX,
	  xmid = DBL_MAX;

	/* TOLER is the relative accuracy with which lim is determined */
	const double TOLER = 1.e-6;

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

	/* sanity check */
	ASSERT( rel_cutoff > 0. && rel_cutoff < 1. );

	if( step == 0. )
	{
#		ifdef DEBUG_FUN
		fputs( " <->search_limit()\n", debug_fp );
#		endif
		return ref;
	}

	/* these need to be set in order for size_distr to work...
	 * NB - since this is a local copy of sd, it will not
	 * upset anything in the calling routine */
	sd.lim[ipBLo] = 0.;
	sd.lim[ipBHi] = DBL_MAX;

	x1 = ref;
	/* previous assert guarantees that f1 > 0. */
	f1 = -log(rel_cutoff);
	renorm = f1 - log(POW4(x1)*size_distr(x1,&sd));

	/* bracket solution */
	f2 = 1.;
	for( i=0; i < 20 && f2 > 0.; ++i )
	{
		x2 = MAX2(ref + step,SMALLEST_GRAIN);
		f2 = log(POW4(x2)*size_distr(x2,&sd)) + renorm;
		if( f2 >= 0. )
		{
			x1 = x2;
			f1 = f2;
		}
		step *= 2.;
	}
	if( f2 > 0. )
	{
		fprintf( ioQQQ, " Could not bracket solution\n" );
		puts( "[Stop in search_limit]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* do bisection search */
	while( 2.*fabs(x1-x2)/(x1+x2) > TOLER )
	{
		xmid = (x1+x2)/2.;
		fmid = log(POW4(xmid)*size_distr(xmid,&sd)) + renorm;

		if( fmid == 0. )
			break;

		if( f1*fmid > 0. )
		{
			x1 = xmid;
			f1 = fmid;
		}
		else
		{
			x2 = xmid;
			f2 = fmid;
		}
	}

#	ifdef DEBUG_FUN
	fputs( " <->search_limit()\n", debug_fp );
#	endif
	return (x1+x2)/2.;
}

/* calculate the inverse attenuation length for given refractive index data */
static void mie_calc_ial(/*@in@*/ grain_data *gd,
			 long int n,
			 /*@out@*/ double invlen[],  /* invlen[n] */
			 /*@in@*/ const char *string,
			 /*@in@*/ int *lgWarning)
{
	/* >>chng 01 aug 22, MALLOC this space */
	int *ErrorIndex/*[NC_ELL]*/,
	  lgErrorOccurred=TRUE,
	  lgOutOfBounds;
	long int i,
	  ind,
	  j;
	double frac,
	  InvDep,
	  nim,
	  wavlen;

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

	/* sanity check */
	ASSERT( gd->rfiType == RFI_TABLE );

	/* >>chng 01 aug 22, MALLOC this space */
	if( (ErrorIndex = (int*)MALLOC(sizeof(int)*(unsigned)rfield.nupper) ) == NULL )
		BadMalloc();

	for( i=0; i < n; i++ ) 
	{
		wavlen = WAVNRYD/rfield.anu[i]*1.e4;

		ErrorIndex[i] = 0;
		lgErrorOccurred = FALSE;
		invlen[i] = 0.;

		for( j=0; j < gd->nAxes; j++ ) 
		{
			/* first interpolate optical constants */
			find_arr(wavlen,gd->wavlen[j],gd->ndata[j],&ind,&lgOutOfBounds);
			if( lgOutOfBounds ) 
			{
				ErrorIndex[i] = 3;
				lgErrorOccurred = TRUE;
				invlen[i] = 0.;
				break;
			}
			frac = (wavlen-gd->wavlen[j][ind])/(gd->wavlen[j][ind+1]-gd->wavlen[j][ind]);
			nim = (1.-frac)*gd->n[j][ind].im + frac*gd->n[j][ind+1].im;
			/* this is the inverse of the photon attenuation depth,
			 * >>refer	grain	physics	Weingartner & Draine, 2000, ApJ, ... */
			InvDep = PI4*nim/wavlen*1.e4;
			ASSERT( InvDep > 0. );

			invlen[i] += InvDep*gd->wt[j];
		}
	}

	if( lgErrorOccurred ) 
	{
		mie_repair(string,n,3,3,rfield.anu,invlen,ErrorIndex,FALSE,lgWarning);
	}

	free( ErrorIndex );

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


/* this is the number of x-values we use for extrapolating functions */
#define NPTS_DERIV 8
#define NPTS_COMB  (NPTS_DERIV*(NPTS_DERIV-1)/2)

/* extrapolate/interpolate mie data to fill in the blanks */
static void mie_repair(/*@in@*/ const char *string,
		       long int n,
		       int val,
		       int del,
		       /*@in@*/ float anu[],      /* anu[n] */
		       double data[],             /* data[n] */
		       /*@in@*/ int ErrorIndex[], /* ErrorIndex[n] */
		       int lgRound,
		       /*@in@*/ int *lgWarning)
{
	int lgExtrapolate,
	  lgVerbose;
	long int i1,
	  i2,
	  ind1,
	  ind2,
	  j;
	double dx,
	  sgn,
	  slp1,
	  xlg1,
	  xlg2,
	  y1lg1,
	  y1lg2;

	/* interpolating over more that this number of points results in a warning */
	const long BIG_INTERPOLATION = 10;

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

	lgVerbose = ( string[0] != '\0' );

	for( ind1=0; ind1 < n; ) 
	{
		if( ErrorIndex[ind1] == val ) 
		{
			/* search for region with identical error index */
			ind2 = ind1;
			while( ind2 < n && ErrorIndex[ind2] == val )
				ind2++;

			if( lgVerbose )
				fprintf( ioQQQ, "    %s", string );

			if( ind1 == 0 ) 
			{
				/* low energy extrapolation */
				i1 = ind2;
				i2 = ind2+NPTS_DERIV-1;
				lgExtrapolate = TRUE;
				sgn = +1.;
				if( lgVerbose ) 
				{
					fprintf( ioQQQ, " extrapolated below %.4e Ryd\n",anu[i1] );
				}
			}
			else if( ind2 == n ) 
			{
				/* high energy extrapolation */
				i1 = ind1-NPTS_DERIV;
				i2 = ind1-1;
				lgExtrapolate = TRUE;
				sgn = -1.;
				if( lgVerbose ) 
				{
					fprintf( ioQQQ, " extrapolated above %.4e Ryd\n",anu[i2] );
				}
			}
			else 
			{
				/* interpolation */
				i1 = ind1-1;
				i2 = ind2;
				lgExtrapolate = FALSE;
				sgn = 0.;
				if( lgVerbose ) 
				{
					fprintf( ioQQQ, " interpolated between %.4e and %.4e Ryd\n",
						 anu[i1],anu[i2] );
				}
				if( i2-i1-1 > BIG_INTERPOLATION )
				{
					if( lgVerbose )
					{
						fprintf( ioQQQ, " ***Warning: extensive interpolation used\n" );
					}
					*lgWarning = TRUE;
				}
			}

			if( i1 < 0 || i2 >= n ) 
			{
				fprintf( ioQQQ, " Insufficient data for extrapolation\n" );
				puts( "[Stop in mie_repair]" );
				cdEXIT(EXIT_FAILURE);
			}

			xlg1 = log(anu[i1]);
			y1lg1 = log(data[i1]);
			/* >>chng 01 jul 30, replace simple-minded extrapolation with more robust routine, PvH */
			if( lgExtrapolate )
				slp1 = mie_find_slope(anu,data,ErrorIndex,i1,i2,val,lgVerbose,lgWarning);
			else
			{
				xlg2 = log(anu[i2]);
				y1lg2 = log(data[i2]);
				slp1 = (y1lg2-y1lg1)/(xlg2-xlg1);
			}
			if( lgRound && lgExtrapolate && sgn > 0. ) 
			{
				/* in low energy extrapolation, 1-g is very close to 1 and almost constant
				 * hence slp1 is very close to 0 and can even be slightly negative
				 * to prevent 1-g becoming greater than 1, the following is necessary */
				slp1 = MAX2(slp1,0.);
			}
			/* >>chng 02 oct 22, changed from sgn*slp1 <= 0. to accomodate grey grains */
			else if( lgExtrapolate && sgn*slp1 < 0. ) 
			{
				fprintf( ioQQQ, " Illegal value for slope in extrapolation %.6e\n", slp1 );
				puts( "[Stop in mie_repair]" );
				cdEXIT(EXIT_FAILURE);
			}

			for( j=ind1; j < ind2; j++ ) 
			{
				dx = log(anu[j]) - xlg1;
				data[j] = exp(y1lg1 + dx*slp1);
				ErrorIndex[j] -= del;
			}

			ind1 = ind2;
		}
		else 
		{
			ind1++;
		}
	}
	/* sanity check */
	for( j=0; j < n; j++ )
	{
		if( ErrorIndex[j] > val-del ) 
		{
			fprintf( ioQQQ, " Internal error in mie_repair, index=%ld, val=%d\n",j,ErrorIndex[j] );
			ShowMe();
			puts( "[Stop in mie_repair]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

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


static double mie_find_slope(/*@in@*/ const float anu[],
			     /*@in@*/ const double data[],
			     /*@in@*/ const int ErrorIndex[],
			     long i1,
			     long i2,
			     int val,
			     int lgVerbose,
			     /*@in@*/ int *lgWarning)
{
	long i,
	  j,
	  k;
	double s1,
	  s2,
	  slope,
	  slp1[NPTS_COMB],
	  stdev;

	/* threshold for standard deviation in the logarithmic derivative to generate warning,
	 * this corresponds to an uncertainty of a factor 10 for a typical extrapolation */
	const double LARGE_STDEV = 0.2;

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

	/* sanity check */
	ASSERT( i2-i1 == NPTS_DERIV-1 );
	for( i=i1; i <= i2; i++ )
	{
		ASSERT( ErrorIndex[i] < val );
		ASSERT( anu[i] > 0. && data[i] > 0. );
	}

	/* this initialization is to keep lint happy, the next statement will do the real initialization */
	for( i=0; i < NPTS_COMB; i++ )
	{
		slp1[i] = -DBL_MAX;
	}

	k = 0;
	/* calculate the logarithmic derivative for every possible combination of data points */
	/* this loop is guaranteed to initialize all members of slp1, the first ASSERT checks this */
	for( i=i1; i < i2; i++ )
	{
		for( j=i+1; j <= i2; j++ )
		{
			slp1[k++] = log(data[j]/data[i])/log(anu[j]/anu[i]);
		}
	}
	/* sort the values; we want the median -> values for i > NPTS_COMB/2 are irrelevant */
	for( i=0; i <= NPTS_COMB/2; i++ )
	{
		for( j=i+1; j < NPTS_COMB; j++ )
		{
			if( slp1[i] > slp1[j] )
			{
				double xxx = slp1[i];
				slp1[i] = slp1[j];
				slp1[j] = xxx;
			}
		}
	}
	/* now calculate the median value */
	slope = ( NPTS_COMB%2 == 1 ) ? slp1[NPTS_COMB/2] : (slp1[NPTS_COMB/2-1] + slp1[NPTS_COMB/2])/2.;

	/* and finally calculate the standard deviation of all slopes */
	s1 = s2 = 0.;
	for( i=0; i < NPTS_COMB; i++ )
	{
		s1 += slp1[i];
		s2 += POW2(slp1[i]);
	}
	stdev = sqrt(s2/(double)NPTS_COMB - POW2(s1/(double)NPTS_COMB));

#if 0
	for( i=i1; i <= i2; i++ )
		printf("input: %ld %.4e %.4e\n",i,anu[i],data[i]);
	for( i=0; i < NPTS_COMB; i++ )
		printf("%.3f ",slp1[i]);
	printf("\n");
	printf("slope %.3f +/- %.3f\n",slope,stdev);
#endif
	
	/* print warning if standard deviation is large */
	if( stdev > LARGE_STDEV )
	{
		if( lgVerbose )
		{
			fprintf( ioQQQ, " ***Warning: slope for extrapolation may be unreliable\n" );
		}
		*lgWarning = TRUE;
	}

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


/* read in the file with optical constants and other relevant information */
static void mie_read_rfi(/*@in@*/  const char *chFile,
			 /*@out@*/ grain_data *gd)
{
	int lgLogData=FALSE;
	long int dl,
	  help,
	  i,
	  nelem,
	  j,
	  nridf,
	  sgn = 0;
	double eps1,
	  eps2,
	  LargestLog,
	  molw,
	  nAtoms,
	  nr,
	  ni,
	  tmp1,
	  tmp2,
	  total = 0.;
	char chLine[FILENAME_PATH_LENGTH_2],
	  chWord[FILENAME_PATH_LENGTH_2];
	FILE *io2;

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

	if( (io2 = fopen(chFile,"r")) == NULL )
	{
		fprintf( ioQQQ, " Could not open %s for reading\n",chFile);
		puts( "[Stop in mie_read_rfi]" );
		cdEXIT(EXIT_FAILURE);
	}

	dl = 0; /* line counter for input file */

	/* first read magic number */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_long(chFile,chLine,&gd->magic,TRUE,dl);
	if( gd->magic != MAGIC_RFI ) 
	{
		fprintf( ioQQQ, " Refractive index file %s has obsolete magic number\n",chFile );
		fprintf( ioQQQ, " I found %ld, but expected %ld on line #%ld\n",gd->magic,MAGIC_RFI,dl );
		fprintf( ioQQQ, " Please replace this file with an up to date version\n" );
		puts( "[Stop in mie_read_rfi]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* get chemical formula of the grain, e.g., Mg0.4Fe0.6SiO3 */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_word(chLine,chWord,FILENAME_PATH_LENGTH_2,FALSE);
	mie_read_form(chWord,gd->elmAbun,&nAtoms,&molw);

	/* molecular weight, in atomic units */
	gd->mol_weight = molw;
	gd->atom_weight = gd->mol_weight/nAtoms;

	/* determine abundance of grain molecule assuming max depletion */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_double(chFile,chLine,&gd->abun,TRUE,dl);

	/* default depletion of grain molecule */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_double(chFile,chLine,&gd->depl,TRUE,dl);
	if( gd->depl > 1. ) 
	{
		fprintf( ioQQQ, " Illegal value for default depletion in %s\n",chFile );
		fprintf( ioQQQ, " Line #%ld, depl=%14.6e\n",dl,gd->depl);
		puts( "[Stop in mie_read_rfi]" );
		cdEXIT(EXIT_FAILURE);
	}

	for( nelem=0; nelem < LIMELM; nelem++ )
		gd->elmAbun[nelem] *= gd->abun*gd->depl;

	/* material density, to get cross section per unit mass: rho in cgs */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_double(chFile,chLine,&gd->rho,FALSE,dl);

	/* material type, determines enthalpy function: 1 -- carbonaceous, 2 -- silicate */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_long(chFile,chLine,&help,TRUE,dl);
	gd->matType = (mat_type)help;
	if( gd->matType >= MAT_TOP ) 
	{
		fprintf( ioQQQ, " Illegal value for material type in %s\n",chFile );
		fprintf( ioQQQ, " Line #%ld, type=%d\n",dl,gd->matType);
		puts( "[Stop in mie_read_rfi]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* work function, in Rydberg */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_double(chFile,chLine,&gd->work,TRUE,dl);

	/* bandgap, in Rydberg */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_double(chFile,chLine,&gd->bandgap,FALSE,dl);
	if( gd->bandgap >= gd->work ) 
	{
		fprintf( ioQQQ, " Illegal value for bandgap in %s\n",chFile );
		fprintf( ioQQQ, " Line #%ld, bandgap=%.4e, work function=%.4e\n",dl,gd->bandgap,gd->work);
		fprintf( ioQQQ, " Bandgap should always be less than work function\n");
		puts( "[Stop in mie_read_rfi]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* efficiency of thermionic emission, between 0 and 1 */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_double(chFile,chLine,&gd->therm_eff,TRUE,dl);
	if( gd->therm_eff > 1.f ) 
	{
		fprintf( ioQQQ, " Illegal value for thermionic efficiency in %s\n",chFile );
		fprintf( ioQQQ, " Line #%ld, value=%.4e\n",dl,gd->therm_eff);
		fprintf( ioQQQ, " Allowed values are 0. < efficiency <= 1.\n");
		puts( "[Stop in mie_read_rfi]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* sublimation temperature in K */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_double(chFile,chLine,&gd->subl_temp,TRUE,dl);

	/* >>chng 02 sep 18, add keyword for special files (grey grains, PAH's, etc.), PvH */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_word(chLine,chWord,WORDLEN,TRUE);

	if( lgMatch( chWord, "RFI_" ) )
		gd->rfiType = RFI_TABLE;
	else if( lgMatch( chWord, "OPC_" ) )
		gd->rfiType = OPC_TABLE;
	else if( lgMatch( chWord, "GREY" ) )
		gd->rfiType = OPC_GREY;
	else if( lgMatch( chWord, "PAH1" ) )
		gd->rfiType = OPC_PAH1;
	else
	{
		fprintf( ioQQQ, " Illegal keyword in %s\n",chFile );
		fprintf( ioQQQ, " Line #%ld, value=%s\n",dl,chWord);
		fprintf( ioQQQ, " Allowed values are: RFI_TBL, OPC_TBL, GREY, PAH1\n");
		puts( "[Stop in mie_read_rfi]" );
		cdEXIT(EXIT_FAILURE);
	}

	switch( gd->rfiType )
	{
	case RFI_TABLE:
		/* nridf is for choosing ref index or diel function input
		 * case 2 allows greater accuracy reading in, when nr is close to 1. */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_long(chFile,chLine,&nridf,TRUE,dl);
		if( nridf > 3 ) 
		{
			fprintf( ioQQQ, " Illegal data code in %s\n",chFile );
			fprintf( ioQQQ, " Line #%ld, data code=%ld\n",dl,nridf);
			puts( "[Stop in mie_read_rfi]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* no. of principal axes, always 1 for amorphous grains,
		 * maybe larger for crystalline grains */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_long(chFile,chLine,&gd->nAxes,TRUE,dl);
		if( gd->nAxes > NAX ) 
		{
			fprintf( ioQQQ, " Illegal no. of axes in %s\n",chFile );
			fprintf( ioQQQ, " Line #%ld, number=%ld\n",dl,gd->nAxes);
			puts( "[Stop in mie_read_rfi]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* now get relative weights of axes */
		mie_next_data(chFile,io2,chLine,&dl);
		switch( gd->nAxes ) 
		{
		case 1:
			mie_read_double(chFile,chLine,&gd->wt[0],TRUE,dl);
			total = gd->wt[0];
			break;
		case 2:
			if( sscanf( chLine, "%lf %lf", &gd->wt[0], &gd->wt[1] ) != 2 ) 
			{
				fprintf( ioQQQ, " Syntax error in %s\n",chFile);
				fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
				puts( "[Stop in mie_read_rfi]" );
				cdEXIT(EXIT_FAILURE);
			}
			if( gd->wt[0] <= 0. || gd->wt[1] <= 0. ) 
			{
				fprintf( ioQQQ, " Illegal data in %s\n",chFile);
				fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
				puts( "[Stop in mie_read_rfi]" );
				cdEXIT(EXIT_FAILURE);
			}
			total = gd->wt[0] + gd->wt[1];
			break;
		case 3:
			if( sscanf( chLine, "%lf %lf %lf", &gd->wt[0], &gd->wt[1], &gd->wt[2] ) != 3 ) 
			{
				fprintf( ioQQQ, " Syntax error in %s\n",chFile);
				fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
				puts( "[Stop in mie_read_rfi]" );
				cdEXIT(EXIT_FAILURE);
			}
			if( gd->wt[0] <= 0. || gd->wt[1] <= 0. || gd->wt[2] <= 0. ) 
			{
				fprintf( ioQQQ, " Illegal data in %s\n",chFile);
				fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
				puts( "[Stop in mie_read_rfi]" );
				cdEXIT(EXIT_FAILURE);
			}
			total = gd->wt[0] + gd->wt[1] + gd->wt[2];
			break;
		default:
			fprintf( ioQQQ, " insane case for gd->nAxes: %ld\n", gd->nAxes );
			ShowMe();
			puts( "[Stop in mie_read_rfi]" );
			cdEXIT(EXIT_FAILURE);
		}
		for( j=0; j < gd->nAxes; j++ ) 
		{
			gd->wt[j] /= total;

			/* read in optical constants for each principal axis. */
			mie_next_data(chFile,io2,chLine,&dl);
			mie_read_long(chFile,chLine,&gd->ndata[j],FALSE,dl);
			if ( gd->ndata[j] < 2 ) 
			{
				fprintf( ioQQQ, " Illegal number of data points in %s\n",chFile );
				fprintf( ioQQQ, " Line #%ld, number=%ld\n",dl,gd->ndata[j]);
				puts( "[Stop in mie_read_rfi]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* allocate space for wavelength and optical constants arrays
			 * the memory is freed in mie_write_opc */
			if( ( gd->wavlen[j] = (double*)MALLOC(sizeof(double)*(unsigned)gd->ndata[j]) ) == NULL )
				BadMalloc();
			if( ( gd->n[j] = (complex*)MALLOC(sizeof(complex)*(unsigned)gd->ndata[j]) ) == NULL )
				BadMalloc();
			if( ( gd->nr1[j] = (double*)MALLOC(sizeof(double)*(unsigned)gd->ndata[j]) ) == NULL )
				BadMalloc();

			for( i=0; i < gd->ndata[j]; i++ ) 
			{
				/* read in the wavelength in microns
				 * and the complex refractive index or dielectric function of material */
				mie_next_data(chFile,io2,chLine,&dl);
				if( sscanf( chLine, "%lf %lf %lf", gd->wavlen[j]+i, &nr, &ni ) != 3 ) 
				{
					fprintf( ioQQQ, " Syntax error in %s\n",chFile);
					fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
					puts( "[Stop in mie_read_rfi]" );
					cdEXIT(EXIT_FAILURE);
				}
				if( gd->wavlen[j][i] <= 0. ) 
				{
					fprintf( ioQQQ, " Illegal value for wavelength in %s\n",chFile);
					fprintf( ioQQQ, " Line #%ld, wavl=%14.6e\n",dl,gd->wavlen[j][i]);
					puts( "[Stop in mie_read_rfi]" );
					cdEXIT(EXIT_FAILURE);
				}
				/* the data in the input file should be sorted on wavelength, either
				 * strictly monotonically increasing or decreasing, check this here... */
				if( i == 1 )
				{
					sgn = SIGN3(gd->wavlen[j][1]-gd->wavlen[j][0]);
					if( sgn == 0 ) 
					{
						fprintf( ioQQQ, " Illegal value for wavelength in %s\n",chFile);
						fprintf( ioQQQ, " Line #%ld, wavl=%14.6e\n",dl,gd->wavlen[j][i]);
						puts( "[Stop in mie_read_rfi]" );
						cdEXIT(EXIT_FAILURE);
					}
				}
				else if( i > 1 ) 
				{
					if( SIGN3(gd->wavlen[j][i]-gd->wavlen[j][i-1]) != sgn ) 
					{
						fprintf( ioQQQ, " Illegal value for wavelength in %s\n",chFile);
						fprintf( ioQQQ, " Line #%ld, wavl=%14.6e\n",dl,gd->wavlen[j][i]);
						puts( "[Stop in mie_read_rfi]" );
						cdEXIT(EXIT_FAILURE);
					}
				}
				/* this version reads in real and imaginary parts of the refractive
				 * index, with imaginary part positive (nridf = 3) or nr-1 (nridf = 2) or
				 * real and imaginary parts of the dielectric function (nridf = 1) */
				switch( nridf ) 
				{
				case 1:
					eps1 = nr;
					eps2 = ni;
					dftori(&nr,&ni,eps1,eps2);
					gd->nr1[j][i] = nr - 1.;
					break;
				case 2:
					gd->nr1[j][i] = nr;
					nr += 1.;
					break;
				case 3:
					gd->nr1[j][i] = nr - 1.;
					break;
				default:
					fprintf( ioQQQ, " insane case for nridf: %ld\n", nridf );
					ShowMe();
					puts( "[Stop in mie_read_rfi]" );
					cdEXIT(EXIT_FAILURE);
				}
				gd->n[j][i] = ftocm(nr,ni);

				/* sanity checks */
				if( nr <= 0. || ni < 0. ) 
				{
					fprintf( ioQQQ, " Illegal value for refractive index in %s\n",chFile);
					fprintf( ioQQQ, " Line #%ld, (nr,ni)=(%14.6e,%14.6e)\n",dl,nr,ni);
					puts( "[Stop in mie_read_rfi]" );
					cdEXIT(EXIT_FAILURE);
				}
				ASSERT( fabs(nr-1.-gd->nr1[j][i]) < 10.*nr*DBL_EPSILON );
			}
		}
		break;
	case OPC_TABLE:
		/* no. of data columns in OPC_TABLE file:
		 * 1: absorption cross sections only
		 * 2: absorption + scattering cross sections
		 * 3: absorption + pure scattering cross sections + asymmetry factor
		 * 4: absorption + pure scattering cross sections + asymmetry factor +
		 *      inverse attenuation length */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_long(chFile,chLine,&gd->nOpcCols,TRUE,dl);
		if( gd->nOpcCols > NDAT ) 
		{
			fprintf( ioQQQ, " Illegal no. of data columns in %s\n",chFile );
			fprintf( ioQQQ, " Line #%ld, number=%ld\n",dl,gd->nOpcCols);
			puts( "[Stop in mie_read_rfi]" );
			cdEXIT(EXIT_FAILURE);
		}
		
		/* keyword indicating whether the table contains linear or logarithmic data */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_word(chLine,chWord,WORDLEN,TRUE);

		if( lgMatch( chWord, "LINE" ) )
			lgLogData = FALSE;
		else if( lgMatch( chWord, "LOG" ) )
			lgLogData = TRUE;
		else
		{
			fprintf( ioQQQ, " Keyword not recognized in %s\n",chFile );
			fprintf( ioQQQ, " Line #%ld, keyword=%s\n",dl,chWord);
			puts( "[Stop in mie_read_rfi]" );
			cdEXIT(EXIT_FAILURE);
		}
			

		/* read in number of data points supplied. */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_long(chFile,chLine,&gd->nOpcData,FALSE,dl);
		if ( gd->nOpcData < 2 ) 
		{
			fprintf( ioQQQ, " Illegal number of data points in %s\n",chFile );
			fprintf( ioQQQ, " Line #%ld, number=%ld\n",dl,gd->nOpcData);
			puts( "[Stop in mie_read_rfi]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* allocate space for frequency and data arrays
		 * the memory is freed in mie_write_opc */
		if( ( gd->opcAnu = (double*)MALLOC(sizeof(double)*(unsigned)gd->nOpcData) ) == NULL )
			BadMalloc();
		for( j=0; j < gd->nOpcCols; j++ ) 
		{
			if( ( gd->opcData[j] = (double*)MALLOC(sizeof(double)*(unsigned)gd->nOpcData) ) == NULL )
				BadMalloc();
		}

		tmp1 = -log10(1.01*DBL_MIN);
		tmp2 = log10(0.99*DBL_MAX);
		LargestLog = MIN2(tmp1,tmp2);

		/* now read the data... each line should contain:
		 *
		 * if gd->nOpcCols == 1, anu, abs_cs
		 * if gd->nOpcCols == 2, anu, abs_cs, sct_cs
		 * if gd->nOpcCols == 3, anu, abs_cs, sct_cs, inv_att_len
		 * 
		 * the frequencies in the table should be either monotonically increasing or decreasing.
		 * frequencies should be in Ryd, cross sections in cm^2/H, and the inverse attenuation length
		 * in cm^-1. If lgLogData is TRUE, each number should be the log10 of the data value */
		for( i=0; i < gd->nOpcData; i++ ) 
		{
			mie_next_data(chFile,io2,chLine,&dl);
			switch( gd->nOpcCols )
			{
			case 1:
				if( sscanf( chLine, "%lf %lf", &gd->opcAnu[i], &gd->opcData[0][i] ) != 2 ) 
				{
					fprintf( ioQQQ, " Syntax error in %s\n",chFile);
					fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
					puts( "[Stop in mie_read_rfi]" );
					cdEXIT(EXIT_FAILURE);
				}
				break;
			case 2:
				if( sscanf( chLine, "%lf %lf %lf", &gd->opcAnu[i], &gd->opcData[0][i],
					    &gd->opcData[1][i] ) != 3 ) 
				{
					fprintf( ioQQQ, " Syntax error in %s\n",chFile);
					fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
					puts( "[Stop in mie_read_rfi]" );
					cdEXIT(EXIT_FAILURE);
				}
				break;
			case 3:
				if( sscanf( chLine, "%lf %lf %lf %lf", &gd->opcAnu[i], &gd->opcData[0][i],
					    &gd->opcData[1][i], &gd->opcData[2][i] ) != 4 ) 
				{
					fprintf( ioQQQ, " Syntax error in %s\n",chFile);
					fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
					puts( "[Stop in mie_read_rfi]" );
					cdEXIT(EXIT_FAILURE);
				}
				break;
			case 4:
				if( sscanf( chLine, "%lf %lf %lf %lf %lf", &gd->opcAnu[i], &gd->opcData[0][i],
					    &gd->opcData[1][i], &gd->opcData[2][i], &gd->opcData[3][i] ) != 5 ) 
				{
					fprintf( ioQQQ, " Syntax error in %s\n",chFile);
					fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
					puts( "[Stop in mie_read_rfi]" );
					cdEXIT(EXIT_FAILURE);
				}
				break;
			default:
				fprintf( ioQQQ, " insane case for gd->nOpcCols: %ld\n", gd->nOpcCols );
				ShowMe();
				puts( "[Stop in mie_read_rfi]" );
				cdEXIT(EXIT_FAILURE);
			}
			if( lgLogData )
			{
				/* this test will guarantee there will be neither under- nor overflows */
				if( fabs(gd->opcAnu[i]) > LargestLog )
				{
					fprintf( ioQQQ, " Illegal value for frequency in %s\n",chFile );
					fprintf( ioQQQ, " Line #%ld, freq=%14.6e\n",dl,gd->opcAnu[i] );
					puts( "[Stop in mie_read_rfi]" );
					cdEXIT(EXIT_FAILURE);
				}
				gd->opcAnu[i] = pow(10.,gd->opcAnu[i]);
				for( j=0; j < gd->nOpcCols; j++ )
				{
					if( fabs(gd->opcData[j][i]) > LargestLog )
					{
						fprintf( ioQQQ, " Illegal data value in %s\n",chFile );
						fprintf( ioQQQ, " Line #%ld, value=%14.6e\n",dl,gd->opcData[j][i] );
						puts( "[Stop in mie_read_rfi]" );
						cdEXIT(EXIT_FAILURE);
					}
					gd->opcData[j][i] = pow(10.,gd->opcData[j][i]);
				}
			}
			if( gd->opcAnu[i] <= 0. )
			{
				fprintf( ioQQQ, " Illegal value for frequency in %s\n",chFile );
				fprintf( ioQQQ, " Line #%ld, freq=%14.6e\n",dl,gd->opcAnu[i] );
				puts( "[Stop in mie_read_rfi]" );
				cdEXIT(EXIT_FAILURE);
			}
			for( j=0; j < gd->nOpcCols; j++ )
			{
				if( gd->opcData[j][i] <= 0. )
				{
					fprintf( ioQQQ, " Illegal data value in %s\n",chFile );
					fprintf( ioQQQ, " Line #%ld, value=%14.6e\n",dl,gd->opcData[j][i] );
					puts( "[Stop in mie_read_rfi]" );
					cdEXIT(EXIT_FAILURE);
				}
			}
			/* the data in the input file should be sorted on frequency, either
			 * strictly monotonically increasing or decreasing, check this here... */
			if( i == 1 )
			{
				sgn = SIGN3(gd->opcAnu[1]-gd->opcAnu[0]);
				if( sgn == 0 ) 
				{
					double dataVal = lgLogData ? log10(gd->opcAnu[1]) : gd->opcAnu[1];
					fprintf( ioQQQ, " Illegal value for frequency in %s\n",chFile );
					fprintf( ioQQQ, " Line #%ld, freq=%14.6e\n",dl,dataVal );
					puts( "[Stop in mie_read_rfi]" );
					cdEXIT(EXIT_FAILURE);
				}
			}
			else if( i > 1 ) 
			{
				if( SIGN3(gd->opcAnu[i]-gd->opcAnu[i-1]) != sgn ) 
				{
					double dataVal = lgLogData ? log10(gd->opcAnu[i]) : gd->opcAnu[i];
					fprintf( ioQQQ, " Illegal value for frequency in %s\n",chFile);
					fprintf( ioQQQ, " Line #%ld, freq=%14.6e\n",dl,dataVal);
					puts( "[Stop in mie_read_rfi]" );
					cdEXIT(EXIT_FAILURE);
				}
			}
		}
		gd->nAxes = 1;
		break;
	case OPC_GREY:
	case OPC_PAH1:
		/* nothing much to be done here, the opacities
		 * will be calculated without any further data */
		gd->nAxes = 1;
		break;
	default:
		fprintf( ioQQQ, " Insane value for gd->rfiType: %d, bailing out\n", gd->rfiType );
		ShowMe();
		puts( "[Stop in mie_read_rfi]" );
		cdEXIT(EXIT_FAILURE);
	}

	fclose(io2);

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


/* construct optical constants for mixed grain using a specific EMT */
static void mie_read_mix(/*@in@*/  const char *chFile,
			 /*@out@*/ grain_data *gd)
{
	emt_type EMTtype;
	long int dl,
	  i,
	  j,
	  k,
	  l,
	  nelem,
	  nMaterial,
	  sumAxes;
	double *delta,
	  *frac,
	  *frac2,
	  *frdelta,
	  maxIndex = DBL_MAX,
	  minIndex = DBL_MAX,
	  nAtoms,
	  sum,
	  sum2,
	  wavHi,
	  wavLo,
	  wavStep;
	complex *eps,
	  eps_eff = {-DBL_MAX,-DBL_MAX};
	char chLine[FILENAME_PATH_LENGTH_2],
	  chWord[FILENAME_PATH_LENGTH_2],
	  *str;
	grain_data *gdArr;
	FILE *io2;

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

	if( (io2 = fopen(chFile,"r")) == NULL )
	{
		fprintf( ioQQQ, " Could not open %s for reading\n",chFile);
		puts( "[Stop in mie_read_mix]" );
		cdEXIT(EXIT_FAILURE);
	}

	dl = 0; /* line counter for input file */

	/* first read magic number */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_long(chFile,chLine,&gd->magic,TRUE,dl);
	if( gd->magic != MAGIC_MIX ) 
	{
		fprintf( ioQQQ, " Mixed grain file %s has obsolete magic number\n",chFile );
		fprintf( ioQQQ, " I found %ld, but expected %ld on line #%ld\n",gd->magic,MAGIC_MIX,dl );
		fprintf( ioQQQ, " Please replace this file with an up to date version\n" );
		puts( "[Stop in mie_read_mix]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* default depletion of grain molecule */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_double(chFile,chLine,&gd->depl,TRUE,dl);
	if( gd->depl > 1. ) 
	{
		fprintf( ioQQQ, " Illegal value for default depletion in %s\n",chFile );
		fprintf( ioQQQ, " Line #%ld, depl=%14.6e\n",dl,gd->depl);
		puts( "[Stop in mie_read_mix]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* read number of different materials contained in this grain */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_long(chFile,chLine,&nMaterial,TRUE,dl);
	if( nMaterial < 2 ) 
	{
		fprintf( ioQQQ, " Illegal number of materials in mixed grain file %s\n",chFile );
		fprintf( ioQQQ, " I found %ld on line #%ld\n",nMaterial,dl );
		fprintf( ioQQQ, " This number should be at least 2\n" );
		puts( "[Stop in mie_read_mix]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( ( frac = (double*)MALLOC(sizeof(double)*(unsigned)nMaterial) ) == NULL )
		BadMalloc();
	if( ( frac2 = (double*)MALLOC(sizeof(double)*(unsigned)nMaterial) ) == NULL )
		BadMalloc();
	if( ( gdArr = (grain_data*)MALLOC(sizeof(grain_data)*(unsigned)nMaterial) ) == NULL )
		BadMalloc();

	sum = 0.;
	sum2 = 0.;
	sumAxes = 0;
	for( i=0; i < nMaterial; i++ )
	{
		char chFile2[FILENAME_PATH_LENGTH_2];

		gdArr[i].nAxes = 0;
		for( j=0; j < NAX; j++ ) 
		{
			gdArr[i].wavlen[j] = NULL;
			gdArr[i].n[j] = NULL;
			gdArr[i].nr1[j] = NULL;
		}
		gdArr[i].nOpcCols = 0;
		gdArr[i].opcAnu = NULL;
		for( j=0; j < NDAT; j++ )
		{
			gdArr[i].opcData[j] = NULL;
		}

		/* each line contains relative fraction of volume occupied by each material,
		 * followed by the name of the refractive index file between double quotes */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_double(chFile,chLine,&frac[i],TRUE,dl);

		sum += frac[i];

		str = strchr( chLine, '\"' );
		if( str != NULL )
		{
			strcpy( chFile2, ++str );
			str = strchr( chFile2, '\"' );
			if( str != NULL )
				*str = '\0';
		}
		if( str == NULL )
		{
			fprintf( ioQQQ, " No pair of double quotes was found on line #%ld of file %s\n",dl,chFile );
			fprintf( ioQQQ, " Please supply the refractive index file name between double quotes\n" );
			puts( "[Stop in mie_read_mix]" );
			cdEXIT(EXIT_FAILURE);
		}

		mie_read_rfi( chFile2, &gdArr[i] );
		if( gdArr[i].rfiType != RFI_TABLE )
		{
			fprintf( ioQQQ, " Input error on line #%ld of file %s\n",dl,chFile );
			fprintf( ioQQQ, " File %s is not of type RFI_TBL, this is illegal\n",chFile2 );
			puts( "[Stop in mie_read_mix]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* this test also guarantees that sum2 > 0 */
		if( i == nMaterial-1 && gdArr[i].mol_weight == 0. )
		{
			fprintf( ioQQQ, " Input error on line #%ld of file %s\n",dl,chFile );
			fprintf( ioQQQ, " Last entry in list of materials is vacuum, this is not allowed\n" );
			fprintf( ioQQQ, " Please move this entry to an earlier position in the file\n" );
			puts( "[Stop in mie_read_mix]" );
			cdEXIT(EXIT_FAILURE);
		}

		frac2[i] = ( gdArr[i].mol_weight > 0. ) ? frac[i] : 0.;
		sum2 += frac2[i];

		sumAxes += gdArr[i].nAxes;
	}

	/* read keyword to determine which EMT to use */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_word(chLine,chWord,WORDLEN,TRUE);

	if( lgMatch( chWord, "FA00" ) )
		EMTtype = FARAFONOV00;
	else if( lgMatch( chWord, "ST95" ) )
		EMTtype = STOGNIENKO95;
	else if( lgMatch( chWord, "BR35" ) )
		EMTtype = BRUGGEMAN35;
	else
	{
		fprintf( ioQQQ, " Keyword not recognized in %s\n",chFile );
		fprintf( ioQQQ, " Line #%ld, keyword=%s\n",dl,chWord);
		puts( "[Stop in mie_read_mix]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* normalize sum of fractional volumes to 1 */
	for( i=0; i < nMaterial; i++ )
	{
		frac[i] /= sum;
		frac2[i] /= sum2;
		/* renormalize elmAbun to chemical formula */
		for( nelem=0; nelem < LIMELM; nelem++ )
		{
			gdArr[i].elmAbun[nelem] /= gdArr[i].abun*gdArr[i].depl;
		}
	}

	wavLo = 0.;
	wavHi = DBL_MAX;
	gd->abun = DBL_MAX;
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		gd->elmAbun[nelem] = 0.;
	}
	gd->mol_weight = 0.;
	gd->rho = 0.;
	gd->work = DBL_MAX;
	gd->bandgap = DBL_MAX;
	gd->therm_eff = 0.;
	gd->subl_temp = DBL_MAX;
	gd->nAxes = 1;
	gd->wt[0] = 1.;
	gd->ndata[0] = MIX_TABLE_SIZE;
	gd->rfiType = RFI_TABLE;

	for( i=0; i < nMaterial; i++ )
	{
		for( k=0; k < gdArr[i].nAxes; k++ )
		{
			double wavMin = MIN2(gdArr[i].wavlen[k][0],gdArr[i].wavlen[k][gdArr[i].ndata[k]-1]);
			double wavMax = MAX2(gdArr[i].wavlen[k][0],gdArr[i].wavlen[k][gdArr[i].ndata[k]-1]);
			wavLo = MAX2(wavLo,wavMin);
			wavHi = MIN2(wavHi,wavMax);
		}
		minIndex = DBL_MAX;
		maxIndex = 0.;
		for( nelem=0; nelem < LIMELM; nelem++ )
		{
			gd->elmAbun[nelem] += frac[i]*gdArr[i].elmAbun[nelem];
			if( gd->elmAbun[nelem] > 0. )
			{
				minIndex = MIN2(minIndex,gd->elmAbun[nelem]);
			}
			maxIndex = MAX2(maxIndex,gd->elmAbun[nelem]);
		}
		gd->mol_weight += frac[i]*gdArr[i].mol_weight;
		gd->rho += frac[i]*gdArr[i].rho;
		/* ignore parameters for vacuum */
		if( gdArr[i].mol_weight > 0. )
		{
			gd->abun = MIN2(gd->abun,gdArr[i].abun/frac[i]);
			switch( EMTtype )
			{
			case FARAFONOV00:
				/* this is appropriate for a layered grain */
				gd->work = gdArr[i].work;
				gd->bandgap = gdArr[i].bandgap;
				gd->therm_eff = gdArr[i].therm_eff;
				break;
			case STOGNIENKO95:
			case BRUGGEMAN35:
				/* this is appropriate for a randomly mixed grain */
				gd->work = MIN2(gd->work,gdArr[i].work);
				gd->bandgap = MIN2(gd->bandgap,gdArr[i].bandgap);
				gd->therm_eff += frac2[i]*gdArr[i].therm_eff;
				break;
			default:
				fprintf( ioQQQ, " Insanity in mie_read_mix\n" );
				ShowMe();
				puts( "[Stop in mie_read_mix]" );
				cdEXIT(EXIT_FAILURE);
			}
			gd->matType = gdArr[i].matType;
			gd->subl_temp = MIN2(gd->subl_temp,gdArr[i].subl_temp);
		}
	}

	if( gd->rho <= 0. )
	{
		fprintf( ioQQQ, " Illegal value for the density: %.3e\n", gd->rho );
		puts( "[Stop in mie_read_mix]" );
		cdEXIT(EXIT_FAILURE);
	}
	if( gd->mol_weight <= 0. )
	{
		fprintf( ioQQQ, " Illegal value for the molecular weight: %.3e\n", gd->mol_weight );
		puts( "[Stop in mie_read_mix]" );
		cdEXIT(EXIT_FAILURE);
	}
	if( maxIndex <= 0. )
	{
		fprintf( ioQQQ, " No atoms were found in the grain molecule\n" );
		puts( "[Stop in mie_read_mix]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* further sanity checks */
	ASSERT( wavLo > 0. && wavHi < DBL_MAX && wavLo < wavHi );
	ASSERT( gd->abun > 0. && gd->abun < DBL_MAX );
	ASSERT( gd->work > 0. && gd->work < DBL_MAX );
	ASSERT( gd->bandgap >= 0. && gd->bandgap < gd->work );
	ASSERT( gd->therm_eff > 0. && gd->therm_eff <= 1. );
	ASSERT( gd->subl_temp > 0. && gd->subl_temp < DBL_MAX );
	ASSERT( minIndex > 0. && minIndex < DBL_MAX );

	/* apply safety margin */
	wavLo *= 1. + 10.*DBL_EPSILON;
	wavHi *= 1. - 10.*DBL_EPSILON;

	/* renormalize the chemical formula such that the lowest index is 1 */
	nAtoms = 0.;
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		gd->elmAbun[nelem] /= minIndex;
		nAtoms += gd->elmAbun[nelem];
	}
	ASSERT( nAtoms > 0. );
	gd->abun *= minIndex;
	gd->mol_weight /= minIndex;
	/* calculate average weight per atom */
	gd->atom_weight = gd->mol_weight/nAtoms;

	mie_write_form(gd->elmAbun,chWord);
	fprintf( ioQQQ, "\n The chemical formula of the new grain molecule is: %s\n", chWord );
	fprintf( ioQQQ, " The abundance wrt H at maximum depletion of this molecule is: %.3e\n",
		 gd->abun );
	fprintf( ioQQQ, " The abundance wrt H at standard depletion of this molecule is: %.3e\n",
		 gd->abun*gd->depl );

	/* finally renormalize elmAbun back to abundance relative to hydrogen */
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		gd->elmAbun[nelem] *= gd->abun*gd->depl;
	}

	if( ( delta = (double*)MALLOC(sizeof(double)*(unsigned)sumAxes) ) == NULL )
		BadMalloc();
	if( ( frdelta = (double*)MALLOC(sizeof(double)*(unsigned)sumAxes) ) == NULL )
		BadMalloc();
	if( ( eps = (complex*)MALLOC(sizeof(complex)*(unsigned)sumAxes) ) == NULL )
		BadMalloc();

	l = 0;
	for( i=0; i < nMaterial; i++ )
	{
		for( k=0; k < gdArr[i].nAxes; k++ )
		{
			frdelta[l] = gdArr[i].wt[k]*frac[i];
			delta[l] = ( l == 0 ) ? frdelta[l] : delta[l-1] + frdelta[l];
			++l;
		}
	}
	ASSERT( l == sumAxes && fabs(delta[l-1]-1.) < 10.*DBL_EPSILON );

	/* allocate space for wavelength and optical constants arrays
	 * the memory is freed in mie_write_opc */
	if( ( gd->wavlen[0] = (double*)MALLOC(sizeof(double)*(unsigned)gd->ndata[0]) ) == NULL )
		BadMalloc();
	if( ( gd->n[0] = (complex*)MALLOC(sizeof(complex)*(unsigned)gd->ndata[0]) ) == NULL )
		BadMalloc();
	if( ( gd->nr1[0] = (double*)MALLOC(sizeof(double)*(unsigned)gd->ndata[0]) ) == NULL )
		BadMalloc();

	wavStep = log(wavHi/wavLo)/(double)(gd->ndata[0]-1);

	switch( EMTtype )
	{
	case FARAFONOV00:
		/* this implements the EMT described in
		 * >>refer Voshchinnikov N.V., Mathis J.S., 1999, ApJ, 526, 257
		 * based on the theory described in
		 * >>refer Farafonov V.G., 2000, Optics & Spectroscopy, 88, 441
		 *
		 * NB - note that Eq. 3 in Voshchinnikov & Mathis is incorrect! */
		for( j=0; j < gd->ndata[0]; j++ )
		{
			double nre,nim;
			complex a1,a2,a1c,a2c,a11,a12,a21,a22,ratio;

			gd->wavlen[0][j] = wavLo*exp((double)j*wavStep);

			init_eps(gd->wavlen[0][j],nMaterial,gdArr,eps);

			ratio = cmdiv(eps[0],eps[1]);

			a1 = ftocm((ratio.re+2.)/3.,ratio.im/3.);
			a2 = ftocm(-(ratio.re-1.)*delta[0],-ratio.im*delta[0]);

			for( l=1; l < sumAxes-1; l++ )
			{
				ratio = cmdiv(eps[l],eps[l+1]);

				a1c = a1;
				a2c = a2;
				a11 = ftocm((ratio.re+2.)/3.,ratio.im/3.);
				a12 = ftocm((2.-2.*ratio.re)/(9.*delta[l]),-2.*ratio.im/(9.*delta[l]));
				a21 = ftocm(-(ratio.re-1.)*delta[l],-ratio.im*delta[l]);
				a22 = ftocm((2.*ratio.re+1.)/3.,2.*ratio.im/3.);

				a1 = cmadd(cmmul(a11,a1c),cmmul(a12,a2c));
				a2 = cmadd(cmmul(a21,a1c),cmmul(a22,a2c));
			}

			a1c = a1;
			a2c = a2;
			a11 = ftocm(1.,0.);
			a12 = ftocm(1./3.,0.);
			a21 = eps[sumAxes-1];
			a22 = ftocm(-2./3.*eps[sumAxes-1].re,-2./3.*eps[sumAxes-1].im);

			a1 = cmadd(cmmul(a11,a1c),cmmul(a12,a2c));
			a2 = cmadd(cmmul(a21,a1c),cmmul(a22,a2c));

			ratio = cmdiv(a2,a1);
			dftori(&nre,&nim,ratio.re,ratio.im);

			gd->n[0][j] = ftocm(nre,nim);
			gd->nr1[0][j] = nre-1.;
		}
		break;
	case STOGNIENKO95:
	case BRUGGEMAN35:
		for( j=0; j < gd->ndata[0]; j++ )
		{
			const double EPS_TOLER = 10.*DBL_EPSILON;
			double nre,nim;
			complex eps0;

			gd->wavlen[0][j] = wavLo*exp((double)j*wavStep);

			init_eps(gd->wavlen[0][j],nMaterial,gdArr,eps);

			/* get initial estimate for effective dielectric function */
			if( j == 0 )
			{
				/* use simple average as first estimate */
				eps0 = ftocm(0.,0.);
				for( l=0; l < sumAxes; l++ )
				{
					eps0.re += frdelta[l]*eps[l].re;
					eps0.im += frdelta[l]*eps[l].im;
				}
			}
			else
			{
				/* use solution from previous wavelength as first estimate */
				eps0 = eps_eff;
			}

			if( EMTtype == STOGNIENKO95 )
				/* this implements the EMT described in
				 * >>refer Stognienko R., Henning Th., Ossenkopf V., 1995, A&A, 296, 797 */
				eps_eff = cnewton( Stognienko, frdelta, eps, sumAxes, eps0, EPS_TOLER );
			else if( EMTtype == BRUGGEMAN35 )
				/* this implements the classical Bruggeman rule
				 * >>refer Bruggeman D.A.G., 1935, Ann. Phys. (5th series), 24, 636 */
				eps_eff = cnewton( Bruggeman, frdelta, eps, sumAxes, eps0, EPS_TOLER );
			else
			{
				fprintf( ioQQQ, " Insanity in mie_read_mix\n" );
				ShowMe();
				puts( "[Stop in mie_read_mix]" );
				cdEXIT(EXIT_FAILURE);
			}

			dftori(&nre,&nim,eps_eff.re,eps_eff.im);

			gd->n[0][j] = ftocm(nre,nim);
			gd->nr1[0][j] = nre-1.;
		}
		break;
	default:
		fprintf( ioQQQ, " Insanity in mie_read_mix\n" );
		ShowMe();
		puts( "[Stop in mie_read_mix]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* clean up memory allocated by mie_read_rfi */
	for( i=0; i < nMaterial; i++ )
	{
		for( j=0; j < gdArr[i].nOpcCols; j++ ) 
		{
			if( gdArr[i].opcData[j] != NULL )
				free(gdArr[i].opcData[j]);
		}
		if( gdArr[i].opcAnu != NULL )
			free(gdArr[i].opcAnu);
		for( j=0; j < gdArr[i].nAxes; j++ ) 
		{
			if( gdArr[i].nr1[j] != NULL )
				free(gdArr[i].nr1[j]);
			if( gdArr[i].n[j] != NULL )
				free(gdArr[i].n[j]);
			if( gdArr[i].wavlen[j] != NULL )
				free(gdArr[i].wavlen[j]);
		}
	}
	/* these were locally allocated */
	free(eps);
	free(frdelta);
	free(delta);
	free(gdArr);
	free(frac2);
	free(frac);

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


/* helper routine for mie_read_mix, initializes the array of dielectric functions */
static void init_eps(double wavlen,
		     long nMaterial,
		     /*@in@*/ grain_data gdArr[],  /* gdArr[nMaterial] */
		     /*@out@*/ complex eps[])      /* eps[sumAxes] */
{
	long i,
	  k;

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

	long l = 0;
	for( i=0; i < nMaterial; i++ )
	{
		for( k=0; k < gdArr[i].nAxes; k++ )
		{
			int lgErr;
			long ind;
			double eps1,eps2,frc,nim,nre;
	  
			find_arr(wavlen,gdArr[i].wavlen[k],gdArr[i].ndata[k],&ind,&lgErr);
			ASSERT( !lgErr );
			frc = (wavlen-gdArr[i].wavlen[k][ind])/(gdArr[i].wavlen[k][ind+1]-gdArr[i].wavlen[k][ind]);
			ASSERT( frc > 0.-10.*DBL_EPSILON && frc < 1.+10.*DBL_EPSILON );
			nre = (1.-frc)*gdArr[i].n[k][ind].re + frc*gdArr[i].n[k][ind+1].re;
			ASSERT( nre > 0. );
			nim = (1.-frc)*gdArr[i].n[k][ind].im + frc*gdArr[i].n[k][ind+1].im;
			ASSERT( nim >= 0. );
			ritodf(nre,nim,&eps1,&eps2);
			eps[l++] = ftocm(eps1,eps2);
		}
	}

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


/* find complex root of fun using the Newton-Raphson algorithm */
static complex cnewton(void(*fun)(complex,double[],complex[],long,complex*,double*,double*),
		       /*@in@*/ double frdelta[], /* frdelta[sumAxes] */
		       /*@in@*/ complex eps[],    /* eps[sumAxes] */
		       long sumAxes,
		       complex x0,
		       double tol)
{
	int i;

	const int LOOP_MAX = 100;
	const double TINY = 1.e-12;

#	ifdef DEBUG_FUN
	fputs( "<+>cnewton()\n", debug_fp );
#	endif
	for( i=0; i < LOOP_MAX; i++ )
	{
		complex x1,y;
		double dudx,dudy,norm2;

		(*fun)(x0,frdelta,eps,sumAxes,&y,&dudx,&dudy);

		norm2 = POW2(dudx) + POW2(dudy);
		/* guard against norm2 == 0 */
		if( norm2 < TINY*(POW2(y.re) + POW2(y.im)) )
		{
			fprintf( ioQQQ, " cnewton - zero divide error\n" );
			ShowMe();
			puts( "[Stop in cnewton]" );
			cdEXIT(EXIT_FAILURE);
		}
		x1.re = x0.re + (y.im*dudy - y.re*dudx)/norm2;
		x1.im = x0.im - (y.im*dudx + y.re*dudy)/norm2;

		/* check for convergence */
		if( fabs(x0.re/x1.re-1.) + fabs(x0.im/x1.im-1.) < tol )
		{
#			ifdef DEBUG_FUN
			fputs( " <->cnewton()\n", debug_fp );
#			endif
			return x1;
		}
			
		x0 = x1;
	}

	fprintf( ioQQQ, " cnewton did not converge\n" );
	ShowMe();
	puts( "[Stop in cnewton]" );
	cdEXIT(EXIT_FAILURE);
}


/* this evaluates the function defined in Eq. 3 of
 * >>refer Stognienko R., Henning Th., Ossenkopf V., 1995, A&A, 296, 797
 * and its derivatives */
static void Stognienko(complex x,
		       /*@in@*/ double frdelta[], /* frdelta[sumAxes] */
		       /*@in@*/ complex eps[],    /* eps[sumAxes] */
		       long sumAxes,
		       /*@out@*/ complex *f,
		       /*@out@*/ double *dudx,
		       /*@out@*/ double *dudy)
{
	long i,
	  l;

	static const double L[4] = { 0., 1./2., 1., 1./3. };
	static const double fl[4] = { 5./9., 2./9., 2./9., 1. };

#	ifdef DEBUG_FUN
	fputs( "<+>Stognienko()\n", debug_fp );
#	endif
	*f = ftocm(0.,0.);
	*dudx = 0.;
	*dudy = 0.;
	for( l=0; l < sumAxes; l++ )
	{
		complex hlp = cmsub(eps[l],x);
		double h1 = eps[l].im*x.re - eps[l].re*x.im;

		for( i=0; i < 4; i++ )
		{
			double f1 = fl[i]*frdelta[l];
			double xx = ( i < 3 ) ? sin(PI*frdelta[l]) : cos(PI*frdelta[l]);
			complex f2 = ftocm(f1*xx*xx,0.);
			complex one = cmadd(x,cmmul(hlp,ftocm(L[i],0.)));
			complex two = cmmul(f2,cmdiv(hlp,one));
			double h2 = POW2(one.re) + POW2(one.im);
			*f = cmadd(*f,two);
			*dudx -= f2.re*(eps[l].re*h2 + h1*2.*one.im*(1.-L[i]))/POW2(h2);
			*dudy -= f2.re*(eps[l].im*h2 - h1*2.*one.re*(1.-L[i]))/POW2(h2);
		}
	}

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


/* this evaluates the classical Bruggeman rule and its derivatives
 * >>refer Bruggeman D.A.G., 1935, Ann. Phys. (5th series), 24, 636 */
static void Bruggeman(complex x,
		      /*@in@*/ double frdelta[], /* frdelta[sumAxes] */
		      /*@in@*/ complex eps[],    /* eps[sumAxes] */
		      long sumAxes,
		      /*@out@*/ complex *f,
		      /*@out@*/ double *dudx,
		      /*@out@*/ double *dudy)
{
	long l;

	static const double L = 1./3.;

#	ifdef DEBUG_FUN
	fputs( "<+>Bruggeman()\n", debug_fp );
#	endif
	*f = ftocm(0.,0.);
	*dudx = 0.;
	*dudy = 0.;
	for( l=0; l < sumAxes; l++ )
	{
		complex hlp = cmsub(eps[l],x);
		double h1 = eps[l].im*x.re - eps[l].re*x.im;
		complex f2 = ftocm(frdelta[l],0.);
		complex one = cmadd(x,cmmul(hlp,ftocm(L,0.)));
		complex two = cmmul(f2,cmdiv(hlp,one));
		double h2 = POW2(one.re) + POW2(one.im);
		*f = cmadd(*f,two);
		*dudx -= f2.re*(eps[l].re*h2 + h1*2.*one.im*(1.-L))/POW2(h2);
		*dudy -= f2.re*(eps[l].im*h2 - h1*2.*one.re*(1.-L))/POW2(h2);
	}

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


/* read in the file with optical constants and other relevant information */
static void mie_read_szd(/*@in@*/  const char *chFile,
			 /*@out@*/ sd_data *sd)
{
	int lgTryOverride = FALSE;
	long int dl,
	  i;
	double mul = 0.,
	  ref_neg = DBL_MAX,
	  ref_pos = DBL_MAX,
	  step_neg = DBL_MAX,
	  step_pos = DBL_MAX;
	char chLine[FILENAME_PATH_LENGTH_2],
	  chWord[WORDLEN];
	FILE *io2;

	/* these constants are used to get initial estimates for the cutoffs (lim)
	 * in the SD_EXP_CUTOFFx and SD_xxx_NORMAL cases, they are iterated by
	 * search_limit such that
	 * lim^4 * dN/da(lim) == FRAC_CUTOFF * ref^4 * dN/da(ref)
	 * where ref as an appropriate reference point for each of the cases */
	const double FRAC_CUTOFF = 1.e-4;
	const double MUL_CO1 = -log(FRAC_CUTOFF);
	const double MUL_CO2 = sqrt(MUL_CO1);
	const double MUL_CO3 = pow(MUL_CO1,1./3.);
	const double MUL_LND = sqrt(-2.*log(FRAC_CUTOFF));
	const double MUL_NRM = MUL_LND;

#	ifdef DEBUG_FUN
	fputs( "<+>mie_read_szd()\n", debug_fp );
#	endif
	if( (io2 = fopen(chFile,"r")) == NULL )
	{
		fprintf( ioQQQ, " Could not open %s for reading\n",chFile);
		puts( "[Stop in mie_read_szd]" );
		cdEXIT(EXIT_FAILURE);
	}

	dl = 0; /* line counter for input file */

	/* first read magic number */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_long(chFile,chLine,&sd->magic,TRUE,dl);
	if( sd->magic != MAGIC_SZD ) 
	{
		fprintf( ioQQQ, " Size distribution file %s has obsolete magic number\n",chFile );
		fprintf( ioQQQ, " I found %ld, but expected %ld on line #%ld\n",sd->magic,MAGIC_SZD,dl );
		fprintf( ioQQQ, " Please replace this file with an up to date version\n" );
		puts( "[Stop in mie_read_szd]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* size distribution case */
	mie_next_data(chFile,io2,chLine,&dl);
	mie_read_word(chLine,chWord,WORDLEN,TRUE);

	if( lgMatch( chWord, "SSIZ" ) )
	{
		sd->sdCase = SD_SINGLE_SIZE;
	}
	else if( lgMatch( chWord, "POWE" ) )
	{
		sd->sdCase = SD_POWERLAW;
	}
	else if( lgMatch( chWord, "EXP1" ) )
	{
		sd->sdCase = SD_EXP_CUTOFF1;
		sd->a[ipAlpha] = 1.;
		mul = MUL_CO1;
	}
	else if( lgMatch( chWord, "EXP2" ) )
	{
		sd->sdCase = SD_EXP_CUTOFF2;
		sd->a[ipAlpha] = 2.;
		mul = MUL_CO2;
	}
	else if( lgMatch( chWord, "EXP3" ) )
	{
		sd->sdCase = SD_EXP_CUTOFF3;
		sd->a[ipAlpha] = 3.;
		mul = MUL_CO3;
	}
	else if( lgMatch( chWord, "LOGN" ) )
	{
		sd->sdCase = SD_LOG_NORMAL;
		mul = MUL_LND;
	}
	/* this one must come after LOGNORMAL */
	else if( lgMatch( chWord, "NORM" ) )
	{
		sd->sdCase = SD_LIN_NORMAL;
		mul = MUL_NRM;
	}
	else if( lgMatch( chWord, "TABL" ) )
	{
		sd->sdCase = SD_TABLE;
	}
	else
	{
		sd->sdCase = SD_ILLEGAL;
	}

	switch( sd->sdCase ) 
	{
	case SD_SINGLE_SIZE:
		/* single sized grain */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_double(chFile,chLine,&sd->a[ipSize],TRUE,dl);
		if( sd->a[ipSize] < SMALLEST_GRAIN || sd->a[ipSize] > LARGEST_GRAIN ) 
		{
			fprintf( ioQQQ, " Illegal value for grain size\n" );
			fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
				 SMALLEST_GRAIN, LARGEST_GRAIN );
			fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
			puts( "[Stop in mie_read_szd]" );
			cdEXIT(EXIT_FAILURE);
		}
		break;
	case SD_POWERLAW:
		/* simple power law distribution, first get lower limit */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_double(chFile,chLine,&sd->a[ipBLo],TRUE,dl);
		if( sd->a[ipBLo] < SMALLEST_GRAIN || sd->a[ipBLo] > LARGEST_GRAIN ) 
		{
			fprintf( ioQQQ, " Illegal value for grain size (lower limit)\n" );
			fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
				 SMALLEST_GRAIN, LARGEST_GRAIN );
			fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
			puts( "[Stop in mie_read_szd]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* upper limit */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_double(chFile,chLine,&sd->a[ipBHi],TRUE,dl);
		if( sd->a[ipBHi] < SMALLEST_GRAIN || sd->a[ipBHi] > LARGEST_GRAIN ||
		    sd->a[ipBHi] <= sd->a[ipBLo] ) 
		{
			fprintf( ioQQQ, " Illegal value for grain size (upper limit)\n" );
			fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
				 SMALLEST_GRAIN, LARGEST_GRAIN );
			fprintf( ioQQQ, " and upper limit should be greater than lower limit\n" );
			fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
			puts( "[Stop in mie_read_szd]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* slope */
		mie_next_data(chFile,io2,chLine,&dl);
		if( sscanf( chLine, "%lf", &sd->a[ipExp] ) != 1 ) 
		{
			fprintf( ioQQQ, " Syntax error in %s\n",chFile);
			fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
			puts( "[Stop in mie_read_szd]" );
			cdEXIT(EXIT_FAILURE);
		}

		sd->a[ipBeta] = 0.;
		sd->a[ipSLo] = 0.;
		sd->a[ipSHi] = 0.;

		sd->lim[ipBLo] = sd->a[ipBLo];
		sd->lim[ipBHi] = sd->a[ipBHi];
		break;
	case SD_EXP_CUTOFF1:
	case SD_EXP_CUTOFF2:
	case SD_EXP_CUTOFF3:
		/* powerlaw with first/second/third order exponential cutoff, inspired by
		 * Greenberg (1978), Cosmic Dust, ed. J.A.M. McDonnell, Wiley, p. 187 */
		/* "lower limit", below this the exponential cutoff sets in */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_double(chFile,chLine,&sd->a[ipBLo],TRUE,dl);

		/* "upper" limit, above this the exponential cutoff sets in */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_double(chFile,chLine,&sd->a[ipBHi],TRUE,dl);

		/* exponent for power law */
		mie_next_data(chFile,io2,chLine,&dl);
		if( sscanf( chLine, "%lf", &sd->a[ipExp] ) != 1 ) 
		{
			fprintf( ioQQQ, " Syntax error in %s\n",chFile);
			fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
			puts( "[Stop in mie_read_szd]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* beta parameter, for extra curvature in the powerlaw region */
		mie_next_data(chFile,io2,chLine,&dl);
		if( sscanf( chLine, "%lf", &sd->a[ipBeta] ) != 1 ) 
		{
			fprintf( ioQQQ, " Syntax error in %s\n",chFile);
			fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
			puts( "[Stop in mie_read_szd]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* scale size for lower exponential cutoff, zero indicates normal cutoff */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_double(chFile,chLine,&sd->a[ipSLo],FALSE,dl);

		/* scale size for upper exponential cutoff, zero indicates normal cutoff */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_double(chFile,chLine,&sd->a[ipSHi],FALSE,dl);

		ref_neg = sd->a[ipBLo];
		step_neg = -mul*sd->a[ipSLo];
		ref_pos = sd->a[ipBHi];
		step_pos = mul*sd->a[ipSHi];
		lgTryOverride = TRUE;
		break;
	case SD_LOG_NORMAL:
		/* log-normal distribution, first get center of gaussian */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_double(chFile,chLine,&sd->a[ipGCen],TRUE,dl);

		/* 1-sigma width */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_double(chFile,chLine,&sd->a[ipGSig],TRUE,dl);

		/* ref_pos, ref_neg is the grain radius at which a^4*dN/da peaks */
		ref_neg = ref_pos = sd->a[ipGCen]*exp(3.*POW2(sd->a[ipGSig]));
		step_neg = sd->a[ipGCen]*(exp(-mul*sd->a[ipGSig]) - 1.);
		step_pos = sd->a[ipGCen]*(exp(mul*sd->a[ipGSig]) - 1.);
		lgTryOverride = TRUE;
		break;
	case SD_LIN_NORMAL:
		/* normal gaussian distribution, first get center of gaussian */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_double(chFile,chLine,&sd->a[ipGCen],TRUE,dl);

		/* 1-sigma width */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_double(chFile,chLine,&sd->a[ipGSig],TRUE,dl);

		/* ref_pos, ref_neg is the grain radius at which a^4*dN/da peaks */
		ref_neg = ref_pos = (sd->a[ipGCen] + sqrt(POW2(sd->a[ipGCen]) + 12.*POW2(sd->a[ipGSig])))/2.;
		step_neg = -mul*sd->a[ipGSig];
		step_pos = mul*sd->a[ipGSig];
		lgTryOverride = TRUE;
		break;
	case SD_TABLE:
		/* user-supplied table of a^4*dN/da vs. a, first get lower limit on a */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_double(chFile,chLine,&sd->a[ipBLo],TRUE,dl);
		if( sd->a[ipBLo] < SMALLEST_GRAIN || sd->a[ipBLo] > LARGEST_GRAIN ) 
		{
			fprintf( ioQQQ, " Illegal value for grain size (lower limit)\n" );
			fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
				 SMALLEST_GRAIN, LARGEST_GRAIN );
			fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
			puts( "[Stop in mie_read_szd]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* upper limit */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_double(chFile,chLine,&sd->a[ipBHi],TRUE,dl);
		if( sd->a[ipBHi] < SMALLEST_GRAIN || sd->a[ipBHi] > LARGEST_GRAIN ||
		    sd->a[ipBHi] <= sd->a[ipBLo] ) 
		{
			fprintf( ioQQQ, " Illegal value for grain size (upper limit)\n" );
			fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
				 SMALLEST_GRAIN, LARGEST_GRAIN );
			fprintf( ioQQQ, " and upper limit should be greater than lower limit\n" );
			fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
			puts( "[Stop in mie_read_szd]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* number of user supplied points */
		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_long(chFile,chLine,&sd->npts,TRUE,dl);
		if( sd->npts < 2 )
		{
			fprintf( ioQQQ, " Illegal value for no. of points in table\n" );
			fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
			puts( "[Stop in mie_read_szd]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* allocate space for the table */
		sd->ln_a = (double *)MALLOC(sizeof(double)*(unsigned)sd->npts);
		sd->ln_a4dNda = (double *)MALLOC(sizeof(double)*(unsigned)sd->npts);
		if( sd->ln_a == NULL || sd->ln_a4dNda == NULL )
		{
			fprintf( ioQQQ, " Could not MALLOC arrays\n" );
			puts( "[Stop in mie_read_szd]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* and read the table */
		for( i=0; i < sd->npts; ++i )
		{
			double help1, help2;

			mie_next_data(chFile,io2,chLine,&dl);
			/* read data pair: a (micron), a^4*dN/da (arbitrary units) */
			if( sscanf( chLine, "%le %le", &help1, &help2 ) != 2 ) 
			{
				fprintf( ioQQQ, " Syntax error in %s\n",chFile);
				fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
				puts( "[Stop in mie_read_szd]" );
				cdEXIT(EXIT_FAILURE);
			}

			if( help1 <= 0. || help2 <= 0. )
			{
				fprintf( ioQQQ, " Reading table failed on line #%ld of %s\n",dl,chFile );
				fprintf( ioQQQ, " Illegal data value %.6e or %.6e\n", help1, help2 );
				puts( "[Stop in mie_read_szd]" );
				cdEXIT(EXIT_FAILURE);
			}

			sd->ln_a[i] = log(help1);
			sd->ln_a4dNda[i] = log(help2);

			if( i > 0 && sd->ln_a[i] <= sd->ln_a[i-1] )
			{
				fprintf( ioQQQ, " Reading table failed on line #%ld of %s\n",dl,chFile );
				fprintf( ioQQQ, " Grain radii should be monotonically increasing\n" );
				puts( "[Stop in mie_read_szd]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		sd->lim[ipBLo] = sd->a[ipBLo];
		sd->lim[ipBHi] = sd->a[ipBHi];
		break;
	case SD_ILLEGAL:
	default:
		fprintf( ioQQQ, " unimplemented case for grain size distribution in file %s\n" , chFile );
		fprintf( ioQQQ, " Line #%ld: value %s\n",dl,chWord);
		puts( "[Stop in mie_read_szd]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* >>chng 01 feb 12, use a^4*dN/da instead of dN/da to determine limits,
	 * this assures that upper limit gives negligible mass fraction, PvH */
	/* in all cases where search_limit is used to determine lim[ipBLo] and lim[ipBHi],
	 * the user can override these values in the last two lines of the size distribution
	 * file. these inputs are mandatory, and should be given in the sequence lower
	 * limit, upper limit. a value <= 0 indicates that search_limit should be used. */
	if( lgTryOverride )
	{
		double help;

		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_double(chFile,chLine,&help,FALSE,dl);
		sd->lim[ipBLo] = ( help <= 0. ) ? search_limit(ref_neg,step_neg,FRAC_CUTOFF,*sd) : help;

		mie_next_data(chFile,io2,chLine,&dl);
		mie_read_double(chFile,chLine,&help,FALSE,dl);
		sd->lim[ipBHi] = ( help <= 0. ) ? search_limit(ref_pos,step_pos,FRAC_CUTOFF,*sd) : help;

		if( sd->lim[ipBLo] < SMALLEST_GRAIN || sd->lim[ipBHi] > LARGEST_GRAIN ||
		    sd->lim[ipBHi] <= sd->lim[ipBLo] ) 
		{
			fprintf( ioQQQ, " Illegal size limits: lower %.5f and/or upper %.5f\n",
				 sd->lim[ipBLo], sd->lim[ipBHi] );
			fprintf( ioQQQ, " Grain sizes should be between %.5f and %.0f micron\n",
				 SMALLEST_GRAIN, LARGEST_GRAIN );
			fprintf( ioQQQ, " and upper limit should be greater than lower limit\n" );
			fprintf( ioQQQ, " Please alter the size distribution file\n" );
			puts( "[Stop in mie_read_szd]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	fclose(io2);

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


static void mie_read_long(/*@in@*/  const char *chFile,
			  /*@in@*/  const char chLine[], /* chLine[FILENAME_PATH_LENGTH_2] */
			  /*@out@*/ long int *data,
			  int lgZeroIllegal,
			  long int dl)
{
#	ifdef DEBUG_FUN
	fputs( "<+>mie_read_long()\n", debug_fp );
#	endif
	if( sscanf( chLine, "%ld", data ) != 1 )
	{
		fprintf( ioQQQ, " Syntax error in %s\n",chFile);
		fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
		puts( "[Stop in mie_read_long]" );
		cdEXIT(EXIT_FAILURE);
	}
	if( *data < 0 || (*data == 0 && lgZeroIllegal) )
	{
		fprintf( ioQQQ, " Illegal data value in %s\n",chFile);
		fprintf( ioQQQ, " Line #%ld: %ld\n",dl,*data);
		puts( "[Stop in mie_read_long]" );
		cdEXIT(EXIT_FAILURE);
	}
#	ifdef DEBUG_FUN
	fputs( " <->mie_read_long()\n", debug_fp );
#	endif
	return;
}


static void mie_read_float(/*@in@*/  const char *chFile,
			   /*@in@*/  const char chLine[], /* chLine[FILENAME_PATH_LENGTH_2] */
			   /*@out@*/ float *data,
			   int lgZeroIllegal,
			   long int dl)
{
#	ifdef DEBUG_FUN
	fputs( "<+>mie_read_float()\n", debug_fp );
#	endif
	if( sscanf( chLine, "%f", data ) != 1 )
	{
		fprintf( ioQQQ, " Syntax error in %s\n",chFile);
		fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
		puts( "[Stop in mie_read_float]" );
		cdEXIT(EXIT_FAILURE);
	}
	if( *data < 0. || (*data == 0. && lgZeroIllegal) )
	{
		fprintf( ioQQQ, " Illegal data value in %s\n",chFile);
		fprintf( ioQQQ, " Line #%ld: %14.6e\n",dl,*data);
		puts( "[Stop in mie_read_float]" );
		cdEXIT(EXIT_FAILURE);
	}
#	ifdef DEBUG_FUN
	fputs( " <->mie_read_float()\n", debug_fp );
#	endif
	return;
}


static void mie_read_double(/*@in@*/  const char *chFile,
			    /*@in@*/  const char chLine[], /* chLine[FILENAME_PATH_LENGTH_2] */
			    /*@out@*/ double *data,
			    int lgZeroIllegal,
			    long int dl)
{
#	ifdef DEBUG_FUN
	fputs( "<+>mie_read_double()\n", debug_fp );
#	endif
	if( sscanf( chLine, "%lf", data ) != 1 )
	{
		fprintf( ioQQQ, " Syntax error in %s\n",chFile);
		fprintf( ioQQQ, " Line #%ld: %s\n",dl,chLine);
		puts( "[Stop in mie_read_double]" );
		cdEXIT(EXIT_FAILURE);
	}
	if( *data < 0. || (*data == 0. && lgZeroIllegal) )
	{
		fprintf( ioQQQ, " Illegal data value in %s\n",chFile);
		fprintf( ioQQQ, " Line #%ld: %14.6e\n",dl,*data);
		puts( "[Stop in mie_read_double]" );
		cdEXIT(EXIT_FAILURE);
	}
#	ifdef DEBUG_FUN
	fputs( " <->mie_read_double()\n", debug_fp );
#	endif
	return;
}


static void mie_read_form(/*@in@*/  const char chWord[], /* chWord[FILENAME_PATH_LENGTH_2] */
			  /*@out@*/ double elmAbun[],    /* elmAbun[LIMELM] */
			  /*@out@*/ double *no_atoms,
			  /*@out@*/ double *mol_weight)
{
	long int nelem,
	  len;
	const char *str;
	char chElmName[3];
	double frac;

#	ifdef DEBUG_FUN
	fputs( "<+>mie_read_form()\n", debug_fp );
#	endif
	*no_atoms = 0.;
	*mol_weight = 0.;
	for( nelem=0; nelem < LIMELM; nelem++ ) 
	{
		frac = 0.;
		strcpy(chElmName,elementnames.chElementSym[nelem]);
		if( chElmName[1] == ' ' )
			chElmName[1] = '\0';
		str = strstr(chWord,chElmName);
		if( str != NULL ) 
		{
			len = (long)strlen(chElmName);
			/* prevent spurious match, e.g. F matches Fe */
			if( !islower((unsigned char)str[len]) ) 
			{
				if( isdigit((unsigned char)str[len]) ) 
				{
					sscanf(str+len,"%lf",&frac);
				}
				else 
				{
					frac = 1.;
				}
			}
		}
		elmAbun[nelem] = frac;
		/* >>chng 02 apr 22, don't count hydrogen in PAH's, PvH */
		if( nelem != ipHYDROGEN )
			*no_atoms += frac;
		*mol_weight += frac*dense.AtomicWeight[nelem];
	}
	/* prevent division by zero when no chemical formula was supplied */
	if( *no_atoms == 0. ) 
		*no_atoms = 1.;

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


static void mie_write_form(/*@in@*/  const double elmAbun[], /* elmAbun[LIMELM] */
			   /*@out@*/ char chWord[])          /* chWord[FILENAME_PATH_LENGTH_2] */
{
	long int len,
	  nelem;

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

	len=0;
	for( nelem=0; nelem < LIMELM; nelem++ ) 
	{
		if( elmAbun[nelem] > 0. )
		{
			char chElmName[3];
			long index100 = NINT(100.*elmAbun[nelem]);

			strcpy(chElmName,elementnames.chElementSym[nelem]);
			if( chElmName[1] == ' ' )
				chElmName[1] = '\0';

			if( index100 == 100 )
				sprintf( &chWord[len], "%s", chElmName );
			else if( index100%100 == 0 )
				sprintf( &chWord[len], "%s%ld", chElmName, index100/100 );
			else
			{
				double xIndex = (double)index100/100.;
				sprintf( &chWord[len], "%s%.2f", chElmName, xIndex );
			}
			len = (long)strlen( chWord );
			ASSERT( len < FILENAME_PATH_LENGTH_2-25 );
		}
	}

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


static void mie_read_word(/*@in@*/  const char chLine[], /* chLine[FILENAME_PATH_LENGTH_2] */
			  /*@out@*/ char chWord[],       /* chWord[n] */
			  long n,
			  int toUpper)
{
	long ip = 0,
	  op = 0;

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

	/* skip leading spaces or double quotes */
	while( chLine[ip] == ' ' || chLine[ip] == '"' )
		ip++;
	/* now copy string until we hit next space or double quote */
	while( op < n-1 && chLine[ip] != ' ' && chLine[ip] != '"' )
		if( toUpper )
			chWord[op++] = (char)toupper(chLine[ip++]);
		else
			chWord[op++] = chLine[ip++];
	/* the output string is always zero terminated */
	chWord[op] = '\0';

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

/*=====================================================================*/
static void mie_next_data(/*@in@*/  const char *chFile,
			  /*@in@*/  FILE *io,
			  /*@out@*/ char chLine[], /* chLine[FILENAME_PATH_LENGTH_2] */
			  /*@in@*/  long int *dl)
{
	char *str;

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

	/* lines starting with a pound sign are considered comments and are skipped,
	 * lines not starting with a pound sign are considered to contain useful data.
	 * however, comments may still be appended to the line and will be erased. */

	chLine[0] = '#';
	while( chLine[0] == '#' )
	{
		mie_next_line(chFile,io,chLine,dl);
	}

	/* erase comment part of the line */
	str = strstr(chLine,"#");
	if( str != NULL )
		*str = '\0';

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


/*=====================================================================*/
static void mie_next_line(/*@in@*/  const char *chFile,
			  /*@in@*/  FILE *io,
			  /*@out@*/ char chLine[], /* chLine[FILENAME_PATH_LENGTH_2] */
			  /*@in@*/  long int *dl)
{
#	ifdef DEBUG_FUN
	fputs( "<+>mie_next_line()\n", debug_fp );
#	endif

	if( fgets( chLine, FILENAME_PATH_LENGTH_2, io ) == NULL ) 
	{
		fprintf( ioQQQ, " Could not read from %s\n",chFile);
		if( feof(io) )
			fprintf( ioQQQ, " EOF reached\n");
		puts( "[Stop in mie_next_line]" );
		cdEXIT(EXIT_FAILURE);
	}
	(*dl)++;

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

/* make next two gauss.c ??? */

/*=====================================================================*/
/* set up Gaussian quadrature for arbitrary interval */
void gauss_init(long int nn,
	    double xbot,
	    double xtop,
	    double x[],  /* x[nn]  */
	    double a[],  /* a[nn]  */
	    double rr[], /* rr[nn] */
	    double ww[]) /* ww[nn] */
{
	long int i;
	double bma,
	  bpa;

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

	bpa = (xtop+xbot)/2.;
	bma = (xtop-xbot)/2.;

	for( i=0; i < nn; i++ ) 
	{
		rr[i] = bpa + bma*x[nn-1-i];
		ww[i] = bma*a[i];
	}

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


/*=====================================================================*/
/* set up abscissas and weights for Gauss-Legendre intergration of arbitrary even order */
void gauss_legendre(long int nn,
	     double x[], /* x[nn] */
	     double a[]) /* a[nn] */
{
	long int i,
	  iter,
	  j;
	double *c,
	  cc,
	  csa,
	  d,
	  dp1,
	  dpn = 0.,
	  dq,
	  fj,
	  fn,
	  pn,
	  pn1 = 0.,
	  q,
	  xt = 0.;

	const double SAFETY = 5.;


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

	if( nn%2 == 1 ) 
	{
		fprintf( ioQQQ, " Illegal number of abcissas\n" );
		puts( "[Stop in gauss_legendre]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( ( c = (double *)MALLOC(sizeof(double)*(unsigned)nn) ) == NULL )
		BadMalloc();

	fn = (double)nn;
	csa = 0.;
	cc = 2.;
	for( j=1; j < nn; j++ ) 
	{
		fj = (double)j;
		/* >>chng 01 apr 10, prevent underflows in cc, pn, pn1, dpn and dp1 for large nn
		 * renormalize c[j] -> 4*c[j],  cc -> 4^(nn-1)*cc,  hence cc = O(1), etc...
		 * Old code: c[j] = POW2(fj)/(4.*(fj-0.5)*(fj+0.5)); */
		c[j] = POW2(fj)/((fj-0.5)*(fj+0.5));
		cc *= c[j];
	}

	for( i=0; i < nn/2; i++ ) 
	{
		switch( i ) 
		{
		case 0:
			xt = 1. - 2.78/(4. + POW2(fn));
			break;
		case 1:
			xt = xt - 4.1*(1. + 0.06*(1. - 8./fn))*(1. - xt);
			break;
		case 2:
			xt = xt - 1.67*(1. + 0.22*(1. - 8./fn))*(x[0] - xt);
			break;
		default:
			xt = 3.*(x[i-1] - x[i-2]) + x[i-3];
		}
		d = 1.;
		for( iter=1; (iter < 20) && (fabs(d) > DBL_EPSILON); iter++ ) 
		{
			/* >>chng 01 apr 10, renormalize pn -> 2^(nn-1)*pn, dpn -> 2^(nn-1)*dpn
			 * pn1 -> 2^(nn-2)*pn1, dp1 -> 2^(nn-2)*dp1
			 * Old code: pn1 = 1.; */
			pn1 = 0.5;
			pn = xt;
			dp1 = 0.;
			dpn = 1.;
			for( j=1; j < nn; j++ )
			{
				/* >>chng 01 apr 10, renormalize pn -> 2^(nn-1)*pn, dpn -> 2^(nn-1)*dpn
				 * Old code: q = xt*pn - c[j]*pn1;  dq = xt*dpn - c[j]*dp1 + pn; */
				q = 2.*xt*pn - c[j]*pn1;
				dq = 2.*xt*dpn - c[j]*dp1 + 2.*pn;
				pn1 = pn;
				pn = q;
				dp1 = dpn;
				dpn = dq;
			}
			d = pn/dpn;
			xt -= d;
		}
		x[i] = xt;
		x[nn-1-i] = -xt;
		/* >>chng 01 apr 10, renormalize dpn -> 2^(nn-1)*dpn, pn1 -> 2^(nn-2)*pn1
		 * Old code: a[i] = cc/(dpn*pn1); */
		a[i] = cc/(dpn*2.*pn1);
		a[nn-1-i] = a[i];
		csa += a[i];
	}

	/* this routine has been tested for every even nn between 2 and 4096
	 * it passed the test for each of those cases with SAFETY < 3.11 */
	if( fabs(1.-csa) > SAFETY*fn*DBL_EPSILON ) 
	{
		fprintf( ioQQQ, " gauss_legendre failed to converge: delta = %.4e\n",fabs(1.-csa) );
		puts( "[Stop in gauss_legendre]" );
		cdEXIT(EXIT_FAILURE);
	}
	free(c);

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


/* find index ind such that min(xa[ind],xa[ind+1]) <= x <= max(xa[ind],xa[ind+1]).
 * xa is assumed to be strictly monotically increasing or decreasing.
 * if x is outside the range spanned by xa, lgOutOfBounds is raised and ind is set to -1
 * n is the number of elements in xa. */
void find_arr(double x,
	      double xa[],
	      long int n,
	      /*@out@*/ long int *ind,
	      /*@out@*/ int *lgOutOfBounds)
{
	long int i1,
		i2,
		i3,
		sgn,
		sgn2;

#	ifdef DEBUG_FUN
	fputs( "<+>find_arr()\n", debug_fp );
#	endif
	/* this routine works for strictly monotically increasing
	 * and decreasing arrays, sgn indicates which case it is */
	if( n < 2 ) 
	{
		fprintf( ioQQQ, " Invalid array\n");
		puts( "[Stop in find_arr]" );
		cdEXIT(EXIT_FAILURE);
	}

	i1 = 0;
	i3 = n-1;
	sgn = SIGN3(xa[i3]-xa[i1]);
	if( sgn == 0 ) 
	{
		fprintf( ioQQQ, " Ill-ordered array\n");
		puts( "[Stop in find_arr]" );
		cdEXIT(EXIT_FAILURE);
	}

	*lgOutOfBounds = x < MIN2(xa[0],xa[n-1]) || x > MAX2(xa[0],xa[n-1]);
	if( *lgOutOfBounds ) 
	{
		*ind = -1;
#		ifdef DEBUG_FUN
		fputs( " <->find_arr()\n", debug_fp );
#		endif
		return;
	}

	i2 = (n-1)/2;
	while( (i3-i1) > 1 ) 
	{
		sgn2 = SIGN3(x-xa[i2]);
		if( sgn2 != 0 )
		{
			if( sgn == sgn2 ) 
			{
				i1 = i2;
			}
			else 
			{
				i3 = i2;
			}
			i2 = (i1+i3)/2;
		}
		else 
		{
			*ind = i2;
#			ifdef DEBUG_FUN
			fputs( " <->find_arr()\n", debug_fp );
#			endif
			return;
		}
	}
	*ind = i1;
#	ifdef DEBUG_FUN
	fputs( " <->find_arr()\n", debug_fp );
#	endif
	return;
}


/* Oct 1988 for UV - X-ray extinction, including anomalous diffraction check
 *     this version reads in real and imaginary parts of the refractive
 *     index, with imaginary part positive (nridf = 3) or nr-1 (nridf = 2) or
 *     real and imaginary parts of the dielectric function (nridf = 1)
 * Dec 1988: added qback; approximation for small x;
 * qphase, better convergence checking
 *
 * in anomalous diffraction: qext and qabs calculated - qscatt by subtraction
 * in rayleigh-gans:         qscatt and qabs calculated
 * in mie:                   qext and qscatt calculated
 * */

/* sinpar.f
 * consistency checks updated july 1999
 * t1 updated mildly 19 oct 1992
 * utility for mieuvx.f and mieuvxsd.f */
#define	NMXLIM	16000

static void sinpar(double nre,
		   double nim,
		   double x,
		   /*@out@*/ double *qext,
		   /*@out@*/ double *qphase,
		   /*@out@*/ double *qscat,
		   /*@out@*/ double *ctbrqs,
		   /*@out@*/ double *qback,
		   /*@out@*/ long int *iflag)
{
	long int n,
	  nmx1,
	  nmx2,
	  nn,
	  nsqbk;
	double ectb,
	  eqext,
	  eqpha,
	  eqscat,
	  error=0.,
	  error1=0.,
	  rx,
	  t1,
	  t2,
	  t3,
	  t4,
	  t5,
	  tx,
	  x3,
	  x5=0.,
	  xcut,
	  xrd;
	/* >>chng 01 sep 11, replace array allocation on stack with
	 * MALLOC to avoid bug in gcc 3.0 on Sparc platforms, PvH */
/*  	complex a[NMXLIM], */
	complex *a;
	complex cdum1,
	  cdum2,
	  ci,
	  eqb,
	  nc,
	  nc2,
	  nc212,
	  qbck,
	  rrf,
	  rrfx,
	  sman,
	  sman1,
	  smbn,
	  smbn1,
	  tc1,
	  tc2,
	  wn,
	  wn1,
	  wn2;

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

	if( ( a = (complex*)MALLOC(sizeof(complex)*(unsigned)NMXLIM) ) == NULL )
		BadMalloc();

	*iflag = 0;
	ci = ftocm(0.,1.);
	nc = ftocm(nre,-nim);
	nc2 = cmmul(nc,nc);
	rrf = cmdiv(ftocm(1.,0.),nc);
	rx = 1.0/x;
	rrfx = cmmul(rrf,ftocm(rx,0.));

	/*  t1 is the number of terms nmx2 that will be needed to obtain convergence
	 *  try to minimize this, because the a(n) downwards recursion has to
	 *  start at nmx1 larger than this
	 *
	 * major loop series is summed to nmx2, or less when converged
	 * nmx1 is used for a(n) only, n up to nmx2.
	 * must start evaluation sufficiently above nmx2 that a(nmx1)=(0.,0.)
	 * is a good approximation
	 *
	 *
	 *orig with slight modification for extreme UV and X-ray, n near 1., large x
	 *orig      t1=x*dmax1( 1.1d0,dsqrt(nr*nr+ni*ni) )*
	 *orig     1(1.d0+0.02d0*dmax1(dexp(-x/100.d0)*x/10.d0,dlog10(x)))
	 *
	 * rules like those of wiscombe 1980 are slightly more efficient */
	xrd = exp(log(x)/3.);
	/* the final number in t1 was 1., 2. for large x, and 3. is needed sometimes
	 * see also idnint use below */
	t1 = x + 4.*xrd + 3.;
	/*      t1=t1+0.05d0*xrd
	 * was 0., then 1., then 2., now 3. for intermediate x
	 * 19 oct 1992 */
	if( !(x <= 8. || x >= 4200.) )
		t1 += 0.05*xrd + 3.;
	t1 *= 1.01;

	/* the original rule of dave for starting the downwards recursion was
	 * to start at 1.1*|mx| + 1, i.e. with the original version of t1
	 *orig      nmx1=1.10d0*t1
	 *
	 * try a simpler, less costly one, as in bohren and huffman, p 478
	 * this is the form for use with wiscombe rules for t1
	 * tests: it produces the same results as the more costly version
	 * */
	t4 = x*sqrt(nre*nre+nim*nim);
	nmx1 = NINT(MAX2(t1,t4)) + 15;

	if( nmx1 < NMXLIM ) 
	{
		nmx2 = NINT(t1);
		/*orig      if ( nmx1  .gt. 150 ) go to 22
		 *orig      nmx1 = 150
		 *orig      nmx2 = 135
		 *
		 * try a more efficient scheme */
		if( nmx2 <= 4 ) 
		{
			nmx2 = 4;
			nmx1 = NINT(MAX2(4.,t4)) + 15;
		}

		/* downwards recursion for logarithmic derivative */
		a[nmx1] = ftocm(0.,0.);

		/* note that with the method of lentz 1976 (appl opt 15, 668), it would be
		 * possible to find a(nmx2) directly, and start the downwards recursion there
		 * however, there is not much in it with above form for nmx1 which uses just */
		for( n=0; n < nmx1; n++ )
		{
			nn = nmx1 - n;
			a[nn-1] = cmsub(cmmul(ftocm((double)(nn+1),0.),rrfx),cmdiv(ftocm(1.,0.),
			  (cmadd(cmmul(ftocm((double)(nn+1),0.),rrfx),a[nn]))));
		}

		t1 = cos(x);
		t2 = sin(x);
		wn2 = ftocm(t1,-t2);
		wn1 = ftocm(t2,t1);
		wn = cmsub(cmmul(ftocm(rx,0.),wn1),wn2);
		tc1 = cmadd(cmmul(a[0],rrf),ftocm(rx,0.));
		tc2 = cmadd(cmmul(a[0],nc),ftocm(rx,0.));
		sman = cmdiv((cmsub(cmmul(tc1,ftocm(wn.re,0.)),ftocm(wn1.re,0.))),
		  (cmsub(cmmul(tc1,wn),wn1)));
		smbn = cmdiv((cmsub(cmmul(tc2,ftocm(wn.re,0.)),ftocm(wn1.re,0.))),
		  (cmsub(cmmul(tc2,wn),wn1)));

		/* small x; above calculations subject to rounding errors
		 * see bohren and huffman p 131
		 * wiscombe 1980 appl opt 19, 1505 gives alternative formulation */
		xcut = 3.e-04;
		if( x < xcut ) 
		{
			nc212 = cmdiv((cmsub(nc2,ftocm(1.,0.))),(cmadd(nc2,ftocm(2.,0.))));
			x3 = POW3(x);
			x5 = x3*POW2(x);
			/* note change sign convention for m = n - ik here */
			sman = cmadd(cmmul(cmmul(cmmul(cmmul(ci,ftocm(2.,0.)),ftocm(x3,0.)),
			  nc212),(cmadd(ftocm(1./3.,0.),cmdiv(cmmul(ftocm(x*x*0.2,0.),
			  (cmsub(nc2,ftocm(2.,0.)))),(cmadd(nc2,ftocm(2.,0.))))))),
			  cmdiv(cmmul(cmmul(ftocm(4.*x5*x,0.),nc212),nc212),ftocm(9.,0.)));
			smbn = cmdiv(cmmul(cmmul(ci,ftocm(x5,0.)),(cmsub(nc2,ftocm(1.,0.)))),
			  ftocm(45.,0.));
		}

		sman1 = sman;
		smbn1 = smbn;
		t1 = 1.5;
		sman.re *= t1;
		sman.im *= t1;
		smbn.re *= t1;
		smbn.im *= t1;
		/* case n=1; note previous multiplication of sman and smbn by t1=1.5 */
		*qext = 2.*(sman.re + smbn.re);
		*qphase = 2.*(sman.im + smbn.im);
		nsqbk = -1;
		qbck.re = -2.*(sman.re - smbn.re);
		qbck.im = -2.*(sman.im - smbn.im);
		*qscat = (POW2(sman.re) + POW2(sman.im) + POW2(smbn.re) + POW2(smbn.im))/.75;
		  
		*ctbrqs = 0.0;
		n = 2;

		/************************* Major loop begins here ************************/
		while( TRUE ) 
		{
			t1 = 2.*(double)n - 1.;
			t3 = 2.*(double)n + 1.;
			wn2 = wn1;
			wn1 = wn;
			wn = cmsub(cmmul(ftocm(t1*rx,0.),wn1),wn2);
			cdum1 = a[n-1];
			cdum2 = ftocm(n*rx,0.);
			tc1 = cmadd(cmmul(cdum1,rrf),cdum2);
			tc2 = cmadd(cmmul(cdum1,nc),cdum2);
			sman = cmdiv((cmsub(cmmul(tc1,ftocm(wn.re,0.)),ftocm(wn1.re,0.))),
			  (cmsub(cmmul(tc1,wn),wn1)));
			smbn = cmdiv((cmsub(cmmul(tc2,ftocm(wn.re,0.)),ftocm(wn1.re,0.))),
			  (cmsub(cmmul(tc2,wn),wn1)));

			/* small x, n=2
			 * see bohren and huffman p 131 */
			if( x < xcut && n == 2 ) 
			{
				/* note change sign convention for m = n - ik here */
				sman = cmdiv(cmmul(cmmul(ci,ftocm(x5,0.)),(cmsub(nc2,ftocm(1.,0.)))),
					 (cmmul(ftocm(15.,0.),(cmadd(cmmul(ftocm(2.,0.),nc2),
					 ftocm(3.,0.))))));
				smbn = ftocm(0.,0.);
			}

			eqext = t3*(sman.re + smbn.re);
			*qext += eqext;
			eqpha = t3*(sman.im + smbn.im);
			*qphase += eqpha;
			nsqbk = -nsqbk;
			eqb.re = t3*(sman.re - smbn.re)*nsqbk;
			qbck.re += eqb.re;
			eqb.im = t3*(sman.im - smbn.im)*nsqbk;
			qbck.im += eqb.im;
			tx = POW2(sman.re) + POW2(sman.im) + POW2(smbn.re) + POW2(smbn.im);
			eqscat = t3*tx;
			*qscat += eqscat;
			t2 = (double)(n - 1);
			t5 = (double)n;
			t4 = t1/(t5*t2);
			t2 = (t2*(t5 + 1.))/t5;
			ectb = t2*(sman1.re*sman.re + sman1.im*sman.im + smbn1.re*smbn.re +
				 smbn1.im*smbn.im) + t4*(sman1.re*smbn1.re + sman1.im*smbn1.im);
			*ctbrqs += ectb;

			/* check convergence
			 * could decrease for large x and small m-1 in UV - X-ray; probably negligible */
			if( tx < 1.e-14 ) 
			{
				/* looks good but check relative convergence */
				eqext = fabs(eqext/ *qext);
				eqpha = fabs(eqpha/ *qphase);
				eqscat = fabs(eqscat/ *qscat);
				ectb = ( n == 2 ) ? 0. : fabs(ectb/ *ctbrqs);
				eqb.re = fabs(eqb.re/qbck.re);
				eqb.im = fabs(eqb.im/qbck.im);
				/* leave out eqb.re/im, which are sometimes least well converged */
				error = MAX4(eqext,eqpha,eqscat,ectb);
				/* put a milder constraint on eqb.re/im */
				error1 = MAX2(eqb.re,eqb.im);
				if( error < 1.e-07 && error1 < 1.e-04 )
					break;

				/* not sufficiently converged
				 *
				 * cut out after n=2 for small x, since approximation is being used */
				if( x < xcut )
					break;
			}

			smbn1 = smbn;
			sman1 = sman;
			n++;
			if( n > nmx2 ) 
			{
				*iflag = 1;
				break;
			}
		}
		/* renormalize */
		t1 = 2.*POW2(rx);
		*qext *= t1;
		*qphase *= t1;
		*qback = (POW2(qbck.re) + POW2(qbck.im))*POW2(rx);
		*qscat *= t1;
		*ctbrqs *= 2.*t1;
	}
	else 
	{
		*iflag = 2;
	}

	free( a );

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


static void anomal(double x,
		   /*@out@*/ double *qext,
		   /*@out@*/ double *qabs,
		   /*@out@*/ double *qphase,
		   /*@out@*/ double *xistar,
		   double delta,
		   double beta)
{
	/*
	 *
	 * in anomalous diffraction: qext and qabs calculated - qscatt by subtraction
	 * in rayleigh-gans:         qscatt and qabs calculated
	 * in mie:                   qext and qscatt calculated
	 *
	 */
	double xi,
	  xii;
	complex cbigk,
	  ci,
	  cw;

#	ifdef DEBUG_FUN
	fputs( "<+>anomal()\n", debug_fp );
#	endif
	/* anomalous diffraction: x>>1 and |m-1|<<1, any xi,xii
	 * original approach see Martin 1970. MN 149, 221 */
	xi = 2.*x*delta;
	xii = 2.*x*beta;
	/* xistar small is the basis for rayleigh-gans, any x, m-1 */
	*xistar = sqrt(POW2(xi)+POW2(xii));
	/* alternative approach see martin 1978 p 23 */
	ci = ftocm(0.,1.);
	cw = cmneg(cmmul(ftocm(xi,xii),ci));
	bigk(cw,&cbigk);
	*qext = 4.*cbigk.re;
	*qphase = 4.*cbigk.im;
	cw = ftocm(2.*xii,0.);
	bigk(cw,&cbigk);
	*qabs = 2.*cbigk.re;
	/* ?? put g in here - analytic version not known */

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


static void bigk(complex cw,
		 /*@out@*/ complex *cbigk)
{
	/*
	 * see martin 1978 p 23
	 */

#	ifdef DEBUG_FUN
	fputs( "<+>bigk()\n", debug_fp );
#	endif
	/* non-vax; use generic function */
	if( cmabs(cw) < 1.e-2 ) 
	{
		/* avoid severe loss of precision for small cw; expand exponential
		 * coefficients are 1/n! - 1/(n+1)! = 1/(n+1)(n-1)! ;n=2,3,4,5,6,7
		 * accurate to  (1+ order cw**6) */
		*cbigk = cmmul(cw,(cmsub(ftocm(1./3.,0.),cmmul(cw,(cmsub(ftocm(1./
		  8.,0.),cmmul(cw,(cmsub(ftocm(1./30.,0.),cmmul(cw,(cmsub(ftocm(1./
		  144.,0.),cmmul(cw,(cmsub(ftocm(1./840.,0.),cmmul(cw,ftocm(1./
		  5760.,0.)))))))))))))))));
	}
	else 
	{
		*cbigk = cmadd(ftocm(0.5,0.),cmdiv((cmsub(cmmul(cmexp(cmneg(cw)),
		  (cmadd(ftocm(1.,0.),cw))),ftocm(1.,0.))),(cmmul(cw,cw))));
	}

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


/* utility for use with mieuvx/sd */
static void ritodf(double nr,
		   double ni,
		   /*@out@*/ double *eps1,
		   /*@out@*/ double *eps2)
{

#	ifdef DEBUG_FUN
	fputs( "<+>ritodf()\n", debug_fp );
#	endif
	/* refractive index to dielectric function */
	*eps1 = nr*nr - ni*ni;
	*eps2 = 2.*nr*ni;

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


/* utility for use with mieuvx/sd */
static void dftori(/*@out@*/ double *nr,
		   /*@out@*/ double *ni,
		   double eps1,
		   double eps2)
{
	double eps;

#	ifdef DEBUG_FUN
	fputs( "<+>dftori()\n", debug_fp );
#	endif
	/* dielectric function to refractive index  */
	eps = sqrt(eps2*eps2+eps1*eps1);
	*nr = sqrt((eps+eps1)/2.);
	ASSERT( *nr > 0. );
	/* >>chng 03 jan 02, old expression for ni suffered
	 * from cancellation error in the X-ray regime, PvH */
	/* *ni = sqrt((eps-eps1)/2.); */
	*ni = eps2/(2.*(*nr));

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