/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*the prototypes are in cddefines.h and so are automatically picked up by all routines */
/*cdEXIT the routine that should be called to exit cloudy */
/*FFmtRead scan input line for free format number */
/*e2 second exponential integral */
/*caps convert input command line (through eol) to ALL CAPS */
/*ShowMe produce request to send information to GJF after a crash */
/*AnuUnit produce continuum energy in arbitrary units */
/*cap4 convert first 4 char of input line chLab into chCAP all in caps, null termination */
/*insane set flag saying that insanity has occurred */
/*lgMatch determine whether match to a keyword occurs on command line,
 * return value is 0 if no match, and position of match within string if hit */
/*fudge enter fudge factors, or some arbitrary number, with fudge command*/
/*GetElem scans line image, finds element. returns atomic number j, on C scale */
/*GetQuote get any name between double quotes off command line
 * return string as chLabel, is null terminated */
/*qip compute pow(x,n) for positive integer n through repeated sqarings */
/*NoNumb general error handler for no numbers on input line */
/*dsexp safe exponential function for doubles */
/*sexp safe exponential function */
/*TestCode set flag saying that test code is in place */
/*CodeReview - placed next to code that needs to be checked */
/*fixit - say that code needs to be fixed */
/*broken set flag saying that the code is broken, */
/*dbg_printf is a debug print routine that was provided by Peter Teuben,
 * as a component from his NEMO package.  It offers run-time specification
 * of the level of debugging */
/*qg32 32 point Gaussian quadrature, original fortran given to Gary F by Jim Lattimer */
/*TotalInsanity general error handler for something that cannot happen */
/*BadRead general error handler for trying to read data, but failing */
/*BadOpen general error handler for trying to open file, but failing */
/*MyMalloc wrapper for malloc().  Returns a good pointer or dies. */
/*MyCalloc wrapper for calloc().  Returns a good pointer or dies. */
/*BadMalloc announce memory allocation failure - put breakpoint here to find where */
/*spsort netlib routine to sort array returning sorted indices */
/*chLineLbl use information in line transfer arrays to generate a line label *
 * this label is null terminated */
/*chIonLbl use information in line array to generate a null terminated ion label in "Fe 2" */
/*csphot returns photoionization cross section from opacity stage using std pointers */
/*MyAssert a version of assert that fails gracefully */
/*RandGauss normal random variate generator */
/*MyGaussRanGauss a wrapper for RandGauss, see below */
#include <ctype.h>
#include <stdarg.h>	/* ANSI variable arg macros */
/* 
 * a set of standard routines that are widely used across the code for various
 * houeskeeping chores.  These do not do any physics, and are unlikely to
 * change over time.
 *
 * the prototypes are in cddefines.h and so are automatically picked up by all routines 
 *
 */
#include "cddefines.h"
#include "physconst.h"
#include "cddrive.h"
#include "opacity.h"
#include "rfield.h"
#include "hextra.h"
#include "hmi.h"
#include "fudgec.h"
#include "broke.h"
#include "trace.h"
#include "input.h"
#include "elementnames.h"
#include "punch.h"
#include "version.h"
#include "warnings.h"
#include "conv.h"
#include "cpu.h"

/* a version of assert that fails gracefully, and does not pause with an open window */
#ifdef _MSC_VER
/* MS compiler directive saying that cdEXIT does not return */
__declspec(noreturn) 
#elif defined(__GNUC__) || ( defined(__INTEL_COMPILER) && defined(__linux) )
 __attribute__ ((noreturn))
#endif
void MyAssert(const char *file, int line)
{
	fprintf(ioQQQ," PROBLEM An assert has been thrown, this is bad.\n");
	fprintf(ioQQQ," It happened in the file %s at line number %i\n", file, line );
	fprintf(ioQQQ," This is iteration %li zone %.2f.\n", iteration , fnzone );

	ShowMe();
	if( lgAssertFPE )
	{
		/* this flag set with SET ASSERT FPE command to throw fpe
		 * when assert is thrown, makes possible to get traceback from
		 * core file */
		float a = 1.f/(float)ZeroNum;
		fprintf(ioQQQ," when I div by 0 I get %f\n", a);
	}

	fflush( ioQQQ );
	puts( "[Stop in MyAssert]" );
	cdEXIT(EXIT_FAILURE);
}

/*AnuUnit produce continuum energy in arbitrary units */
double AnuUnit(float energy_ryd)
/*static double AnuUnit(long int ip)*/
{
	double AnuUnit_v;

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

	/* energy comes in in Ryd */
	if( energy_ryd <=0. )
	{
		/* this is insanity */
		AnuUnit_v = 0.;
	}
	else if( strcmp(punch.chConPunEnr[punch.ipConPun],"ryd ") == 0 )
	{
		/* energy in Ryd */
		AnuUnit_v = energy_ryd;
	}
	else if( strcmp(punch.chConPunEnr[punch.ipConPun],"micr") == 0 )
	{
		/* wavelength in microns */
		AnuUnit_v = RYDLAM/energy_ryd*1e-4;
	}
	else if( strcmp(punch.chConPunEnr[punch.ipConPun]," kev") == 0 )
	{
		/* energy in keV */
		AnuUnit_v = energy_ryd*EVRYD*1.e-3;
	}
	else if( strcmp(punch.chConPunEnr[punch.ipConPun]," ev ") == 0 )
	{
		/* energy in eV */
		AnuUnit_v = energy_ryd*EVRYD;
	}
	else if( strcmp(punch.chConPunEnr[punch.ipConPun],"angs") == 0 )
	{
		/* wavelength in Angstroms */
		AnuUnit_v = RYDLAM/energy_ryd;
	}
	else if( strcmp(punch.chConPunEnr[punch.ipConPun],"cent") == 0 )
	{
		/* wavelength in centimeters */
		AnuUnit_v = RYDLAM/energy_ryd*1e-8;
	}
	else if( strcmp(punch.chConPunEnr[punch.ipConPun],"wave") == 0 )
	{
		/* wavenumbers */
		AnuUnit_v = RYD_INF*energy_ryd;
	}
	else if( strcmp(punch.chConPunEnr[punch.ipConPun]," mhz") == 0 )
	{
		/* MHz */
		AnuUnit_v = RYD_INF*energy_ryd*SPEEDLIGHT/1e6;
	}
	else
	{
		fprintf( ioQQQ, " insane units in AnuUnit =%4.4s\n", 
		  punch.chConPunEnr[punch.ipConPun] );
		puts( "[Stop in AnuUnit]" );
		cdEXIT(EXIT_FAILURE);
	}

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

	return( AnuUnit_v );
}

/*ShowMe produce request to send information to GJF after a crash */
void ShowMe(void)
{

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

	/* print info if output unit is defined */
	if( ioQQQ != NULL )
	{
		/* >>chng 06 mar 02 - check if molecular but cosmic rays are ignored */
		if( (hextra.cryden == 0.) && hmi.BiggestH2 > 0.1 )
		{
			fprintf( ioQQQ, " >>> \n >>> \n >>> Cosmic rays are not included and the gas is molecular.  "
				"THIS IS KNOWN TO BE UNSTABLE.  Add cosmic rays and try again.\n >>> \n >>>\n\n");
		}
		else
		{
			fprintf( ioQQQ, "\n\n" );
			fprintf( ioQQQ, "\n\n" );
			fprintf( ioQQQ, "           vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv \n" );
			fprintf( ioQQQ, "          > Sorry, Something bad has happened.      <\n" );
			fprintf( ioQQQ, "          > Please post this on the Cloudy web site <\n" );
			fprintf( ioQQQ, "          > discussion board at www.nublado.org     <\n" );
			fprintf( ioQQQ, "          > Please send all following information:  <\n" );
			fprintf( ioQQQ, "           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ \n" );
			fprintf( ioQQQ, "\n\n" );


			fprintf( ioQQQ, " Cloudy version number is %s\n", 
				version.chVersion );
			fprintf( ioQQQ, " %s\n\n", version.chInfo );

			fprintf( ioQQQ, "%5ld warnings,%3ld cautions,%3ld temperature failures.  Messages follow.\n", 
			  warnings.nwarn, warnings.ncaun, conv.nTeFail );

			/* print the warnings first */
			cdWarnings(ioQQQ);

			/* now print the cautions */
			cdCautions(ioQQQ);

			/* now output the commands */
			cdPrintCommands(ioQQQ);

			/* if init command was present, this is the number of lines in it -
			 * if no init then still set to zero as done in cdInit */
			if( input.nSaveIni )
			{
				fprintf(ioQQQ," This input stream included an init file.\n");
				fprintf(ioQQQ," If this init file is not part of the standard Cloudy distribution\n"); 
				fprintf(ioQQQ," then I will need a copy of it too.\n");
			}
		}
	}

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

	return;
}

/*cap4 convert first 4 char of input line chLab into chCAP all in caps, null termination */
void cap4(
		char *chCAP ,	/* output string, cap'd first 4 char of chLab, */
						/* with null terminating */
		char *chLab)	/* input string ending with eol*/
{
	long int /*chr,*/ 
	  i ;

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

	/* convert 4 character string in chLab to ALL CAPS in chCAP */
	for( i=0; i < 4; i++ )
	{
		/* toupper is function in ctype that converts to upper case */
		chCAP[i] = (char)toupper( chLab[i] );
	}

	/* now end string with eol */
	chCAP[4] = '\0';


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

	return;
}

/*caps convert input command line (through eol) to ALL CAPS */
void caps(char *chCard )
{
	long int i;

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

	/* convert full character string in chCard to ALL CAPS */
	i = 0;
	while( chCard[i]!= '\0' )
	{
		chCard[i] = (char)toupper( chCard[i] );
		++i;
	}

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

	return;
}

/*e2 second exponential integral */
double e2(
	/* the argument to E2 */
	double t, 
	/* ln(x) - previously stored */
	double tln)
{
	double e2_v, 
	  remain;

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

	/* fit of second exponential integral;
	 * T is optical depth, and TLN is EXP(-t) */
	if( t < 0.3 )
	{
		remain = (1.998069357 + t*(66.4037741 + t*107.2041376))/(1. + 
		  t*(37.4009646 + t*105.0388805));

	}
	else if( t < 20. )
	{
		remain = (1.823707708 + t*2.395042899)/(1. + t*(2.488885899 - 
		  t*0.00430538));
	}

	else
	{
		remain = 1.;
	}

	e2_v = remain*tln/(2. + t);

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

	return( e2_v );
}

/*ee1 first exponential integral */
double ee1(double x)
{
	double ans, 
	  bot, 
	  top;
	static double a[6]={-.57721566,.99999193,-.24991055,.05519968,-.00976004,
	  .00107857};
	static double b[4]={8.5733287401,18.0590169730,8.6347608925,.2677737343};
	static double c[4]={9.5733223454,25.6329561486,21.0996530827,3.9584969228};

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

	/* computes the exponential integral E1(x),
	 * from abramowitz and stegun
	 * returns error condition for negative argument,
	 * returns zero when limit on sexp exceeded
	 * */

	/* error - this does not do complex numbers */
	if( x <= 0 )
	{
		fprintf( ioQQQ, " negative argument in function ee1, x<0\n" );
		puts( "[Stop in ee1]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* branch for x less than 1 */
	else if( x < 1. )
	{
		/* abs. accuracy better than 2e-7 */
/* 		ans = a[0] + a[1]*x + a[2]*x*x + a[3]*powi(x,3) + a[4]*powi(x,4) +  */
/* 		  a[5]*powi(x,5) - log(x); */
		ans = ((((a[5]*x + a[4])*x + a[3])*x + a[2])*x + a[1])*x + a[0] - log(x);
	}

	/* branch for x greater than, or equal to, one */
	else
	{
		/* abs. accuracy better than 2e-8 */
/* 		top = powi(x,4) + b[0]*powi(x,3) + b[1]*x*x + b[2]*x + b[3]; */
		top = (((x + b[0])*x + b[1])*x + b[2])*x + b[3];
/* 		bot = powi(x,4) + c[0]*powi(x,3) + c[1]*x*x + c[2]*x + c[3]; */
		bot = (((x + c[0])*x + c[1])*x + c[2])*x + c[3];
		/* 
		 *>>>chng 98 dec 17
		 * Jason's original implementation did not require the exp since
		 * the routines that used this knew if was not present.  
		 * put in for safely.
		 
		ans = top/bot/x;*/
		ans = top/bot/x*sexp(x);
	}

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

	return( ans );
}

/* same as ee1, except without factor of exp(x) in returned value	*/
double ee1_safe(double x)
{
	double ans, 
	  bot, 
	  top;
	/*static double a[6]={-.57721566,.99999193,-.24991055,.05519968,-.00976004,
	  .00107857};*/
	static double b[4]={8.5733287401,18.0590169730,8.6347608925,.2677737343};
	static double c[4]={9.5733223454,25.6329561486,21.0996530827,3.9584969228};

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

	ASSERT( x > 1. );

	/* abs. accuracy better than 2e-8 */
	/*	top = powi(x,4) + b[0]*powi(x,3) + b[1]*x*x + b[2]*x + b[3]; */
	top = (((x + b[0])*x + b[1])*x + b[2])*x + b[3];
	/*	bot = powi(x,4) + c[0]*powi(x,3) + c[1]*x*x + c[2]*x + c[3]; */
	bot = (((x + c[0])*x + c[1])*x + c[2])*x + c[3];

	ans = top/bot/x;

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

	return( ans );
}

/*FFmtRead scan input line for free format number */

/********************************************************************
 * start of nxtchr, used only by FFmtRead                           *
 ********************************************************************/
static void nxtchr(char *chr, 
	  long int *ichr, 
	  long int ipnt, 
	  char *chCard, 
	  long int last, 
	  int *lgEOL)
{

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

	/* read character chCard(IPNT:IPNT), return charac as CHR, octal as ICHR
	 * return lgEOL=.TRUE. if END reached
	 * */
	*chr = chCard[ipnt - 1];
	*ichr = (*chr);

	/* 
	 * check for end of line, either last as passed to FFmtRead,
	 * or (most likely) the end of line character 
	 */
	/* >>chng 01 feb 16, last will actually be read, had been
	 * first not read */
	if( ipnt > last || *chr=='\0' )
	{
		*lgEOL = TRUE;
	}
	else
	{
		*lgEOL = FALSE;
	}

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

	return;
}

/********************************************************************
 * start of FFmtRead                                                *
 ********************************************************************/
/*lint -e801 use of goto is deprecated */
double FFmtRead(char *chCard, 
	  long int *ipnt, 
	  /* the contents of this array element is the last that will be read */
	  long int last, 
	  int *lgEOL)
{
	int lgNFound, 
	  lgNumber;
	char chr;
	long int ichr, 
	  l1, 
	  l2;
	double FFmtRead_v, 
	  expn, 
	  PM1, 
	  value;

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

	/* FFmtRead= first number encountered after column IPNT in card image */

	/* ipnt is pointer to position within array, must be in range 1 to 
	 * longest possible line image */
	ASSERT(*ipnt > 0 );

	l1 = '0';
	l2 = '9';

L_999:
	lgNumber = FALSE;
	lgNFound = FALSE;
	PM1 = 1.;
	value = 0.;
	FFmtRead_v = 0.;
	/* have we hit the end of line? */
	nxtchr(&chr,&ichr,*ipnt,chCard,last,lgEOL);
	if( *lgEOL )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->FFmtRead()\n", debug_fp );
#		endif
		return( FFmtRead_v );
	}

	/*************************************************
	 * find start of number                          *
	 *************************************************/
L_4:
	if( chr == '.' )
		goto L_5;

	/* is the character a number? */
	if( (ichr >= l1) && (ichr <= l2) )
		goto L_1;

	*ipnt += 1;
	nxtchr(&chr,&ichr,*ipnt,chCard,last,lgEOL);
	if( *lgEOL )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->FFmtRead()\n", debug_fp );
#		endif
		return( FFmtRead_v );
	}
	goto L_4;

	/*************************************************
	 * found start of number, check if neg           *
	 *************************************************/
L_1:
	/* only do this test if we are deep in the line */
	if( *ipnt>1 )
	{
		if( chCard[*ipnt - 2] == '-' )
		  PM1 = -1.;
		lgNFound = TRUE;
	}

	/*************************************************
	 * find numerical value                          *
	 *************************************************/
L_6:
	value = 10.*value + (double)(labs(ichr-l1));
	lgNumber = TRUE;
L_7:
	*ipnt += 1;
	nxtchr(&chr,&ichr,*ipnt,chCard,last,lgEOL);
	if( *lgEOL )
	{
		FFmtRead_v = value*PM1;
		/* will be encountered on next call */
		*lgEOL = FALSE;
		
#		ifdef DEBUG_FUN
		fputs( " <->FFmtRead()\n", debug_fp );
#		endif
		return( FFmtRead_v );
	}
	if( chr == ',' )
		goto L_7;

	/* is the character a number? */
	if( (ichr >= l1) && (ichr <= l2) )
		goto L_6;
	/***************************************************/
	if( *ipnt > 0 )
	{
		if( chCard[*ipnt - 1] != '.' )
		{
			FFmtRead_v = value*PM1;
			
#			ifdef DEBUG_FUN
			fputs( " <->FFmtRead()\n", debug_fp );
#			endif
			return( FFmtRead_v );
		}
	}

	/* decimal point encountered */
L_5:
	expn = 1.;
	if( *ipnt > 1 )
	{
		if( chCard[*ipnt - 2] == '-' )
			PM1 = -1.;
	}

	/*****************************************************/
L_3:
	*ipnt += 1;
	nxtchr(&chr,&ichr,*ipnt,chCard,last,lgEOL);
	if( chr == ',' )
		goto L_3;
	if( (*lgEOL || ichr < l1) || ichr > l2 )
	{
		/* following trips if no number found (only . on line) */
		/* >>chng 01 sep 29, did not work when number started in very first col,
		 * added test for non-zero value */
		if( !lgNFound && value==0. )
			goto L_999;
		FFmtRead_v = value*PM1;
		if( lgNumber )
			*lgEOL = FALSE;
		
#		ifdef DEBUG_FUN
		fputs( " <->FFmtRead()\n", debug_fp );
#		endif
		return( FFmtRead_v );
	}
	expn *= 0.1;
	lgNFound = TRUE;
	lgNumber = TRUE;
	value += expn*(double)(ichr-l1);
	goto L_3;
	/*****************************************************/
}
/* +e801 use of goto deprecated */

static int istrstr( char *str, const char *substr )	/* return the int location of substr in str */
/*char *str, *substr;*/
{
	char *s;

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

	if( substr==NULL || *substr=='\0' )
		return( (int)strlen(str)+1 ); /*from Chuck Ditter of Infinity Graphics,6/93*/
	else if( (s=strstr(str,substr)) == NULL )
		return( 0 );	/* not present */

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

	return( (int)(s-str) + 1 );
}

/*lgMatch determine whether match to a keyword occurs on command line,
 * return value is 0 if no match, and position of match within string if hit */
int lgMatch(
	  const char *chKey, 
	  char *chCard)
{
	int lgMatch_v;
	long int i;

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

	/* cobalt blue routine that returns 0 if string is not in other,
	 * else returns position of chKey within chCard */
	i = istrstr(chCard,chKey);
	if( i == 0 )
	{
		/* did not find match, returned false */
		lgMatch_v = FALSE;
	}
	else
	{
		/* value was position with chCard, return this to lgMatch */
		lgMatch_v = (int)i;
	}

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

	return( lgMatch_v );
}

/* fudge enter fudge factors, or some arbitrary number, with fudge command
 * other sections of the code access these numbers by calling fudge
 * fudge(0) returns the first number that was entered
 * prototype for this function is in cddefines.h so that it can be used without
 * declarations 
 * fudge(-1) queries the routine for the number of fudge parameters that were entered,
 * zero returned if none */
double fudge(long int ipnt)
{
	double fudge_v;

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

	if( ipnt < 0 )
	{
		/* this is special case, return number of args */
		fudge_v = fudgec.nfudge;
	}
	else if( ipnt >= fudgec.nfudge )
	{
		fprintf( ioQQQ, " FUDGE factor not entered for array number %3ld\n", 
		  ipnt );
		puts( "[Stop in fudge]" );
		cdEXIT(EXIT_FAILURE);
	}
	else
	{
		fudge_v = fudgec.fudgea[ipnt];
	}

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

	return( fudge_v );
}

/*GetElem scans line image, finds element. returns atomic number j, on C scale, -1 if no hit */
long int GetElem(char *chCard )
{
	int i;

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

	/* find which element */

	/* >>>chng 99 apr 17, lower limit to loop had been 1, so search started with helium,
	 * change to 0 so we can pick up hydrogen.  needed for parseasserts command */
	/* find match with element name, start with helium */
	for( i=0; i<(int)LIMELM ; ++i )
	{
		if( lgMatch( elementnames.chElementNameShort[i], chCard ) )
		{

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

			/* return value is in C counting, hydrogen would be 0*/
			return( i );
		}
	}

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

	/* fall through, did not hit, return -1 as error condition */
	return (-1 );
}

/* GetQuote get any name between double quotes off command line
 * return string as chLabel, is null terminated */
void GetQuote(
		/* we will generate a label and stuff it here */
		char *chLabel,	
		/* local capd line, we will blank this out */
		char *chCard ,
		/* if true then abort if no double quotes found, 
		 * if false then return null string in this case */
		int lgABORT )
{
	char *i0,       /* pointer to first " */
	  *i1,          /* pointer to second ", name is in between */
	  *iLoc;        /* pointer to first " in local version of card in calling routine */
	size_t len;

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

	/*label within quotes returned within chLabel
	 *label in line image is set to blanks when complete */

	/* find start of string, string must begin and end with double quotes */
	/* get pointer to first quote */
	i0 = strchr( input.chOrgCard,'\"' );

	if( i0 != NULL ) 
	{
		/* get pointer to next quote */
		/*i1 = strrchr( input.chOrgCard,'\"' );*/
		i1 = strchr( i0+1 ,'\"' );
	}
	else 
	{
		i1 = NULL;
	}

	/* check that pointers were not NULL */
	/* >>chng 00 apr 27, check for i0 and i1 equal not needed anymore, by PvH */
	if( i0 == NULL || i1 == NULL )
	{
		if( lgABORT )
		{
			/* this case really do need this to work - it did not, so must abort */
			fprintf( ioQQQ, 
				" A filename or label must be specified within double quotes, but no quotes were encountered on this command.\n" );
			fprintf( ioQQQ, " Name must be surrounded by exactly two double quotes, like \"name.txt\". \n" );
			fprintf( ioQQQ, " The line image follows:\n" );
			fprintf( ioQQQ, " %s\n", input.chOrgCard);
			fprintf( ioQQQ, " Sorry\n" );
			puts( "[Stop in getquote]" );
			cdEXIT(EXIT_FAILURE);
		}
		else
		{
			/* this branch, ok if not present, return null string in that case */
			chLabel[0] = '\0';
			return;
		}
	}

	/* now copy the text inbetween quotes */
	len = (size_t)(i1-i0-1);
	strncpy(chLabel,i0+1,len);
	/* strncpy doesn't terminate the label */
	chLabel[len] = '\0';

	/* get pointer to first quote in local copy of line image in calling routine */
	iLoc = strchr( chCard ,'\"' );
	if( iLoc == NULL )
	{
		fprintf( ioQQQ, " Insanity in GetQuote - line image follows:\n" );
		fprintf( ioQQQ, " %s\n", input.chOrgCard);
		puts( "[Stop in getquote]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* >>chng 97 jul 01, blank out label once finished, to not be picked up later */
	/* >>chng 00 apr 27, erase quotes as well, so that we can find second label, by PvH */
	while( i0 <= i1 )
	{
		*i0 = ' ';
		*iLoc = ' ';
		++i0;
		++iLoc;
	}

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

	return;
}

/* another version of double ^ pos int */
/* C/C++ users journal, March 2003, p 53 */
/*qip compute pow(x,n) for positive integer n through repeated sqarings */
double qip( double x , unsigned int n)
{
	double h = 1.;

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

	while( n )
	{
		if( n&1 )
			h *= x;
		x *= x;
		n >>=1;
	}

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

	return h;
}


/* powi.c - calc x^n, where n is an integer! */

/* Very slightly modified version of power() from Computer Language, Sept. 86,
	pg 87, by Jon Snader (who extracted the binary algorithm from Donald Knuth,
	"The Art of Computer Programming", vol 2, 1969).
	powi() will only be called when an exponentiation with an integer
	exponent is performed, thus tests & code for fractional exponents were 
	removed.	-- cw lightfoot, COBALT BLUE
  COBALT BLUE has NO copyright on the functions in this file.  The source may 
  be considered public domain.
 */

/* this will only define powi if this is not an alpha running the native os */
#if !defined(__alpha) || defined(__linux)

/* #ifndef __alpha   this fcn is on an alpha, only defn if not alpha */

double powi( double x, long int n )	/* returns:  x^n */
/* x;	 base */
/* n;	 exponent */
{
	double p;	/* holds partial product */

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

	if( x == 0 )
		return( 0. );

	/* test for negative exponent */
	if( n < 0 )
	{	
		n = -n;
		x = 1/x;
	}

	p = IS_ODD(n) ? x : 1;	/* test & set zero power */

	/*lint -e704 shift right of signed quantity */
	/*lint -e720 Boolean test of assignment */
	while( n >>= 1 )
	{	/* now do the other powers */
		x *= x;			/* sq previous power of x */
		if( IS_ODD(n) )	/* if low order bit set */
			p *= x;		/*	then, multiply partial product by latest power of x */
	}
	/*lint +e704 shift right of signed quantity */
	/*lint +e720 Boolean test of assignment */

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

	return( p );
}

#	endif /*alpha unix*/
	/*		*		*		*		*/

long ipow( long m, long n )	/* returns:  m^n */
/* m;		 base */
/* n;		 exponent */
{
	long p;	/* holds partial product */

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

	if( m == 0 || (n < 0 && m > 1) )
		return( 0L );
	/* NOTE: negative exponent always results in 0 for integers!
	 * (except for the case when m==1 or -1) */

	if( n < 0 )
	{	/* test for negative exponent */
		n = -n;
		m = 1/m;
	}

	p = IS_ODD(n) ? m : 1;	/* test & set zero power */

	/*lint -e704 shift right of signed quantity */
	/*lint -e720 Boolean test of assignment */
	while( n >>= 1 )
	{	/* now do the other powers */
		m *= m;			/* sq previous power of m */
		if( IS_ODD(n) )	/* if low order bit set */
			p *= m;		/*	then, multiply partial product by latest power of m */
	}
	/*lint +e704 shift right of signed quantity */
	/*lint +e720 Boolean test of assignment */

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

	return( p );
}

/*PrintE82 - series of routines to mimic 1p, e8.2 fortran output */
/***********************************************************
 * contains the following sets of routines to get around   *
 * the MS C++ compilers unusual exponential output.        *
 * PrintEfmt <= much faster, no overhead with unix         *
 * PrintE93                                                *
 * PrintE82                                                *
 * PrintE71                                                *
 **********************************************************/

/**********************************************************/
/*
* Instead of printf'ing with %e or %.5e or whatever, call
* efmt("%e", val) and print the result with %s.  This lets
* us work around bugs in VS C 6.0.
*/
char *PrintEfmt(const char *fmt, double val /* , char *buf */) 
{
	static char buf[30]; /* or pass it in */

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

	/* create the string */
	sprintf(buf, fmt, val);

	/* we need to fix e in format if ms vs */
#	ifdef _MSC_VER
	{
		/* code to fix incorrect ms v e format.  works only for case where there is
		 * a leading space in the format - for formats like 7.1, 8.2, 9.3, 10.4, etc
		 * result will have 1 too many characters */
		char *ep , buf2[30];

		/* msvc behaves badly in different ways for positive vs negative sign vals,
		 * if val is positive must create a leading space */
		if( val >= 0.)
		{
			strcpy(buf2 , " " );
			strcat(buf2 , buf);
			strcpy( buf , buf2);
		}

		/* allow for both e and E formats */
		if ((ep = strchr(buf, 'e')) == NULL)
		{
			ep = strchr(buf, 'E');
		}

		/* ep can still be NULL if val is Inf or NaN */
		if (ep != NULL) 
		{
			/* move pointer two char past the e, to pick up the e and sign */
			ep += 2;

			/* terminate buf where the e is, *ep points to this location */
			*ep = '\0';

			/* skip next char, */
			++ep;

			/* copy resulting string to return string */
			strcat( buf, ep );
		}
	}
#	endif

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

	return buf;
}
/**********************************************************/
void PrintE82( FILE* ioOut, double value )
{
	double frac , xlog , xfloor , tvalue;
	int iExp;

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

	if( value < 0 )
	{
		fprintf(ioOut,"********");
	}
	else if( value == 0 )
	{
		fprintf(ioOut,"0.00E+00");
	}
	else
	{
		/* round number off for 8.2 format (not needed since can't be negative) */
		tvalue = value;
		xlog = log10( tvalue );
		xfloor = floor(xlog);
		/* this is now the fractional part */
		frac = tvalue*pow(10.,-xfloor);
		/*this is the possibly signed exponential part */
		iExp = (int)xfloor;
		if( frac>9.9945 )
		{
			frac /= 10.;
			iExp += 1 ;
		}
		/* print the fractional part*/
		fprintf(ioOut,"%.2f",frac);
		/* E for exponent */
		fprintf(ioOut,"E");
		/* if positive throw a + sign*/
		if(iExp>=0 )
		{
			fprintf(ioOut,"+");
		}
		fprintf(ioOut,"%.2d",iExp);
	}

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

	return;
}
/*
 *==============================================================================
 */
void PrintE71( FILE* ioOut, double value )
{
	double frac , xlog , xfloor , tvalue;
	int iExp;
	if( value < 0 )
	{
		fprintf(ioOut,"*******");
	}
	else if( value == 0 )
	{
		fprintf(ioOut,"0.0E+00");
	}
	else
	{
		/* round number off for 8.2 format (not needed since can't be negative) */
		tvalue = value;
		xlog = log10( tvalue );
		xfloor = floor(xlog);
		/* this is now the fractional part */
		frac = tvalue*pow(10.,-xfloor);
		/*this is the possibly signed exponential part */
		iExp = (int)xfloor;
		if( frac>9.9945 )
		{
			frac /= 10.;
			iExp += 1 ;
		}
		/* print the fractional part*/
		fprintf(ioOut,"%.1f",frac);
		/* E for exponent */
		fprintf(ioOut,"E");
		/* if positive throw a + sign*/
		if(iExp>=0 )
		{
			fprintf(ioOut,"+");
		}
		fprintf(ioOut,"%.2d",iExp);
	}

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

	return;
}
 
/*
 *==============================================================================
 */
void PrintE93( FILE* ioOut, double value )
{
	double frac , xlog , xfloor, tvalue;
	int iExp;
	if( value < 0 )
	{
		fprintf(ioOut,"*********");
	}
	else if( value == 0 )
	{
		fprintf(ioOut,"0.000E+00");
	}
	else
	{
		/* round number off for 9.3 format, neg numb not possilbe */
		tvalue = value;
		xlog = log10( tvalue );
		xfloor = floor(xlog);
		/* this is now the fractional part */
		frac = tvalue*pow(10.,-xfloor);
		/*this is the possibly signed exponential part */
		iExp = (int)xfloor;
		if( frac>9.99949 )
		{
			frac /= 10.;
			iExp += 1 ;
		}
		/* print the fractional part*/
		fprintf(ioOut,"%5.3f",frac);
		/* E for exponent */
		fprintf(ioOut,"E");
		/* if positive throw a + sign*/
		if(iExp>=0 )
		{
			fprintf(ioOut,"+");
		}
		fprintf(ioOut,"%.2d",iExp);
	}

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

	return;
}

/*TotalInsanity general error handler for something that cannot happen */
#ifdef _MSC_VER
/* MS compiler directive saying that TotalInsanity does not return */
__declspec(noreturn) 
#elif defined(__GNUC__) || ( defined(__INTEL_COMPILER) && defined(__linux) )
 __attribute__ ((noreturn))
#endif
void TotalInsanity(void)
{

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

	/* something that cannot happen, happened,
	 * if this message is triggered, simply place a breakpoint here
	 * and debug the error */
	fprintf( ioQQQ, " Something that cannot happen, has happened.\n" );
	fprintf( ioQQQ, " This is TotalInsanity, I live in service.c.\n" );
	ShowMe();

	puts( "[Stop in TotalInsanity]" );

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

	cdEXIT(EXIT_FAILURE);

}


/*BadRead general error handler for trying to read data, but failing */
#ifdef _MSC_VER
/* MS compiler directive saying that cdEXIT does not return */
__declspec(noreturn) 
#elif defined(__GNUC__) || ( defined(__INTEL_COMPILER) && defined(__linux) )
 __attribute__ ((noreturn))
#endif
void BadRead(void)
{

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

	/* read failed */
	fprintf( ioQQQ, " A read of internal input data has failed.\n" );
	fprintf( ioQQQ, " This is BadRead, I live in service.c.\n" );
	ShowMe();

	puts( "[Stop in BadRead]" );

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

	cdEXIT(EXIT_FAILURE);

}

/*BadOpen general error handler for trying to open file, but failing */
#ifdef _MSC_VER
/* MS compiler directive saying that cdEXIT does not return */
__declspec(noreturn) 
#elif defined(__GNUC__) || ( defined(__INTEL_COMPILER) && defined(__linux) )
 __attribute__ ((noreturn))
#endif
void BadOpen(void)
{

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

	/* read failed */
	fprintf( ioQQQ, " An attempt at opening a files has failed.\n" );
	fprintf( ioQQQ, " This is BadOpen, I live in service.c.\n" );
	ShowMe();

	puts( "[Stop in BadOpen]" );

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

	cdEXIT(EXIT_FAILURE);

}

/*NoNumb general error handler for no numbers on input line */
#ifdef _MSC_VER
/* MS compiler directive saying that cdEXIT does not return */
__declspec(noreturn) 
#elif defined(__GNUC__) || ( defined(__INTEL_COMPILER) && defined(__linux) )
 __attribute__ ((noreturn))
#endif
void NoNumb(char *chCard)
{

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

	/* general catch-all for no number when there should have been */
	fprintf( ioQQQ, " There is a problem on the following command line:\n" );
	fprintf( ioQQQ, " %s\n", chCard );
	fprintf( ioQQQ, " There should have been a number on this line.   Sorry.\n" );
	puts( "[Stop in nonumb]" );

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

	cdEXIT(EXIT_FAILURE);

 }

/*sexp safe exponential function */
double sexp(double x)
{
	double sexp_v;

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

	/* SEXP_LIMIT is 84 in cddefines.h */
	if( x < SEXP_LIMIT )
	{
		sexp_v = exp(-x);
	}
	else
	{
		sexp_v = 0.;
	}

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

	return( sexp_v );
}


/*dsexp safe exponential function for doubles */
double dsexp(double x)
{
	double dsexp_v;

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

	if( x < 680. )
	{
		dsexp_v = exp(-x);
	}
	else
	{
		dsexp_v = 0.;
	}

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

	return( dsexp_v );
}

/*TestCode set flag saying that test code is in place 
 * prototype in cddefines.h */
void TestCode(void)
{

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

	/* called if test code is in place */
	lgTestCodeCalled = TRUE;

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

	return;
}

/* sign.c - support for F77 sign() function translations by FOR_C (TM) */
/*
 *  copyright Cobalt Blue, Inc., 1988 - 1996
 *     ALL RIGHTS RESERVED
 
 * includes ssign (commented out(
 * isign (commented out)
 * 
 */
#define _ABS_(v)	( (v) < 0 ? -(v) : (v) )

#ifdef sign
#	undef sign
#endif

/* floating point sign transfer,
 * returns sign of b times value of a */
double sign(double a,double b )		
{
	return( b < 0 ? -_ABS_(a) : _ABS_(a) );
}

/* min.c - support for F77 min() function variable arg transls by FOR_C (TM) */
/*
 *  copyright Cobalt Blue, Inc., 1988 - 1995
 *     ALL RIGHTS RESERVED
 */
#ifdef _MSC_VER
#	pragma warning( disable : 4244 )/* disable bogus conv int to short warning in MS VS*/
#	pragma warning( disable : 4056 )/* disable bogus underflow warning in MS VS*/
#endif

#define VA_START(a,v)	va_start(a,v)

/* NOTE:
		All of the functions in the rest of this file are,
		VARIABLE NUMBER OF ARGRUMENT FUNCTIONS!
 */


/*lint -e777 float test equality */
double vfmin( double mn, ... )
{
	double a;
	va_list args;	/* arg ptr */

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

	if( mn == FEND )
		return( 0. );	/* NULL arg list */

	VA_START(args,mn);	/* initialize variable args */

	while( (a=va_arg(args,double)) != FEND )
		mn = (mn < a) ? mn : a;

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

	return( mn );
}
/*lint +e777 float test equality */

/* max.c - support for F77 max() function variable arg transls by FOR_C (TM) */
/*
 *  copyright Cobalt Blue, Inc., 1988 - 1995
 *     ALL RIGHTS RESERVED
 */

/* NOTE: These functions guarantee max() problems won't occur due to
		side effects - at the cost of an extra function call. */
#ifdef smax
#	undef smax
#endif
short smax( short a, short b)
{
	return( a > b ? a : b );
}

#if 0
#ifdef max
#	undef max
#endif
long max( long a, long b )
{
	return( a > b ? a : b );
}
#endif

/*lint -e777 float test equality */
double vfmax( double mx, ... )
{
	double a;

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

	va_list args;	/* arg ptr */

	if( mx == FEND )
		return( 0. );	/* NULL arg list */

	VA_START(args,mx);	/* initialize variable args */

	while( (a=va_arg(args,double)) != FEND )
		mx = (mx > a) ? mx : a;

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

	return( mx );
}
/*lint +e777 float test equality */

/*broken set flag saying that the code is broken, */
void broken(void)
{

	broke.lgBroke = TRUE;

	return;
}

/*fixit say that code needs to be fixed */
void fixit(void)
{

	broke.lgFixit = TRUE;

	return;
}

/*CodeReview placed next to code that needs to be checked */
void CodeReview(void)
{

	broke.lgCheckit = TRUE;

	return;
}

/* dbg_printf is a debug print routine that was provided by Peter Teuben,
 * as a component from his NEMO package.  It offers run-time specification
 * of the level of debugging */
void dbg_printf(int debug, const char *fmt, ...)
{
	va_list ap;

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

	/* print this debug message? (debug_level not currently used)*/
	if (debug <= trace.debug_level) 
	{		
		va_start(ap, fmt);	

		/*lint -e534 ignore value of return */
		vfprintf(ioQQQ, fmt, ap);
		/*lint +e534 ignore value of return */
		/* drain ioQQQ */
		fflush(ioQQQ);
		va_end(ap);
	}
#	ifdef DEBUG_FUN
	fputs( " <->dbg_printf()\n", debug_fp );
#	endif

	return;
}


/*qg32 32 point Gaussian quadrature, originally given to Gary F by Jim Lattimer */
double qg32(
	double xl, /*lower limit to integration range*/
	double xu, /*upper limit to integration range*/
	/*following is the pointer to the function that will be evaulated*/
	double (*fct)(double) )
{
	double a, 
	  b, 
	  c, 
	  y;

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


	/********************************************************************************
	 *                                                                              *
	 *  32-point Gaussian quadrature                                                *
	 *  xl  : the lower limit of integration                                        *
	 *  xu  : the upper limit                                                       *
	 *  fct : the (external) function                                               *
	 *  returns the value of the integral                                           *
	 *                                                                              *
	 * simple call to integrate sine from 0 to pi                                   *
	 * double agn = qg32( 0., 3.141592654 ,  sin );                                 *
	 *                                                                              *
	 *******************************************************************************/

	a = .5*(xu + xl);
	b = xu - xl;
	c = .498631930924740780*b;
	y = .35093050047350483e-2*((*fct)(a+c) + (*fct)(a-c));
	c = b*.49280575577263417;
	y += .8137197365452835e-2*((*fct)(a+c) + (*fct)(a-c));
	c = b*.48238112779375322;
	y += .1269603265463103e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.46745303796886984;
	y += .17136931456510717e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.44816057788302606;
	y += .21417949011113340e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.42468380686628499;
	y += .25499029631188088e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.3972418979839712;
	y += .29342046739267774e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.36609105937014484;
	y += .32911111388180923e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.3315221334651076;
	y += .36172897054424253e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.29385787862038116;
	y += .39096947893535153e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.2534499544661147;
	y += .41655962113473378e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.21067563806531767;
	y += .43826046502201906e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.16593430114106382;
	y += .45586939347881942e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.11964368112606854;
	y += .46922199540402283e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.7223598079139825e-1;
	y += .47819360039637430e-1*((*fct)(a+c) + (*fct)(a-c));
	c = b*.24153832843869158e-1;
	y = b*(y + .482700442573639e-1*((*fct)(a+c) + (*fct)(a-c)));

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

	/* the answer */
	return( y );
}

/*spsort netlib routine to sort array returning sorted indices */
void spsort(
	  /* input array to be sorted */
	  float x[], 
	  /* number of values in x */
	  long int n, 
	  /* permutation output array */
	  long int iperm[], 
	  /* flag saying what to do - 1 sorts into increasing order, not changing
	   * the original vector, -1 sorts into decreasing order. 2, -2 change vector */
	  int kflag, 
	  /* error condition, should be 0 */
	  int *ier)
#if 0
void spsort(float x[], 
	  long int n, 
	  long int iperm[], 
	  int kflag, 
	  int *ier)
#endif
{
	/*
	 ****BEGIN PROLOGUE  SPSORT
	 ****PURPOSE  Return the permutation vector generated by sorting a given
	 *            array and, optionally, rearrange the elements of the array.
	 *            The array may be sorted in increasing or decreasing order.
	 *            A slightly modified quicksort algorithm is used.
	 ****LIBRARY   SLATEC
	 ****CATEGORY  N6A1B, N6A2B
	 ****TYPE      SINGLE PRECISION (SPSORT-S, DPSORT-D, IPSORT-I, HPSORT-H)
	 ****KEY WORDS NUMBER SORTING, PASSIVE SORTING, SINGLETON QUICKSORT, SORT
	 ****AUTHOR  Jones, R. E., (SNLA)
	 *           Rhoads, G. S., (NBS)
	 *           Wisniewski, J. A., (SNLA)
	 ****DESCRIPTION
	 *
	 *   SPSORT returns the permutation vector IPERM generated by sorting
	 *   the array X and, optionally, rearranges the values in X.  X may
	 *   be sorted in increasing or decreasing order.  A slightly modified
	 *   quicksort algorithm is used.
	 *
	 *   IPERM is such that X(IPERM(I)) is the Ith value in the rearrangement
	 *   of X.  IPERM may be applied to another array by calling IPPERM,
	 *   SPPERM, DPPERM or HPPERM.
	 *
	 *   The main difference between SPSORT and its active sorting equivalent
	 *   SSORT is that the data are referenced indirectly rather than
	 *   directly.  Therefore, SPSORT should require approximately twice as
	 *   long to execute as SSORT.  However, SPSORT is more general.
	 *
	 *   Description of Parameters
	 *      X - input/output -- real array of values to be sorted.
	 *          If ABS(KFLAG) = 2, then the values in X will be
	 *          rearranged on output; otherwise, they are unchanged.
	 *      N - input -- number of values in array X to be sorted.
	 *      IPERM - output -- permutation array such that IPERM(I) is the
	 *              index of the value in the original order of the
	 *              X array that is in the Ith location in the sorted
	 *              order.
	 *      KFLAG - input -- control parameter:
	 *            =  2  means return the permutation vector resulting from
	 *                  sorting X in increasing order and sort X also.
	 *            =  1  means return the permutation vector resulting from
	 *                  sorting X in increasing order and do not sort X.
	 *            = -1  means return the permutation vector resulting from
	 *                  sorting X in decreasing order and do not sort X.
	 *            = -2  means return the permutation vector resulting from
	 *                  sorting X in decreasing order and sort X also.
	 *      IER - output -- error indicator:
	 *          =  0  if no error,
	 *          =  1  if N is zero or negative,
	 *          =  2  if KFLAG is not 2, 1, -1, or -2.
	 ****REFERENCES  R. C. Singleton, Algorithm 347, An efficient algorithm
	 *                 for sorting with minimal storage, Communications of
	 *                 the ACM, 12, 3 (1969), pp. 185-187.
	 ****ROUTINES CALLED  XERMSG
	 ****REVISION HISTORY  (YYMMDD)
	 *   761101  DATE WRITTEN
	 *   761118  Modified by John A. Wisniewski to use the Singleton
	 *           quicksort algorithm.
	 *   870423  Modified by Gregory S. Rhoads for passive sorting with the
	 *           option for the rearrangement of the original data.
	 *   890620  Algorithm for rearranging the data vector corrected by R.
	 *           Boisvert.
	 *   890622  Prologue upgraded to Version 4.0 style by D. Lozier.
	 *   891128  Error when KFLAG.LT.0 and N=1 corrected by R. Boisvert.
	 *   920507  Modified by M. McClain to revise prologue text.
	 *   920818  Declarations section rebuilt and code restructured to use
	 *           IF-THEN-ELSE-ENDIF.  (SMR, WRB)
	 ****END PROLOGUE  SPSORT
	 *     .. Scalar Arguments ..
	 */
	long int i, 
	  ij, 
	  il[21], 
	  indx, 
	  indx0, 
	  istrt, 
	  istrt_, 
	  iu[21], 
	  j, 
	  k, 
	  kk, 
	  l, 
	  lm, 
	  lmt, 
	  m, 
	  nn;
	float r, 
	  ttemp;

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

	/*     .. Array Arguments .. */
	/*     .. Local Scalars .. */
	/*     .. Local Arrays .. */
	/*     .. External Subroutines .. */
	/*     .. Intrinsic Functions .. */
	/****FIRST EXECUTABLE STATEMENT  SPSORT */
	*ier = 0;
	nn = n;
	if( nn < 1 )
	{
		*ier = 1;
		
#		ifdef DEBUG_FUN
		fputs( " <->spsort()\n", debug_fp );
#		endif
		return;
	}
	else
	{
		kk = labs(kflag);
		if( kk != 1 && kk != 2 )
		{
			*ier = 2;
			
#			ifdef DEBUG_FUN
			fputs( " <->spsort()\n", debug_fp );
#			endif
			return;
		}
		else
		{

			/* Initialize permutation vector to index on f scale
			 * */
			for( i=0; i < nn; i++ )
			{
				iperm[i] = i+1;
			}

			/*     Return if only one value is to be sorted
			 * */
			if( nn == 1 )
			{ 
				--iperm[0];
#				ifdef DEBUG_FUN
				fputs( " <->spsort()\n", debug_fp );
#				endif
				return;
			}

			/*     Alter array X to get decreasing order if needed
			 * */
			if( kflag <= -1 )
			{
				for( i=0; i < nn; i++ )
				{
					x[i] = -x[i];
				}
			}

			/*     Sort X only
			 * */
			m = 1;
			i = 1;
			j = nn;
			r = .375e0;
		}
	}

	while( TRUE )
	{
		if( i == j )
			goto L_80;
		if( r <= 0.5898437e0 )
		{
			r += 3.90625e-2;
		}
		else
		{
			r -= 0.21875e0;
		}

L_40:
		k = i;

		/*     Select a central element of the array and save it in location L
		 * */
		ij = i + (long)((j-i)*r);
		lm = iperm[ij-1];

		/*     If first element of array is greater than LM, interchange with LM
		 * */
		if( x[iperm[i-1]-1] > x[lm-1] )
		{
			iperm[ij-1] = iperm[i-1];
			iperm[i-1] = lm;
			lm = iperm[ij-1];
		}
		l = j;

		/*     If last element of array is less than LM, interchange with LM
		 * */
		if( x[iperm[j-1]-1] < x[lm-1] )
		{
			iperm[ij-1] = iperm[j-1];
			iperm[j-1] = lm;
			lm = iperm[ij-1];

			/*        If first element of array is greater than LM, interchange
			 *        with LM
			 * */
			if( x[iperm[i-1]-1] > x[lm-1] )
			{
				iperm[ij-1] = iperm[i-1];
				iperm[i-1] = lm;
				lm = iperm[ij-1];
			}
		}

		/*     Find an element in the second half of the array which is smaller
		 *     than LM
		 * */
		while( TRUE )
		{
			l -= 1;
			if( x[iperm[l-1]-1] <= x[lm-1] )
			{

				/*     Find an element in the first half of the array which is greater
				 *     than LM
				 * */
				while( TRUE )
				{
					k += 1;
					if( x[iperm[k-1]-1] >= x[lm-1] )
						break;
				}

				/*     Interchange these elements
				 * */
				if( k > l )
					break;
				lmt = iperm[l-1];
				iperm[l-1] = iperm[k-1];
				iperm[k-1] = lmt;
			}
		}

		/*     Save upper and lower subscripts of the array yet to be sorted
		 * */
		if( l - i > j - k )
		{
			il[m-1] = i;
			iu[m-1] = l;
			i = k;
			m += 1;
		}
		else
		{
			il[m-1] = k;
			iu[m-1] = j;
			j = l;
			m += 1;
		}

L_90:
		if( j - i >= 1 )
			goto L_40;
		if( i == 1 )
			continue;
		i -= 1;

		while( TRUE )
		{
			i += 1;
			if( i == j )
				break;
			lm = iperm[i];
			if( x[iperm[i-1]-1] > x[lm-1] )
			{
				k = i;

				while( TRUE )
				{
					iperm[k] = iperm[k-1];
					k -= 1;

					if( x[lm-1] >= x[iperm[k-1]-1] )
						break;
				}
				iperm[k] = lm;
			}
		}

		/*     Begin again on another portion of the unsorted array
		 * */
L_80:
		m -= 1;
		if( m == 0 )
			break;
		/*lint -e644 not explicitly initialized */
		i = il[m-1];
		j = iu[m-1];
		/*lint +e644 not explicitly initialized */
		goto L_90;
	}

	/*     Clean up
	 * */
	if( kflag <= -1 )
	{
		for( i=0; i < nn; i++ )
		{
			x[i] = -x[i];
		}
	}

	/*     Rearrange the values of X if desired
	 * */
	if( kk == 2 )
	{

		/*        Use the IPERM vector as a flag.
		 *        If IPERM(I) < 0, then the I-th value is in correct location
		 * */
		for( istrt=1; istrt <= nn; istrt++ )
		{
			istrt_ = istrt - 1;
			if( iperm[istrt_] >= 0 )
			{
				indx = istrt;
				indx0 = indx;
				ttemp = x[istrt_];
				while( iperm[indx-1] > 0 )
				{
					x[indx-1] = x[iperm[indx-1]-1];
					indx0 = indx;
					iperm[indx-1] = -iperm[indx-1];
					indx = labs(iperm[indx-1]);
				}
				x[indx0-1] = ttemp;
			}
		}

		/*        Revert the signs of the IPERM values
		 * */
		for( i=0; i < nn; i++ )
		{
			iperm[i] = -iperm[i];
		}
	}

	/* gary add move iperm from f to c scale */
	for( i=0; i < nn; i++ )
	{
		--iperm[i];
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->spsort()\n", debug_fp );
#	endif

	return;
}

/*MyMalloc wrapper for malloc().  Returns a good pointer or dies. 
 * memory is filled with NaN
 * >>chng 05 dec 14, do not set to NaN since tricks debuggerws 
 * routines within code do not call this or malloc, but rather MALLOC
 * which is resolved into MyMalloc or malloc depending on whether 
 * NDEBUG is set by the compiler to indicate "not debugging",
 * in typical negative C style */
void *MyMalloc( 
	/*use same type as library function MALLOC*/ 
	size_t size ,
	const char *chFile, 
	int line
	)
{
	void *ptr;

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

	ASSERT( size > 0 );

	/* debug branch for printing malloc args */
	{
		/*@-redef@*/
		enum{DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC)
		{
			static long int kount=0, nTot=0;
			nTot += (long)size;
			fprintf(ioQQQ,"%li\t%li\t%li\n", 
				kount , 
				(long)size , 
				nTot );
			++kount;
		}
	}

	if( ( ptr = malloc( size ) ) == NULL )
	{
		fprintf(ioQQQ,"DISASTER MyMalloc could not allocate %lu bytes.  Exit in MyMalloc.",
			(unsigned long)size );
		fprintf(ioQQQ,"MyMalloc called from file %s at line %i.\n",
			chFile , line );
		puts( "[Stop in MyMalloc]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* flag -DNOINIT will turn off this initialization which can fool valgrind/purify */
#	if !defined(NDEBUG) && !defined(NOINIT)
{
#	ifdef __mips

	/* definition of signaling and quiet NaN is reversed on MIPS
	 * hence this will generate signaling NaN for both float and double */
	memset( ptr, 0xff, size );

#	else

	size_t nFloat = size/4;
	size_t nDouble = size/8;
	int32 *ptr32 = (int32 *)ptr;
	/* >>chng 04 feb 03, fill memory with invalid numbers, PvH */
	/* on IA32/AMD64 processors this will generate NaN's for both float and double;
	 * on most other (modern) architectures it is likely to do the same... */
	/* >>chng 05 dec 14, change code to generate signaling NaN's for most cases (but not all!) */
	if( size == nDouble*8 )
	{
		/* this could be an array of doubles as well as floats -> we will hedge our bets
		 * we will fill the array with a pattern that is interpreted as all signaling
		 * NaN's for doubles, and alternating signaling and quiet NaN's for floats:
		 * byte offset:  0      4      8     12     16
		 * double        | SNaN        | SNan        |
		 * float         | SNaN | QNaN | SNan | QNaN |  (little-endian, e.g. Intel, AMD, alpha)
		 * float         | QNaN | SNaN | QNan | SNaN |  (big-endian, e.g. Sparc, Mac, MIPS) */
		if( CPU_IS_BIG_ENDIAN )
			while( nFloat > 0 )
			{
				--nFloat;
				ptr32[nFloat] = ( (nFloat%2) == 1 ) ? 0xffbfffff : 0xfff7ffff;
			}
		else if( CPU_IS_LITTLE_ENDIAN )
			while( nFloat > 0 )
			{
				--nFloat;
				ptr32[nFloat] = ( (nFloat%2) == 0 ) ? 0xffbfffff : 0xfff7ffff;
			}
		else /* unusual CPU type, pattern for SNaN is unknown... */
			memset( ptr, 0xff, size );
	}
	else if( size == nFloat*4 )
	{
		/* this could be an arrays of floats, but not doubles -> init to all float SNaN */
		while( nFloat > 0 )
			ptr32[--nFloat] = 0xffbfffff;
	}
	else
	{
		memset( ptr, 0xff, size );
	}

#	endif /* defined(__mips) */
}
#	endif /* !defined(NDEBUG) && !defined(NOINIT) */

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

	return ptr;
}


/* wrapper for calloc().  Returns a good pointer or dies. 
 * routines within code do not call this or malloc, but rather CALLOC
 * which is resolved into MyCalloc or calloc depending on whether 
 * NDEBUG is set in cddefines. \h */
void *MyCalloc( 
	/*use same type as library function CALLOC*/ 
	size_t num ,
	size_t size )
{
	void *ptr;

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

	ASSERT( size > 0 );

	/* debug branch for printing malloc args */
	{
		/*@-redef@*/
		enum{DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC)
		{
			static long int kount=0;
			fprintf(ioQQQ,"%li\tcall\t%li\tbytes\n", kount, 
				(long)size );
			++kount;
		}
	}

	if( ( ptr = calloc( num , size ) ) == NULL )
	{
		fprintf(ioQQQ,"MyCalloc could not allocate %lu bytes.  Exit in MyCalloc.",
			(unsigned long)size );
		puts( "[Stop in MyCalloc]" );
		cdEXIT(EXIT_FAILURE);
	}

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

	return ptr;
}

/* announce memory allocation failure - put breakpoint here to find where */
#ifdef _MSC_VER
/* MS compiler directive saying that cdEXIT does not return */
__declspec(noreturn) 
#elif defined(__GNUC__) || ( defined(__INTEL_COMPILER) && defined(__linux) )
 __attribute__ ((noreturn))
#endif
void BadMalloc(void)
{
	fprintf(ioQQQ,"BadMalloc: a memory allocation failed.");
	puts( "[Stop in BadMalloc]" );
	cdEXIT(EXIT_FAILURE);
}

/* function to facilitate addressing opacity array */
double csphot(
	/* INU is array index pointing to frequency where opacity is to be evaluated
	 * on f not c scale */
	long int inu, 
	/* ITHR is pointer to threshold*/
	long int ithr, 
	/* IOFSET is offset as defined in opac0*/
	long int iofset)
{
	double csphot_v;

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

	csphot_v = opac.OpacStack[inu-ithr+iofset-1];

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

	return( csphot_v );
}

/*RandGauss normal Gaussian random number generator
 * the user must call srand to set the seed before using this routine.
 * the return value is the log10 of the normal distribution scale factor,
 * so this routine returns log normal.
 * The convention is for srand to be called when the command setting 
 * the noise is parsed */
/*
NB - following from Ryan Porter:
I discovered that I unintentionally created an antisymmetric skew in my 
monte carlo.  RandGauss is symmetric in log space, which means it is not 
symmetric in linear space.  But to get the right standard deviation you 
have to take 10^x, where x is the return from RandGauss.  The problem is 
10^x will happen less frequently than 10^-x, so without realizing it, the 
average "tweak" to every piece of atomic data in my monte carlo run was 
not 1.0 but something greater than 1.0, causing every single line to have 
an average monte carlo emissivity greater than the regular value.  Not 
sure how much this will affect the results I presented yesterday, but 
Gargi wasn't aware of this general issue.  Anyplace that takes 
10^RandGauss() needs to be adjusted if what is intended is +/- x, for 
example.
*/
double RandGauss(
	/* mean value */
	double xMean, 
	/*standard deviation s */
	double s )
{
	double  x1, x2, w, yy1;
	static double yy2=-BIGDOUBLE;
	static int use_last = FALSE;

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

	if( use_last )
	{
		yy1 = yy2;
		use_last = FALSE;
	}
	else
	{
		do {
			x1 = 2.*((double)rand()/ (double)RAND_MAX) - 1.;
			x2 = 2.*((double)rand() / (double)RAND_MAX) - 1.;
			w = x1 * x1 + x2 * x2;
		} while ( w >= 1.0 );

		w = sqrt((-2.0*log(w))/w);
		yy1 = x1 * w;
		yy2 = x2 * w;
		use_last = TRUE;
	}

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

	return( xMean + yy1 * s );
}

double MyGaussRand( double PctUncertainty )
{
	double StdDev;
	double result;

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

	ASSERT( PctUncertainty < 0.5 );
	/* We want this "percent uncertainty to represent two standard deviations */
	StdDev = 0.5 * PctUncertainty;

	do
	{
		/*result = pow( 10., RandGauss( 0., logStdDev ) );*/
		result = 1.+RandGauss( 0., StdDev );
	}
	/* This will give us a result that is less than or equal to the
	 * percent uncertainty about 98% of the time. */
	while( (result < 1.-PctUncertainty) || (result > 1+PctUncertainty) );
	
	ASSERT( result>0. && result<2. );

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

	return result;
}

/*plankf evaluate Planck function for any cell at current electron temperature */
double plankf(long int ip)
{
	double plankf_v;

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

	/* evaluate Planck function
	 * argument is pointer to cell energy in ANU
	 * return photon flux for cell IP */
	if( rfield.ContBoltz[ip] <= 0. )
	{
		plankf_v = 1e-35;
	}
	else
	{
		plankf_v = 6.991e-21*POW2(FR1RYD*rfield.anu[ip])/
			(1./rfield.ContBoltz[ip] - 1.)*FR1RYD*4.;
	}

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

	return( plankf_v );
}

