/* 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 */
/* these headers are needed by all files */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <string.h>
#include <float.h>
#include <limits.h>
#include <time.h>

#ifdef  __cplusplus
#include <string>
using namespace std;
#endif

/* some definitions for integers with a guaranteed number of bits */

#ifndef INT16_MAX
#define INT16_MAX 32767
#endif
#ifndef INT16_MIN
#define INT16_MIN (-INT16_MAX - 1)
#endif

#if INT_MAX == INT16_MAX
typedef int int16;
#elif SHRT_MAX == INT16_MAX
typedef short int int16;
#else
#error failed to define int16, please report this to gary@pa.uky.edu
#endif

#ifndef UINT16_MAX
#define UINT16_MAX 65535
#endif

#if UINT_MAX == UINT16_MAX
typedef unsigned int uint16;
#elif USHRT_MAX == UINT16_MAX
typedef unsigned short int uint16;
#else
#error failed to define uint16, please report this to gary@pa.uky.edu
#endif

#ifndef INT32_MAX
#define INT32_MAX 2147483647L
#endif
#ifndef INT32_MIN
#define INT32_MIN (-INT32_MAX - 1)
#endif

#if LONG_MAX == INT32_MAX
typedef long int int32;
#elif INT_MAX == INT32_MAX
typedef int int32;
#else
#error failed to define int32, please report this to gary@pa.uky.edu
#endif

#ifndef UINT32_MAX
#define UINT32_MAX 4294967295UL
#endif

#if ULONG_MAX == UINT32_MAX
typedef unsigned long int uint32;
#elif UINT_MAX == UINT32_MAX
typedef unsigned int uint32;
#else
#error failed to define uint32, please report this to gary@pa.uky.edu
#endif

/* to avoid errors introduced by C's infamous double-negative logic,
 * this uses NDEBUG (the ANSI std macro used to tell assert that
 * we are not debugging) to define DEBUG */
#ifndef NDEBUG
#	define DEBUG
#else
#	undef DEBUG
#endif

#ifdef DMALLOC
#include <dmalloc.h>
#endif

/* as the code is developed this should become part of the design for future
 * true parallel processing.  This flag should be false if it is not necessary
 * for a given instance of a routine to keep separate memory.  This is the case in
 * a serial machine, or an MPI run where there are n instances of the code running
 * on n processors.  This should be TRUE if an instance of a routine must have its
 * own unique memory, so that other instances of the routine, which may be running
 * in parallel on another processor, does not share the same memory.  */
#ifndef PARALLEL_MODE
#define PARALLEL_MODE FALSE
#endif

/* use special debug version of malloc for debug runs - it includes assert arg is positive
 * and dies if cannot allocate space*/
/* use regular malloc with the dmalloc library so that it can see where memory is allocated */
#if defined(NDEBUG) || defined(DMALLOC)
#	define MALLOC(exp) (malloc(exp))
#else
#	define MALLOC(exp) (MyMalloc(exp,__FILE__, __LINE__))
#endif

/* now special version of calloc for debug runs - it includes assert arg is positive
 * and dies if cannot allocate space*/
/* use regular calloc with the dmalloc library so that it can see where memory is allocated */
#if defined(NDEBUG) || defined(DMALLOC)
#	define CALLOC calloc
#else
#	define CALLOC MyCalloc 
#endif


/* the number of iso sequences now in the code */
#define	NISO	2

/* stopping temperature */
#define	TEND	4000.

#ifndef TRUE
#define TRUE	1
#endif

#ifndef FALSE
#define FALSE	0
#endif

/* a version of assert that fails gracefully */
#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
/*@noreturn@*/
void MyAssert(const char *file, int line);

/* a cloudy version of the assert macro */
#undef  ASSERT
#ifndef STD_ASSERT
#	ifdef NDEBUG
#		define ASSERT(exp) ((void)0)
#	else
/* a version of assert that fails gracefully */
#		define ASSERT(exp) if (!(exp)) \
		MyAssert(__FILE__, __LINE__)
#	endif
#else
#	define ASSERT(exp)	(assert(exp))
#endif

/**************************************************************************
 * the following are prototypes for some routines that are part of the
 * debugging process - they come and go in any particular sub.  
 * it is not necessary to declare them when used since they are defined here
 **************************************************************************/

void cdDefines(void);

/* fudge enter fudge factors, or some arbitrary number, with fudge command*/
double fudge(long int ipnt);

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

/*fixit set flag saying that this code needs attention, but is not broken,
 * code is in broken.c */
void fixit(void);

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

/*testcde set flag saying that test code is in place */
void TestCode(void);

/*csphot returns photoionization cross section from opacity stage using std pointers */
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);

/*MyMalloc wrapper for malloc().  Returns a good pointer or dies. */
void *MyMalloc( 
	/*use same type as library function malloc*/ 
	size_t size ,
	const char *file, int line
	);

/* normal random variate generator */
double RandGauss(
	/* mean value */
	double xMean, 
	/*standard deviation s */
	double s );

/*A custom wrapper for RandGauss than truncates at two standard deviations. */
double MyGaussRand( double PctUncertainty );

/*MyCalloc wrapper for calloc().  Returns a good pointer or dies. */
void *MyCalloc( 
	/*use same type as library function CALLOC*/ 
   size_t num ,
   size_t size
	);

/* announce memory allocation failure and exit - 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
/*@noreturn@*/
void BadMalloc(void);

/* FILENAME_PATH_LENGTH is the size of the string that holds the path.  The longest
 * string that can be held is one less than this, due to the end
 * end of string sentinal in C.  Increase this is a larger string
 * is needed to hold the path on your system */
#define FILENAME_PATH_LENGTH 200 

/* twice the above, so that we can add file name to end of full path */
#define FILENAME_PATH_LENGTH_2 (FILENAME_PATH_LENGTH*2) 

/* this is limit to longest line of information that is scanned in,
 * end of line char is actually at this +1, dim of vector is this + 1 
 * all routines that scan in information should use this to dim vars */ 
#define INPUT_LINE_LENGTH 200

/*******************************************************************
 *
 * this structure defines the main emission line class.  
 *
 *********************************************************************/
typedef struct 
{
	/* all of these are set to dangerous values by EmLineJunk */

	/* type of redistribution function, 
	-1 complete redis Dopper core only, 
	 0 insanity
	 1 incomplete redistirbution with damping wings, 
	 2 complete redistribution with damping wings,
	*/
	int iRedisFun ;

	/* index for line within continuum array,
	 * this is on the f, not c, scale,
	 * negative ipCont means this is not a real line
	 * sentinal for fake line */
	int ipCont;

	/* index for line within fine continuum array, for line overlap */
	long int ipFine;

	 /* ion stage of element, 1 for atom, 2 ion, etc */
	int IonStg;

	 /* atomic number of element, 1 for H, 2 for He, etc */
	int nelem;

	/* line optical depth [Napier] to continuum source through illuminated face */
	float TauCon ;

	/* collision rate coefficient, [cm^3 s-1], upper to lower */
	float ColUL;

	/* inward line optical depth [Napier], from illuminated face of cloud to current position */
	float TauIn;

	/* total line optical depth [Napier] through cloud */
	float TauTot ;

	/* inward fraction [dimensionless] of total line emission*/
	double FracInwd;

	/* continuum pumping rate [s-1] from lower to upper level, A*occ num * g_up/g_lo,
	 * this is evaluated in RTMakeStat and RTMakeWind, which are called by HydroPEsc, RT_line_all */
	double pump ;

	/* line intensity per unit time and vol [erg s-1 cm-3] */
	double xIntensity ;

	/* number of photons emitted in the line per unit time and vol [photons s-1 cm-3] */
	double phots;

	/* gf value [dimensionless] */
	float gf;

	/* escape prob [dimensionless] */
	double Pesc;

	/* electron scattering escape prob [dimensionless] */
	double Pelec_esc;

	/* destruction probs [dimensionless] */
	double Pdest;

	/* damping constant is dampXvel divided by line width 
	 * units are velocity, since becomes dimensionless when div by line width in cm/2
	 * [cm s-1] */
	float dampXvel ;

	/* [dimensionless] damping constant */
	float damp;

	/* cooling and heating due to collisional excitation [erg s-1 cm-3] */
	double cool , heat ;

	/* [dimensionless] ratio of collisional to radiative excitation, C_lu / ( C_lu + pump )*/
	float ColOvTot ;

	/* [dimensionless] collision strength or rates for transition */
	double cs;

	/* wavelentgh, usually in Angstroms, used for printout, can be any units */
	float WLAng;

	/* transition energy in degrees kelvin*/
	float EnergyK;

	/* transition energy in ergs */
	float EnergyErg;

	/* transition energy in wavenumbers */
	float EnergyWN;

	/* atomic constant part of line opacity per atom, 
	 * divide by line width in cm/s, DoppVel.doppler[nelem],
	 * to get line center opacity per atom, or absorption cross section, with units cm^2
	 * multiply by PopOpc to get PopOpc/dopper width, the true opacity (cm-1),
	 * then by length to get optical depth */
	float opacity;

	/* lower and upper statistical weights [dimensionless] */
	float gLo , gHi ;

	/* populations of lower and upper levels [cm-3] */
	/* >>chng 01 nov 24, from float to double, need the precision for hyperfine structure lines */
	double PopLo , PopHi;

	/* population that enters net opacity after correction for stimulated emission [cm-3] */
	double PopOpc ;

	/* transition prob, Einstein A upper to lower [s-1] */
	double Aul ;

	/* ots rate [cm-3 s-1] */
	double ots;

} EmLine;

/*********************************************************************/

/* the Cloudy exit handler */
#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
/*@noreturn@*/
void cdEXIT(int);

/**************************************************************************
 * these are variables and pointers for output from the code, used everywhere
 * declared extern here, and definition is in cddefines.c
 **************************************************************************/

/* ioQQQ is the file hand to the output file itself, 
 * ioQQQ is set to stdout by default in cdInit, 
 * and is reset to anything else by calling cdOutp */
extern FILE *ioQQQ;

extern FILE *ioStdin;

extern FILE *ioMAP;

/* we shall write errors to this file, it is set
 * to stderr in cdInit */
extern FILE* ioPrnErr; 

/* this says whether (true) or not (false) this is a 32 bit cpu */
extern int lgBit32;

/* this set true when abort sequence is initiated - serious meltdown is happening */
extern int lgAbort;

/* flag lgTestIt turned on if routine TestCode ever called, only generates
 * comment that test code is in place */
extern int lgTestCodeCalled; 

/* flag lgTestOn set true with SET TEST command
 * for some test code to be run somewhere */
extern int lgTestCodeEnabled;

/* this is flag saying whether to print errors to 
 * standard error output */
extern int lgPrnErr;

/* this tells the code to throw an fpe when an assert is thrown,
 * set true with SET ASSERT FPE command */
extern int lgAssertFPE;

/* nzone is zone counter, incremented in routine cloudy */
extern long int nzone;

/* this is nzone + conv.nPres2Ioniz/100 in ConvBase */
extern double fnzone;

/* the iteration counter, set and incremented in routine cloudy,
 * ==1 during first iteration, 2 during second, etc */
extern long int iteration;

/* following three are defined in cddefines */
/* pointers within recombination coefficient array */
/* ipRecEsc is state specific escape probability*/
extern const int ipRecEsc /*= 2*/;
/* the net escape prob, including destruction by background and optical deepth*/
extern const int ipRecNetEsc /*= 1*/;
/* ipRecRad is state specific radiative recombination rate*/
extern const int ipRecRad /*= 0*/;
/* with the above, the total radiative rec per ion is
 * iso.RadRecomb[ipISO][nelem][n][ipRecRad]*
   iso.RadRecomb[ipISO][nelem][n][ipRecNetEsc]*dense.eden; */

/* these specify the form of the line redistribution function */
/* partial redistribution with wings */
extern const int ipPRD; /* = 1 */
/* complete redistribution, core only, no wings, Hummer's K2 function */
extern const int ipCRD; /* = -1 */
/* complete redistribution with wings */
extern const int ipCRDW; /* = 2 */
/* redistribution function for Lya, calls Hummer routine for H-like series only */
extern const int ipLY_A; /* = -1 */
/* core function for K2 destruction */
extern const int ipDEST_K2; /*= 1*/
/* core function for complete redist destruction */
extern const int ipDEST_INCOM; /*= 2*/
/* core function for simple destruction */
extern const int ipDEST_SIMPL; /*= 3*/

/* some levels for hydrogenic species, defined in cddefines.c,
 * where set to 0, 1, and 2  */
extern const int ipH1s;
extern const int ipH2s;
extern const int ipH2p;

/* some levels for he-like species, defined in cddefines.c,
 * where set to 0, 1, 2, 3, and 4  */
/* some levels for he-like species, defined in cddefines.c */

/* level 1 */
extern const int ipHe1s1S;

/* level 2 */
extern const int ipHe2s3S;
extern const int ipHe2s1S;
extern const int ipHe2p3P0;
extern const int ipHe2p3P1;
extern const int ipHe2p3P2;
extern const int ipHe2p1P;

/* level 3 */
extern const int ipHe3s3S;
extern const int ipHe3s1S;
extern const int ipHe3p3P;
extern const int ipHe3d3D;
extern const int ipHe3d1D;
extern const int ipHe3p1P;

/* level 4 */
extern const int ipHe4s3S;
extern const int ipHe4s1S;
extern const int ipHe4p3P;
extern const int ipHe4d3D;
extern const int ipHe4d1D;
extern const int ipHe4f3F;
extern const int ipHe4f1F;
extern const int ipHe4p1P;

/* level 5 */
extern const int ipHe5s3S;
extern const int ipHe5s1S;
extern const int ipHe5p3P;
extern const int ipHe5d3D;
extern const int ipHe5d1D;
extern const int ipHe5f3F;
extern const int ipHe5f1F;
extern const int ipHe5g3G;
extern const int ipHe5g1G;
extern const int ipHe5p1P;

/* these are array indices for isoelectronic sequences,
 * same as element but used for array addressing to make
 * context totally clear */
extern const int ipH_LIKE/*=0*/;
extern const int ipHE_LIKE/*=1*/;
extern const int ipLI_LIKE/*=2*/;
extern const int ipBE_LIKE/*=3*/;
extern const int ipB_LIKE/*=4*/;
extern const int ipC_LIKE/*=5*/;
extern const int ipN_LIKE/*=6*/;
extern const int ipO_LIKE/*=7*/;
extern const int ipF_LIKE/*=8*/;

/* these are incides to some elements, defined in cddefines.c, on the C scale */
extern const int ipHYDROGEN/*=0*/;
extern const int ipHELIUM/*=1*/;
extern const int ipLITHIUM/*=2*/;
extern const int ipBERYLLIUM/*=3*/;
extern const int ipBORON/*=4*/;
extern const int ipCARBON/*=5*/;
extern const int ipNITROGEN/*=6*/;
extern const int ipOXYGEN/*=7*/;
extern const int ipFLUORINE/*=8*/;
extern const int ipNEON/*=9*/;
extern const int ipSODIUM/*=10*/;
extern const int ipMAGNESIUM/*=11*/;
extern const int ipALUMINIUM/*=12*/;
extern const int ipSILICON/*=13*/;
extern const int ipPHOSPHORUS/*=14*/;
extern const int ipSULPHUR/*=15*/;
extern const int ipCHLORINE/*=16*/;
extern const int ipARGON/*=17*/;
extern const int ipPOTASSIUM/*=18*/;
extern const int ipCALCIUM/*=19*/;
extern const int ipSCANDIUM/*=20*/;
extern const int ipTITANIUM/*=21*/;
extern const int ipVANADIUM/*=22*/;
extern const int ipCHROMIUM/*=23*/;
extern const int ipMANGANESE/*=24*/;
extern const int ipIRON/*=25*/;
extern const int ipCOBALT/*=26*/;
extern const int ipNICKEL/*=27*/;
extern const int ipCOPPER/*=28*/;
extern const int ipZINC/*=29*/;

/**************************************************************************
 *
 * these are constants used to dimension several vectors
 *
 **************************************************************************/

/* following is real limit to how many levels can be computed for
 * model hydrogen atom - this has the one extra included, so at most
 * iso.numLevels_max[ipH_LIKE] can be NHYDRO_MAX_LEVEL-1 and vectors should be dim NHYDRO_MAX_LEVEL*/
#define NHYDRO_MAX_LEVEL 401L

/* This is the number of elements included in the code,
 * is used to set lengths of many vectors */
#define	LIMELM	30L

/* dimension of vector to do with number of iterations stored
 * so, this is the the most that can possibly be performed */
#define	ITRDIM	200L

/* number of hydrogen species in hmole */
#define	N_H_MOLEC	8

/* this is used to add to depth to prevent div or log of zero */
#define DEPTH_OFFSET	1E-30
/* safe, small, numbers for the float and double */
/* set something that is too small to max of quantity and SMALLFLOAT,
 * but then compare with SMALLFLOAT */
/*FLT_MAX is 3.40e38 on wintel, so BIGFLOAT is 3.40e36 */
#define BIGFLOAT (FLT_MAX/100.f)
/*FLT_MIN is 1.175494351e-38 on wintel, so SMALLFLOAT is 1.175e-36 */
#define SMALLFLOAT (FLT_MIN*100.f)
/*DBL_MAX is 1.797e308 on wintel so BIGDOUBLE is 1.797e306 */
#define BIGDOUBLE (DBL_MAX/100.)
#define SMALLDOUBLE (DBL_MIN*100.)

/*
 * this is the number zero, used to trick clever compilers when 
 * dividing by it to crash program
 * there is a routine called zero - this name cannot overlap
 * definition is in 
 */
extern const double ZeroNum;

/**************************************************************************
 *
 * disable some bogus errors in the ms c compiler
 *
 **************************************************************************/

/* */
#ifdef _MSC_VER
#	pragma warning( disable : 4127 )/* disable warning that conditional expression*/
                                 /* is constant, TRUE or FALSE in if        */
#	pragma warning( disable : 4056 )/* disable bogus underflow warning in MS VS*/
#	pragma warning( disable : 4514 )/* disable "inline function removed since not used", MS VS*/
#endif
#ifdef __INTEL_COMPILER
#	pragma warning( disable : 1572 )/* disable warning that floating-point comparisons are unreliable */
#endif
/* */

/**************************************************************************
 *
 * various macros used by the code
 *
 **************************************************************************/

/* file ptr for optional C Debug code,
 * entry and exit of each routine will go here,
 * macros enabled if compiler-set flag DEBUG_FUN is true */
#ifdef DEBUG_FUN
extern FILE *debug_fp;
#endif

/* will be redefined as nothing in place where variables are defined */
#define EXTERN extern 

#ifndef TorF
/* TorF(l) returns a 'T' or 'F' depending on the 'logical' expr 'l' */
/* >>chng 01 aug 06, add par, as per pvh comment */
/*#define TorF(l)	( l ? 'T' : 'F' )*/
#define TorF(l) ( (l) ? 'T' : 'F' )
#endif
/* */

#ifndef IS_ODD
/* checks whether argument is odd */
#define IS_ODD(j)	((j) & 1 )
#endif

#ifndef MIN2
/* MIN2 takes two arguments, returns the smaller of the two */
#define MIN2(a,b) (((a)<(b)) ? (a) : (b))
#endif
/* */

#ifndef MIN3
/* MIN3 takes 3 arguments, returns the smallest of the 3 */
#define MIN3(a,b,c) (MIN2(MIN2(a,b),c))
#endif
/* */

#ifndef MIN4
/* MIN2 takes 4 arguments, returns the smallest of the 4 */
#define MIN4(a,b,c,d) (MIN2(MIN2(a,b),MIN2(c,d)))
#endif
/* */

#ifndef MAX2
/* MAX2 takes two arguments, returns the larger of the two */
#define MAX2(a,b) (((a)>(b)) ? (a) : (b))
#endif
/* */

#ifndef SDIV
/* SDIV safe division - div by SDIV(x) - if abs val of arg >SMALLFLOAT,
 * returns arg, if < SMALLFLOAT, returns SMALLFLOAT - with negative arg
 * return is +SMALLFLOAT so sign changes */
#define SDIV(b) (((SMALLFLOAT)>(fabs(b))) ? (SMALLFLOAT) : (b))
#endif
/* */

#ifndef MAX3
/* MAX3 takes 3 arguments, returns the largest of the 3 */
#define MAX3(a,b,c) (MAX2(MAX2(a,b),c))
#endif
/* */

#ifndef MAX4
/* MAX2 takes 4 arguments, returns the largest of the 4 */
#define MAX4(a,b,c,d) (MAX2(MAX2(a,b),MAX2(c,d)))
#endif
/* */

#ifndef POW2
/* POW2 takes 1 argument, and squares it */
#define POW2(a) ((a)*(a))
#endif
/* */

#ifndef POW3
/* POW3 takes 1 argument, and cubes it */
#define POW3(a) ((a)*(a)*(a))
#endif
/* */

/* macro to determine size of an array -
 * from p 22, Practice of Programming, 
 * by B.W. Kernighan & Rob Pike, Addison Wesley.
 * use: 
 * double dbug[100];
 for( i=0; i<NELEMS(dbug); i++)
  ....
 */
#ifdef NELEMS
#undef NELEMS
#endif
#define NELEMS(array) (sizeof(array)/sizeof(array[0]))

#ifndef _ERR_
#define _ERR_	-1
#endif
/* */

/* all vendors supply predefined preprocessor symbols to
 * help identify their hardware/operating system/compiler, the
 * following symbols will be used to bracket hardware/OS/compiler
 * specific code:
 *
 * Sun Sparc:      __sun
 * DEC Alpha:      __alpha
 * SGI Iris:       __sgi
 * HP Unix:        __hpux
 * Cray:           __cray
 * IA32:           __i386
 * AMD64/EM64T:    __amd64
 * UNIX:           __unix (includes Linux)
 * Linux:          __linux
 * MS Vis C:       _MSC_VER
 * Intel compiler: __ICC, __INTEL_COMPILER
 * gcc/g++/icc:    __GNUC__ (also set by icc!)
 * gcc/g++:        __GNUC_EXCL__ (excludes icc)
 *
 * NOTE: the user should NOT define these symbols at compile time.
 */
#ifdef cray
#ifndef __cray
#define __cray 1
#endif
#endif

/* the Intel EM64T compiler does not set the __amd64 flag... */
#ifdef __x86_64
#ifndef __amd64
#define __amd64 1
#endif
#endif

/* on some UNIX systems only the preprocessor symbol "unix"
 * is predefined (e.g. DEC alpha), on others only "__unix"
 * (e.g. Cray), and on yet others both... This ensures
 * "__unix" is allways defined on all UNIX systems.
 */
#if defined(unix) || defined(__unix__)
#ifndef __unix
#define __unix 1
#endif
#endif

/* on SGI IA64 systems, icc calls itself ecc... */
#ifdef __ECC
#ifndef __ICC
#define __ICC __ECC
#endif
#endif

/* this is needed to discriminate between gcc/g++ and icc */
#undef __GNUC_EXCL__
#if defined(__GNUC__) && !defined(__ICC)
#define __GNUC_EXCL__ 1
#endif

/*
 * a series of Cloudy service routines, used throughout code,
 *
 */

/*AnuUnit produce continuum energy in arbitrary units, ip is on C scale */
double AnuUnit(float energy);

/* print comment asking to show output to me */
void ShowMe(void);

/*TotalInsanity general error handler for something that cannot happen, exits */
#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
/*@noreturn@*/
void TotalInsanity(void);

/*BadRead tried to read internal data and failed */
#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
/*@noreturn@*/
void BadRead(void);

/*BadOpen tried to open data files and failed */
#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
/*@noreturn@*/
void BadOpen(void);

/*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*/

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

/*e2 second exponential integral */
double e2(double t, 
	  double tln);

/*ee1 first exponential integral */
double ee1(double x);

/*this one same as ee1, except is divided by a factor of exp(x), and is
 * only to be used for x>1.	*/
double ee1_safe(double x);

/* FFmtRead - thre free format reader */
double FFmtRead(char *chCard, 
	  long int *ipnt, 
	  long int last, 
	  int *lgEOL);

/*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);

/*GetElem scans line image, finds element. returns atomic number j, on C scale */
long int GetElem( char *chCard );

/* 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 );

/*prototype for powi function, which exists in the alpha
 * math library */
/* #ifndef __alpha
 * double powi( double , long int );
 * #endif */
/* want to define this for all cases except alphas with native os */
#if !defined(__alpha) || defined(__linux)
double powi( double , long int );
#endif

/*qip compute pow(x,n) for positive integer n through repeated sqarings */
double qip( double x , unsigned int n);

/* returns:  m^n */
long int ipow( long, long );

/* print with 1p,e8.2 format onto stream FILE 
 * all are located in printe82.c */
void PrintE82( FILE*, double );
/* print with 1p,e8.1 format onto stream FILE */
void PrintE71( FILE*, double );
/* print with 1p,e9.3 format onto stream FILE */
void PrintE93( FILE*, double );
/* create string with val and format, to print with %s,
 * much faster than above, totally native on non-MS systems */
char *PrintEfmt(const char *fmt, double val );

/*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
/*@noreturn@*/
void NoNumb(char *chCard);

 /* this is -ln of smallest number sexp can handle */
#define SEXP_LIMIT (84.)
/*sexp safe exponential function */
double sexp(double x);

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

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

/* floating point sign transfer - signe of b times abs value of a */
double sign( double a, double b );		

/* >>chng 01 mar 08, remove following, and replace with next */
#if 0
#ifndef __HP_CXD_SPP /* MIN and MAX are defined on examplars */ 
#define 	MAX(a,b)	( (a) > (b) ? (a) : (b) )
#endif
#ifndef __HP_CXD_SPP
#define 	MIN(a,b)	( (a) < (b) ? (a) : (b) )
#endif
#endif

#ifndef MAX /* MIN and MAX are defined on some machines */ 
#define 	MAX(a,b)	( (a) > (b) ? (a) : (b) )
#endif
#ifndef MIN
#define 	MIN(a,b)	( (a) < (b) ? (a) : (b) )
#endif

/*long max(long,long);*/
/* Variable no. of arg definitions for max & min functs */
double vfmax( double,... );	/* FLOATING POINT MAXIMUM */
/*long vmax( long, ... );*/ 	/* LONG INTEGER MAXIMUM */
/* short vsmax(short,...);		short INTEGER MAXIMUM */
double vfmin( double,... );	/* FLOATING POINT MINIMUM */
/*long vmin( long, ... );*/ 	/* LONG INTEGER MINIMUM */
/* short vsmin(short,...);		short INTEGER MINIMUM */

/* Variable argument list terminators (used by max & min functs) */
#ifndef INT_MAX
#  include <limits.h>
#endif	/* INT_MAX */
#define IEND	(LONG_MAX-1)
#define FEND	(1.e+38-1.)

#undef HMRATE
/* >>chng 04 aug 14, PvH chng HMRATE into macro to speed up co.c */
/*HMRATE compile molecular rates using Hollenbach and McKee fits */
#define HMRATE(a,b,c) ( ((b) == 0 && (c) == 0) ? (a) : \
	( ((c) == 0) ? (a)*pow(phycon.te/300.,(b)) : \
	( ((c)/phycon.te > 50.) ? 0. : ( ((b) == 0) ?  (a)*exp(-(c)/phycon.te) : \
					 (a)*pow(phycon.te/300.,(b))*exp(-(c)/phycon.te) ) ) ) )


/* 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, ...);

double qg32(			/* 32 point gaussian quadrature integration */
				double,	/* lower limit to integration */
				double,	/* upper limit to integration */
				double(*)(double)	/* pointer to routine to be integrated, arg is x val */
				);
/* declar of optimize_func, the last arg, changed from double(*)() to above,
 * seemed to fix flags that were raised */


/*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 routine */
	  int kflag, 
	  /* error condition, should be 0 */
	  int *ier);
