/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */

#ifndef CDDEFINES_H_
#define CDDEFINES_H_

#include "cdstd.h"

#ifdef _MSC_VER
/* disable warning that conditional expression is constant, true or false in if */
#	pragma warning( disable : 4127 )
/* we are not using MS foundation class */
#	ifndef WIN32_LEAN_AND_MEAN
#		define WIN32_LEAN_AND_MEAN
#	endif
#endif

#ifdef __clang__
// this would generate lots of warnings about mismatched tags in the STL valarray definition
#pragma clang diagnostic ignored "-Wmismatched-tags"
#endif

/* these headers are needed by all files */
/*lint -e129 these resolve several issues pclint has with my system headers */
/*lint -e78 */
/*lint -e830 */
/*lint -e38 */
/*lint -e148 */
/*lint -e114 */
/*lint -e18 */
/*lint -e49 */
// C++ versions of C headers
#include <cstdio>
#include <cstdlib>
#include <cctype>
#ifdef _MSC_VER
// MSVC needs this before cmath in order provide numeric constants 
// (M_PI etc.) defined by C99 but not C++ standards to date.
#define _USE_MATH_DEFINES 
#endif
#include <cmath>
#include <cassert>
#include <cstring>
#include <cfloat>
#include <climits>
#include <ctime>
#if defined(__sun) && defined(__SUNPRO_CC)
// with Solaris Studio 12.2 under Sparc Solaris, csignal doesn't define sigaction...
#include <signal.h>
#else
#include <csignal>
#endif
// C++ headers
#include <limits>
#include <string>
#include <sstream>
#include <iomanip>
#include <vector>
#include <valarray>
#include <complex>
#include <map>
#include <memory>
#include <stdexcept>
#include <algorithm>
#include <fstream>
#include <bitset>
#ifdef DMALLOC
#include <dmalloc.h>
#endif

// Workaround for Windows...
#if defined(_MSC_VER) && !defined(SYS_CONFIG)
#define SYS_CONFIG "cloudyconfig_vs.h"
#endif

// platform specific configuration; generated by configure.sh
#ifdef SYS_CONFIG
#include SYS_CONFIG
#else
#include "cloudyconfig.h"
#endif

#ifdef MPI_GRID_RUN
#define MPI_ENABLED
#endif

/*lint +e18 */
/*lint +e49 */
/*lint +e38 */
/*lint +e148 */
/*lint +e830 */
/*lint +e78 */
/*lint -e129 */

using namespace std;

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

#ifdef FLT_IS_DBL
typedef double realnum;
#else
typedef float realnum;
#endif

typedef float sys_float;
// prevent explicit float's from creeping back into the code
#define float PLEASE_USE_REALNUM_NOT_FLOAT

//Compile-time assertion after Alexandrescu
template<bool> struct StaticAssertFailed;
template<> struct StaticAssertFailed<true> {};
#define STATIC_ASSERT(x) ((void)StaticAssertFailed< (x) == true >())

typedef enum {
	ES_SUCCESS=0,            // everything went fine...
	ES_FAILURE=1,            // general failure exit
	ES_WARNINGS,             // warnings were present
	ES_BOTCHES,              // botched monitors were present
	ES_CLOUDY_ABORT,         // Cloudy aborted
	ES_BAD_ASSERT,           // an assert in the code failed
	ES_BAD_ALLOC,            // a memory allocation failed
	ES_OUT_OF_RANGE,         // an out-of-range exception was thrown
	ES_USER_INTERRUPT,       // the user terminated Cloudy (with ^C)
	ES_TERMINATION_REQUEST,  // Cloudy received a termination request
	ES_ILLEGAL_INSTRUCTION,  // the CPU encountered an illegal instruction
	ES_FP_EXCEPTION,         // a floating point exception was caught
	ES_SEGFAULT,             // a segmentation fault occurred
	ES_BUS_ERROR,            // a bus error occurred
	ES_UNKNOWN_SIGNAL,       // an unknown signal was caught
	ES_UNKNOWN_EXCEPTION,    // an unknown exception was caught
	ES_TOP                   // NB NB -- this should always be the last entry
} exit_type;

// make sure the system definitions are on par with ours
// especially EXIT_FAILURE does not have a guaranteed value!
#undef EXIT_SUCCESS
#define EXIT_SUCCESS ES_SUCCESS
#undef EXIT_FAILURE
#define EXIT_FAILURE ES_FAILURE

/* make sure this is globally visible as well! */
/* This must be done at the start of every file, to ensure that policy
	for FPE handling, etc., is guaranteed to be set up before the
	construction of file-statics and globals. */
#include "cpu.h"

//*************************************************************************
//
//! Singleton: template to construct classes where a single instance can
//! be accessed across the entire program. Construct your class as follows:
//!
//! class t_myclass : public Singleton<t_myclass>
//! {
//! 	friend class Singleton<t_myclass>;
//! protected:
//! 	t_myclass(); // make sure the contructor is protected!
//! public:
//! 	long myfunc() const { return 43; }
//! };
//!
//! and use as follows throughout the code:
//!
//! long l = t_myclass::Inst().myfunc();
//!
//!   NB NB - This implementation is not threadsafe !!
//
// This implementation has been obtained from Wikipedia
//
//*************************************************************************

template<typename T> class Singleton
{
public:
	static T& Inst()
	{
		static T instance;  // assumes T has a protected default constructor
		return instance;
	}
};

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

/** ioQQQ is the file handle to the output file itself,
 * ioQQQ is set to stdout by default, 
 * and is reset to anything else by calling cdOutput */
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 set true when abort sequence is initiated - serious meltdown is happening */
extern bool lgAbort;

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

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

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

/** nzone is zone counter, incremented in routine cloudy 
 * is zero during search phase, 1 for first zone at illuminated face */
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;

/**
 * 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 cddefines.cpp
 */
extern const double ZeroNum;

/**************************************************************************
 *
 * these are constants used to dimension several vectors and index arrays
 *
 **************************************************************************/

/** 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 sentinel in C.  Increase this is a larger string
 * is needed to hold the path on your system */
const int FILENAME_PATH_LENGTH = 200; 

/** twice the above, so that we can add file name to end of full path */
const int 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 */ 
const int INPUT_LINE_LENGTH = 2000;

/** This is the number of elements included in the code,
 * is used to set lengths of many vectors */
const int LIMELM = 30;

/** the number of iso sequences now in the code */
const int NISO = 2;

/** 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*/
const int NHYDRO_MAX_LEVEL = 401;

/** this is the maximum particle density allowed in cm^-3 */
const double MAX_DENSITY = 1.e24;

/** this is used to add to depth to prevent div or log of zero */
const double DEPTH_OFFSET = 1.e-30;

enum {CHARS_SPECIES=10};
enum {CHARS_ISOTOPE_SYM	= 6};

/* indices within recombination coefficient array */
/* ipRecEsc is state specific escape probability*/
const int ipRecEsc = 2;
/* the net escaping, including destruction by background and optical deepth*/
const int ipRecNetEsc = 1;
/* ipRecRad is state specific radiative recombination rate*/
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 */
const int ipPRD = 1;
/* complete redistribution, core only, no wings, Hummer's K2 function */
const int ipCRD = -1;
/* complete redistribution with wings */
const int ipCRDW = 2;
/* redistribution function for Lya, calls Hummer routine for H-like series only */
const int ipLY_A = -2;
/* core function for K2 destruction */
const int ipDEST_K2 = 1;
/* core function for complete redist destruction */
const int ipDEST_INCOM = 2;
/* core function for simple destruction */
const int ipDEST_SIMPL = 3;

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

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

/**
 fudge enter fudge factors, or some arbitrary number, with fudge command
 return value is the fudge factor 
 fudge(-1) queries the routine for the number of fudge parameters that 
 were entered, zero returned if none
 \param ipnt integer saying which of the possible numbers on the fudge
 command to use - 0 would be the first
 */ 
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 service.cpp */
void fixit(void);

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

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

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

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

/**MyRealloc wrapper for realloc().  Returns a good pointer or dies. 
\param num use same type as library function REALLOC
\param size
*/
void *MyRealloc(void *p, size_t size);

/** MyAssert a version of assert that fails gracefully
\param *file
\param line
*/ 
void MyAssert(const char *file, int line, const char *comment);

/** prepare termination of the code, but do not terminate yet */
void cdPrepareExit(exit_type);

class cloudy_exit
{
	const char* p_routine;
	const char* p_file;
	long p_line;
	exit_type p_exit;
public:
	cloudy_exit(const char* routine, const char* file, long line, exit_type exit_code)
	{
		p_routine = routine;
		p_file = file;
		p_line = line;
		p_exit = exit_code;
	}
	virtual ~cloudy_exit() throw()
	{
		p_routine = NULL;
		p_file = NULL;
	}
	const char* routine() const throw()
	{
		return p_routine;
	}
	const char* file() const throw()
	{
		return p_file;
	}
	long line() const
	{
		return p_line;
	}
	exit_type exit_status() const
	{
		return p_exit;
	}
};

// workarounds for __func__ are defined in cpu.h
#define cdEXIT( FAIL ) throw cloudy_exit( __func__, __FILE__, __LINE__, FAIL )

// calls like puts( "[Stop in MyRoutine]" ) have been integrated in cdEXIT above
#define puts( STR ) Using_puts_before_cdEXIT_is_no_longer_needed

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

/**TotalInsanity general error handler for something that cannot happen, exits */
NORETURN void TotalInsanity(void);

/* TotalInsanityAsStub always calls TotalInsanity(), but in such a way that
 * it can be used as a stub for another routine without generating warnings
 * about unreachable code after the stub. Hence this should NOT be NORETURN */
template<class T>
T TotalInsanityAsStub()
{
	// this is always true...
	if( ZeroNum == 0. )
		TotalInsanity();
	else
		return T();
}

/**BadRead tried to read internal data and failed */
NORETURN void BadRead(void);

/** 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 */ 
int dbg_printf(int debug, const char *fmt, ...);

/** dprintf -- version of fprintf which prepends DEBUG */
int dprintf(FILE *fp, const char *format, ...);

/**read_whole_line safe version of fgets - read a line, 
 * return null if cannot read line or if input line is too long 
 \param char *chLine - previously allocated string where the line
 image will be stored
 \param int nChar size of chLine, we will return NULL if input line is
 longer than this
 \param FILE *ioIN a previously opened file handle, will read from from here
*/
char *read_whole_line( char *chLine , int nChar , FILE *ioIN );

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

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

/** use special version of malloc - it tests result and dies if cannot allocate space */
#if defined(malloc)
/* ...but if malloc is a macro, assume it is instrumented by a memory debugging tool
 * (e.g. dmalloc) */
#	define MALLOC(exp) (malloc(exp))
#else
/* Otherwise instrument and protect it ourselves */
#	define MALLOC(exp) (MyMalloc(exp,__FILE__, __LINE__))
#endif

/** now special version of calloc - it dies if cannot allocate space. */
#if defined(calloc)
/* ...but if calloc is a macro, assume it is instrumented by a memory debugging tool */
#	define CALLOC calloc
#else
/* Otherwise instrument and protect it ourselves */
#	define CALLOC MyCalloc 
#endif

/** now special version of calloc - it dies if cannot allocate space. */
#if defined(realloc)
/* ...but if calloc is a macro, assume it is instrumented by a memory debugging tool */
#	define REALLOC realloc
#else
/* Otherwise instrument and protect it ourselves */
#	define REALLOC MyRealloc 
#endif

class bad_signal
{
	int p_sig;
public:
	explicit bad_signal(int sig)
	{
		p_sig = sig;
	}
	virtual ~bad_signal() throw() {}
	int sig() const throw()
	{
		return p_sig;
	}
};

class bad_assert
{
	const char* p_file;
	long p_line;
	const char* p_comment;
public:
	bad_assert(const char* file, long line, const char* comment);
	void print(void) const
	{
		fprintf(ioQQQ,"DISASTER Assertion failure at %s:%ld\n%s\n",
				  p_file, p_line, p_comment);
	}
	virtual ~bad_assert() throw()
	{
		p_file = NULL;
	}
	const char* file() const throw()
	{
		return p_file;
	}
	long line() const throw()
	{
		return p_line;
	}
	const char *comment() const throw()
	{
		return p_comment;
	}
};

/* the do { ... } while ( 0 ) construct prevents bugs in code like this:
 * if( test )
 * 	ASSERT( n == 10 );
 * else
 *	do something else...
 */
#undef  ASSERT
#ifndef OLD_ASSERT
#	if NDEBUG
#		define ASSERT(exp) ((void)0)
#	else
#		define ASSERT(exp) \
			do { \
				if (UNLIKELY(!(exp)))					\
				{ \
					bad_assert aa(__FILE__,__LINE__,"Failed: " #exp); \
					if( cpu.i().lgAssertAbort() ) \
					{ \
						aa.print();	\
						abort();	\
					}  \
					else \
						throw aa; \
				} \
			} while( 0 )
#	endif
#else
/** the old version of the assert macro that terminates safely */
#	ifdef NDEBUG
#		define ASSERT(exp) ((void)0)
#	else
#		define ASSERT(exp) \
			do { \
				if (!(exp)) \
					MyAssert(__FILE__, __LINE__, "Failed: " #exp); \
			} while( 0 )
#	endif
#endif

#define MESSAGE_ASSERT(msg, exp) ASSERT( (msg) ? (exp) : false )

inline NORETURN void OUT_OF_RANGE(const char* str)
{
	if( cpu.i().lgAssertAbort() )
		abort();
	else
		throw out_of_range( str );
}

/* Windows does not define isnan */
/* use our version on all platforms since the isnanf
 * function does not exist under Solaris 9 either */
#undef isnan
#define isnan MyIsnan

/** entry and exit of each routine will go here,
 * macros enabled if compiler-set flag DEBUG_FUN is defined */
class t_debug : public Singleton<t_debug>
{
	friend class Singleton<t_debug>;
	FILE *p_fp;
	int p_callLevel;
protected:
	t_debug() : p_fp(stderr)
	{
		p_callLevel = 0;
	}
public:
	void enter(const char *name)
	{
		++p_callLevel;
		fprintf(p_fp,"%*c%s\n",p_callLevel,'>',name);
	}
	void leave(const char *name)
	{
		fprintf(p_fp,"%*c%s\n",p_callLevel,'<',name);
		--p_callLevel;
	}		
};

/** entry and exit of each routine will go here,
 * stub routines for when compiler-set flag DEBUG_FUN is not defined */
class t_nodebug : public Singleton<t_nodebug>
{
	friend class Singleton<t_nodebug>;
protected:
	t_nodebug()	{}
public:
	void enter(const char *) const {}
	void leave(const char *) const {}		
};

template<class Trace>
class debugtrace
{
	const char *p_name;
public:
	explicit debugtrace(const char *funcname)
	{
		p_name = funcname;
		Trace::Inst().enter(p_name);
	}
	~debugtrace()
	{
		Trace::Inst().leave(p_name);
		p_name = NULL;
	}
	const char* name() const
	{
		return p_name;
	}
};

#ifdef DEBUG_FUN
#define DEBUG_ENTRY( funcname ) debugtrace<t_debug> DEBUG_ENTRY( funcname )
#else
#ifdef HAVE_FUNC
#define DEBUG_ENTRY( funcname ) ((void)0)
#else
#define DEBUG_ENTRY( funcname ) debugtrace<t_nodebug> DEBUG_ENTRY( funcname )
#endif
#endif

// overload the character manipulation routines
inline char tolower(char c)
{
	return static_cast<char>( tolower( static_cast<int>(c) ) );
}
inline unsigned char tolower(unsigned char c)
{
	return static_cast<unsigned char>( tolower( static_cast<int>(c) ) );
}

inline char toupper(char c)
{
	return static_cast<char>( toupper( static_cast<int>(c) ) );
}
inline unsigned char toupper(unsigned char c)
{
	return static_cast<unsigned char>( toupper( static_cast<int>(c) ) );
}

/* TorF(l) returns a 'T' or 'F' depending on the 'logical' expr 'l' */
inline char TorF( bool l ) { return l ? 'T' : 'F'; }
/* */

/** checks whether argument is odd */
inline bool is_odd( int j ) { return (j&1) == 1; }
inline bool is_odd( long j ) { return (j&1L) == 1L; }
/* */

/** nint rounds to the nearest long int */
inline long nint( double x ) { return static_cast<long>( (x < 0.) ? x-0.5 : x+0.5 ); }
/* */

/* define min for mixed arguments, the rest already exists */
inline long min( int a, long b ) { long c = a; return ( (c < b) ? c : b ); }
inline long min( long a, int b ) { long c = b; return ( (a < c) ? a : c ); }
inline double min( sys_float a, double b ) { double c = a; return ( (c < b) ? c : b ); }
inline double min( double a, sys_float b ) { double c = b; return ( (a < c) ? a : c ); }

/* want to define this only if no native os support exists */
#ifndef HAVE_POWI
/** powi raise x to integer power */
double powi( double , long int );
#endif

/* avoid ambiguous overloads */
#ifndef HAVE_POW_DOUBLE_INT
inline double pow( double x, int i ) { return powi( x, long(i) ); }
#endif

#ifndef HAVE_POW_DOUBLE_LONG
inline double pow( double x, long i ) { return powi( x, i ); }
#endif

#ifndef HAVE_POW_FLOAT_INT
inline sys_float pow( sys_float x, int i ) { return sys_float( powi( double(x), long(i) ) ); }
#endif

#ifndef HAVE_POW_FLOAT_LONG
inline sys_float pow( sys_float x, long i ) { return sys_float( powi( double(x), i ) ); }
#endif

#ifndef HAVE_POW_FLOAT_DOUBLE
inline double pow( sys_float x, double y ) { return pow( double(x), y ); }
#endif

#ifndef HAVE_POW_DOUBLE_FLOAT
inline double pow( double x, sys_float y ) { return pow( x, double(y) ); }
#endif

#undef MIN2
/** MIN2 takes two arguments, returns the smaller of the two */
#define MIN2 min
/* */

#undef MIN3
/** MIN3 takes 3 arguments, returns the smallest of the 3 */
#define MIN3(a,b,c) (min(min(a,b),c))
/* */

#undef MIN4
/** MIN4 takes 4 arguments, returns the smallest of the 4 */
#define MIN4(a,b,c,d) (min(min(a,b),min(c,d)))
/* */

/* define max for mixed arguments, the rest already exists */
inline long max( int a, long b ) { long c = a; return ( (c > b) ? c : b ); }
inline long max( long a, int b ) { long c = b; return ( (a > c) ? a : c ); }
inline double max( sys_float a, double b ) { double c = a; return ( (c > b) ? c : b ); }
inline double max( double a, sys_float b ) { double c = b; return ( (a > c) ? a : c ); }

#undef MAX2
/** MAX2 takes two arguments, returns the larger of the two */
#define MAX2 max
/* */

#undef MAX3
/** MAX3 takes 3 arguments, returns the largest of the 3 */
#define MAX3(a,b,c) (max(max(a,b),c))
/* */

#undef MAX4
/** MAX4 takes 4 arguments, returns the largest of the 4 */
#define MAX4(a,b,c,d) (max(max(a,b),max(c,d)))
/* */

/** FP sign transfer (fortran sign function) - sign of y times abs value of x 
\param x
\param y
*/
template<class T>
inline T sign( T x, T y )
{
	return ( y < T() ) ? -abs(x) : abs(x);
}
/* */

/** sign3 returns -1 for negative arguments, +1 for positive, and 0 for zero (pascal sign function) */
template<class T>
inline int sign3( T x ) { return ( x < T() ) ? -1 : ( ( x > T() ) ? 1 : 0 ); }
/* */

/** checks whether two FP numbers are "equal" (differ no more than n epsilon) */
inline bool fp_equal( sys_float x, sys_float y, int n=3 )
{
#ifdef _MSC_VER
	/* disable warning that conditional expression is constant, true or false in if */
#	pragma warning( disable : 4127 )
#endif
	ASSERT( n >= 1 );
	// mimic IEEE behavior
	if( isnan(x) || isnan(y) )
		return false;
	int sx = sign3(x);
	int sy = sign3(y);
	// treat zero cases first to avoid division by zero below
	if( sx == 0 && sy == 0 )
		return true;
	// either x or y is zero (but not both), or x and y have different sign
	if( sx*sy != 1 )
		return false;
	x = abs(x);
	y = abs(y);
	return ( 1.f - min(x,y)/max(x,y) < ((sys_float)n+0.1f)*FLT_EPSILON );
}

inline bool fp_equal( double x, double y, int n=3 )
{
	ASSERT( n >= 1 );
	// mimic IEEE behavior
	if( isnan(x) || isnan(y) )
		return false;
	int sx = sign3(x);
	int sy = sign3(y);
	// treat zero cases first to avoid division by zero below
	if( sx == 0 && sy == 0 )
		return true;
	// either x or y is zero (but not both), or x and y have different sign
	if( sx*sy != 1 )
		return false;
	x = abs(x);
	y = abs(y);
	return ( 1. - min(x,y)/max(x,y) < ((double)n+0.1)*DBL_EPSILON );
}

inline bool fp_equal_tol( sys_float x, sys_float y, sys_float tol )
{
	ASSERT( tol > 0.f );
	// mimic IEEE behavior
	if( isnan(x) || isnan(y) )
		return false;
	// make sure the tolerance is not too stringent
	ASSERT( tol >= FLT_EPSILON*max(abs(x),abs(y)) );
	return ( abs( x-y ) <= tol );
}

inline bool fp_equal_tol( double x, double y, double tol )
{
	ASSERT( tol > 0. );
	// mimic IEEE behavior
	if( isnan(x) || isnan(y) )
		return false;
	// make sure the tolerance is not too stringent
	ASSERT( tol >= DBL_EPSILON*max(abs(x),abs(y)) );
	return ( abs( x-y ) <= tol );
}

/** checks whether a number is within bounds */
inline bool fp_bound( sys_float lo, sys_float x, sys_float hi, int n=3 )
{
	ASSERT( n >= 1 );
	// mimic IEEE behavior
	if( isnan(x) || isnan(lo) || isnan(hi) )
		return false;
	if( fp_equal(lo,hi,n) )
		return fp_equal(0.5f*(lo+hi),x,n);
	if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -((sys_float)n+0.1f)*FLT_EPSILON )
		return false;
	return true;
}
inline bool fp_bound( double lo, double x, double hi, int n=3 )
{
	ASSERT( n >= 1 );
	// mimic IEEE behavior
	if( isnan(x) || isnan(lo) || isnan(hi) )
		return false;
	if( fp_equal(lo,hi,n) )
		return fp_equal(0.5*(lo+hi),x,n);
	if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -((double)n+0.1)*DBL_EPSILON )
		return false;
	return true;
}
inline bool fp_bound_tol( sys_float lo, sys_float x, sys_float hi, sys_float tol )
{
	ASSERT( tol > 0.f );
	// mimic IEEE behavior
	if( isnan(x) || isnan(lo) || isnan(hi) )
		return false;
	if( fp_equal_tol(lo,hi,tol) )
		return fp_equal_tol(0.5f*(lo+hi),x,tol);
	if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -tol )
		return false;
	return true;
}
inline bool fp_bound_tol( double lo, double x, double hi, double tol )
{
	ASSERT( tol > 0. );
	// mimic IEEE behavior
	if( isnan(x) || isnan(lo) || isnan(hi) )
		return false;
	if( fp_equal_tol(lo,hi,tol) )
		return fp_equal_tol(0.5*(lo+hi),x,tol);
	if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -tol )
		return false;
	return true;
}


#undef POW2
/** POW2 takes 1 argument, and squares it */
#define POW2 pow2
template<class T>
inline T pow2(T a) { return a*a; }
/* */

#undef POW3
/** POW3 takes 1 argument, and cubes it */
#define POW3 pow3
template<class T>
inline T pow3(T a) { return a*a*a; }
/* */

#undef POW4
/** POW4 takes 1 argument, and raises it to the power 4 */
#define POW4 pow4
template<class T>
inline T pow4(T a) { T b = a*a; return b*b; }
/* */

#undef 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 */
inline sys_float SDIV( sys_float x ) { return ( fabs((double)x) < (double)SMALLFLOAT ) ? (sys_float)SMALLFLOAT : x; }
/* \todo should we use SMALLDOUBLE here ? it produces overflows now... PvH */
inline double SDIV( double x ) { return ( fabs(x) < (double)SMALLFLOAT ) ? (double)SMALLFLOAT : x; }
// inline double SDIV( double x ) { return ( fabs(x) < SMALLDOUBLE ) ? SMALLDOUBLE : x; }
/* */

/** safe_div( x, y ) - do a really safe division x/y
 * returns +/-FLT_MAX if the division would have overflowed (includes div by 0)
 * returns NaN (i.e. crashes) when x or y are NaN, returns res_0by0 when evaluating 0/0 */
inline sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
{
	// this should crash...
	if( isnan(x) || isnan(y) )
		return x/y;
	int sx = sign3(x);
	int sy = sign3(y);
	// 0/0 -> NaN, this should crash as well...
	if( sx == 0 && sy == 0 )
	{
		if( isnan(res_0by0) )
			return x/y;
		else
			return res_0by0;
	}
	if( sx == 0 )
		return 0.;
	if( sy == 0 )
		return ( sx < 0 ) ? -FLT_MAX : FLT_MAX;
	// at this stage x != 0. and y != 0.
	sys_float ay = abs(y);
	if( ay >= 1.f )
		return x/y;
	else
	{
		// multiplication is safe since ay < 1.
		if( abs(x) < ay*FLT_MAX )
			return x/y;
		else
			return ( sx*sy < 0 ) ? -FLT_MAX : FLT_MAX;
	}
}

inline sys_float safe_div(sys_float x, sys_float y)
{
	return safe_div( x, y, numeric_limits<sys_float>::quiet_NaN() );
}

/** safe_div( x, y ) - do a really safe division x/y
 * returns +/-DBL_MAX if the division would have overflowed (includes div by 0)
 * returns NaN (i.e. crashes) when x or y are NaN, returns res_0by0 when evaluating 0/0 */
inline double safe_div(double x, double y, double res_0by0)
{
	// this should crash...
	if( isnan(x) || isnan(y) )
		return x/y;
	int sx = sign3(x);
	int sy = sign3(y);
	// 0/0 -> NaN, this should crash as well...
	if( sx == 0 && sy == 0 )
	{
		if( isnan(res_0by0) )
			return x/y;
		else
			return res_0by0;
	}
	if( sx == 0 )
		return 0.;
	if( sy == 0 )
		return ( sx < 0 ) ? -DBL_MAX : DBL_MAX;
	// at this stage x != 0. and y != 0.
	double ay = abs(y);
	if( ay >= 1. )
		return x/y;
	else
	{
		// multiplication is safe since ay < 1.
		if( abs(x) < ay*DBL_MAX )
			return x/y;
		else
			return ( sx*sy < 0 ) ? -DBL_MAX : DBL_MAX;
	}
}

inline double safe_div(double x, double y)
{
	return safe_div( x, y, numeric_limits<double>::quiet_NaN() );
}

#undef HMRATE
/*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) ) ) ) ) */
#define HMRATE(a,b,c) hmrate4(a,b,c,phycon.te)

inline double hmrate4( double a, double b, double c, double te )
{
	if( b == 0. && c == 0. )
		return a;
	else if( c == 0. )
		return a*pow(te/300.,b);
	else if( b == 0. )
		return ( c/te <= 50. ) ? a*exp(-c/te) : 0.;
	else
		return ( c/te <= 50. ) ? a*pow(te/300.,b)*exp(-c/te) : 0.;
}

template<class T>
inline void invalidate_array(T* p, size_t size)
{
	if( size > 0 )
		memset( p, -1, size );
}

inline void invalidate_array(double* p, size_t size)
{
	set_NaN( p, (long)(size/sizeof(double)) );
}

inline void invalidate_array(sys_float* p, size_t size)
{
	set_NaN( p, (long)(size/sizeof(sys_float)) );
}

/**get_ptr attribute shim to get raw pointer to contained data with
 * correct type */
template<class T> inline T* get_ptr(T *v)
{
	return v;
}
template<class T> inline T* get_ptr(valarray<T> &v)
{
	return &v[0];
}
template<class T> inline T* get_ptr(vector<T> &v)
{
	return &v[0];
}
template<class T> inline const T* get_ptr(const valarray<T> &v)
{
	return const_cast<const T*>(&const_cast<valarray<T>&>(v)[0]);
}
template<class T> inline const T* get_ptr(const vector<T> &v)
{	
	return const_cast<const T*>(&const_cast<vector<T>&>(v)[0]);
}

/** auto_vec: a smart pointer to allocated arrays which maintains strict ownership.
 * the implementation here largely follows the ISO/ANSI C++98 standard for auto_ptr
 *
 * NB1 - the operator[] and data() methods were added and are not part of the standard
 *
 * NB2 - operator* and operator-> have been omitted as they are not useful for vectors
 *
 * NB3 - the class has been altered to disallow assigning pointers of a different type,
 *       i.e., auto_vec<BaseClass> p( new DerivedClass[10] ) is NOT allowed. This is
 *       done because BaseClass and DerivedClass often have different sizes, and hence
 *       using them like this in vectors is error-prone.
 *
 * NB4 - operator auto_vec<T>() has been omitted because it is only needed to convert
 *       an auto_vec<DerivedClass> to auto_vec<BaseClass>.
 *
 * This is written in the standards draft (with auto_ptr replaced by auto_vec):
 *
 * 1    Template auto_vec holds a pointer to an array object obtained via
 *      new and deletes that object when it itself is destroyed (such
 *      as when leaving block scope 6.7).
 *
 * 2    Template auto_vec_ref holds a reference to an auto_vec. It is
 *      used by the auto_vec conversions to allow auto_vec objects to
 *      be passed to and returned from functions.
 */
template<class T>
class auto_vec
{
	T* ptr;

	template<class U>
	struct auto_vec_ref
	{
		U* ptr;

		explicit auto_vec_ref( U* p )
		{
			ptr = p;
		}
	};

public:
	typedef T element_type;

	// 20.4.5.1 construct/copy/destroy:

	explicit auto_vec( element_type* p = NULL ) throw()
	{
		ptr = p;
	}
	auto_vec( auto_vec& p ) throw()
	{
		ptr = p.release();
	}
	auto_vec& operator= ( auto_vec& p ) throw()
	{
		reset( p.release() );
		return *this;
	}
	~auto_vec() throw()
	{
		delete[] ptr;
	}

	// 20.4.5.2 members:

	element_type& operator[] ( ptrdiff_t n ) const throw()
	{
		return *(ptr+n);
	}
	element_type* get() const throw()
	{
		return ptr;
	}
	// for consistency with other container classes
	element_type* data() const throw()
	{
		return ptr;
	}
	element_type* release() throw()
	{
		element_type* p = ptr;
		ptr = NULL;
		return p;
	}
	void reset( element_type* p = NULL ) throw()
	{
		if( p != ptr )
		{
			delete[] ptr;
			ptr = p;
		}
	}

	// 20.4.5.3 conversions:

	auto_vec( auto_vec_ref<element_type> r ) throw()
	{
		ptr = r.ptr;
	}
	auto_vec& operator= ( auto_vec_ref<element_type> r ) throw()
	{
		if( r.ptr != ptr )
		{
			delete[] ptr;
			ptr = r.ptr;
		}
		return *this;
	}
	operator auto_vec_ref<element_type>() throw()
	{
		return auto_vec_ref<element_type>( this->release() );
	}
};

#include "container_classes.h"
#include "iter_track.h"

/*Many structure were introduced by Humeshkar B Nemala as a part of his Thesis
 *The structures were designed to read in transition,radiative and collisional data
 *from two major databases:LEIDEN and CHIANTI

 * these structures define the emission, collision, state, and transition classes*/

typedef struct t_species species;

#include "lines_service.h"

/*The species structure is used to hold information about a particular atom,ion or molecule
mentioned in the species.ini file.The name of the atom/ion/molecule is used to obtain the density
of molecules in the case of the Leiden Database and along with atomic number and ion stage the 
density of atoms/ions in the case of the CHIANTI database */
struct t_species
{
	/*Name of the atom/ion/ molecule*/
	char *chLabel;
	// index in chmeistry
	long index;
	/*Actual Number of energy levels in the data file*/
	long numLevels_max;
	/*Number of energy levels used locally*/
	long numLevels_local;
	/*Molecular weight*/
	realnum fmolweight;
	/* is molecular? */
	bool lgMolecular;
	// intrepret data as LAMDA or CHIANTI?
	bool lgLAMDA;
	/* fraction in this "type" (e.g. para, ortho) */
	realnum fracType;
	/* chemical fractionation */
	realnum fracIsotopologue;
	/** total cooling due to this species */
	double CoolTotal;
	/** are populations currently defined? */
	bool lgActive;
	/* maximum wavenumber in chianti */
	double maxWN;
	/** should level pops be done in LTE? */
	bool lgLTE;
};

/*This structure is specifically used to hold the collision data in the format given in the LEIDEN Database
The data is available as collision rate coefficients(cm3 s-1) over different temperatures*/
typedef struct t_CollRatesArray
{
	/*Array of temps*/
	vector<double> temps;
	/*Matrix of collision rates(temp,up,lo)*/
	multi_arr<double,3> collrates;
	
} CollRateCoeffArray ;

/*This structure is specifically used to hold the collision data in the format given in the CHIANTI Database
The data is available as spline fits to the Maxwellian averaged collision strengths */
typedef struct t_CollSplinesArray
{
	/*Matrix of spline fits(hi,lo,spline index)*
	 *The first five columns gives the no of spline pts,transition type,gf value,delta E
	 *& Scaling parameter ,in the specified order*/
	/*The transition type basically tells how the temperature and collision
	strengths have been scaled*/
	double *collspline;
	double *SplineSecDer;

	long nSplinePts; 
	long intTranType;
	double EnergyDiff;
	double ScalingParam;

} CollSplinesArray ;

/*This structure is specifically used to hold the collision data in the format given in the STOUT Database
The data are available as collision strengths and rates over different temperatures*/
typedef struct t_StoutColls
{
	/*Number of temps*/
	long ntemps;
	/*Array of temps*/
	double *temps;
	/*Array of collision strengths*/
	double *collstrs;
	/*Is this a deexcitation rate or collision strength*/
	bool lgIsRate;

} StoutColls ;

#include "physconst.h"

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

/** split_mode defines how the routine Split generates substrings
 * SPM_RELAX: multiple adjacent separators will be coalesced into one
 *            this way you can never get an empty substring
 * SPM_KEEP_EMPTY: multiple adjacent separators will result in empty
 *                 substrings to be added to the list
 * SPM_STRICT: empty substrings are illegal */
typedef enum { SPM_RELAX, SPM_KEEP_EMPTY, SPM_STRICT } split_mode;

/** Split: split a string into substrings using "sep" as separator */
void Split(const string& str,   // input string
	   const string& sep,   // separator, may be multiple characters
	   vector<string>& lst, // the separated items will be appended here
	   split_mode mode);    // see above

/** in string str, replace the first instance of substr with newstr
 *  returns true if an instance of substr was found and replaced */
inline bool FindAndReplace(string& str,
			   const string& substr,
			   const string& newstr)
{
	string::size_type ptr = str.find( substr );
	if( ptr != string::npos )
		str.replace( ptr, substr.length(), newstr );
	return ptr != string::npos;
}

/** in string str, erase the first instance of substr
 *  returns true if an instance of substr was erased */
inline bool FindAndErase(string& str,
			 const string& substr)
{
	return FindAndReplace( str, substr, "" );
}

/**csphot returns photoionization cross section from opacity stage using std pointers 
\param inu INU is array index pointing to frequency where opacity is to be evaluated on f not c scale
\param ithr ITHR is pointer to threshold
\param iofset IOFSET is offset as defined in opac0
*/
double csphot(long int inu, long int ithr, long int iofset);

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

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

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

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

/**uncaps convert input command line (through eol) to all lowercase 
\param chCard - line image as string of characters */
void uncaps(char *chCard );

/**caps convert input command line (through eol) to ALL CAPS 
\param chCard - line image as string of characters */
void caps(char *chCard );

/**e2 second exponential integral 
 \param t optical depth argument */
double e2(
	  double t );

/**ee1 first exponential integral
	\param x optical depth argument, returns e1(tau) */
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.	
 \param x  optical depth argument, returns e1(tau) * exp(x) */
double ee1_safe(double x);

 /**FFmtRead - the free-format number reader
  \param *chCard string giving the line image
  \param *ipnt the index for the character in the string where we shall start
  \param last the number of characters in the string - do not search beyond it
  \param *lgEOL true if hit end of line with no number 
 */ 
double FFmtRead(const char *chCard, 
		long int *ipnt, 
		long int last, 
		bool *lgEOL);

 /**nMatch 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 
	  \param *chKey
	  \param *chCard
 */ 
long nMatch(const char *chKey, 
	    const char *chCard);

/** GetQuote get any name between double quotes off command line
 * return string as chLabel, is null terminated 
 * returns zero for success, 1 for did not find double quotes 
\param *chLabel null terminated string between quotes
\param *chCard input line, imagae, we set string between quotes to spaces
\param lgABORT if true then abort if no double quotes found, if false then
 return null string in this case,
\return 0 if found double quotes, 1 if did not, string between quotes set to spaces
*/ 
int GetQuote( char *chLabel, char *chCard, char *chCardRaw, bool lgABORT );

// these are safe versions of strstr, strchr, etc to work around a deficiency in glibc
inline const char *strstr_s(const char *haystack, const char *needle)
{
	return const_cast<const char *>(strstr(haystack, needle));
}

inline char *strstr_s(char *haystack, const char *needle)
{
	return const_cast<char *>(strstr(haystack, needle));
}

inline const char *strchr_s(const char *s, int c)
{
	return const_cast<const char *>(strchr(s, c));
}

inline char *strchr_s(char *s, int c)
{
	return const_cast<char *>(strchr(s, c));
}

/** ipow
\return  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 
\param *fmt
\param val
*/
// prevent compiler warnings on non-MS systems
#ifdef _MSC_VER
char *PrintEfmt(const char *fmt, double val );
#else
#define PrintEfmt( F, V ) F, V
#endif

/** this is -ln of smallest number sexp can handle */
const double SEXP_LIMIT = 84.;
/** this is -ln of smallest number dsexp can handle */
const double DSEXP_LIMIT = 680.;

/** sexp safe exponential function */
sys_float sexp(sys_float x);
double sexp(double x);

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

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

// safe version of getline() that correctly handles all types of EOL lf, crlf and cr...
istream& SafeGetline(istream& is, string& t);

// Define integration methods
typedef enum { Gaussian32, Legendre } methods;

// define an integrator class.  Currently hard-wired to 32-point Gaussian
template<typename Integrand, methods Method>
class Integrator
{
public:
	double numPoints, weights[16], c[16];

	Integrator( void )
	{
		numPoints = 16;
		double weights_temp[16] = {
			.35093050047350483e-2, .81371973654528350e-2, .12696032654631030e-1, .17136931456510717e-1,
			.21417949011113340e-1, .25499029631188088e-1, .29342046739267774e-1, .32911111388180923e-1,
			.36172897054424253e-1, .39096947893535153e-1, .41655962113473378e-1, .43826046502201906e-1,
			.45586939347881942e-1, .46922199540402283e-1, .47819360039637430e-1, .48270044257363900e-1};

		double c_temp[16] = {
			.498631930924740780, .49280575577263417, .4823811277937532200, .46745303796886984000,
			.448160577883026060, .42468380686628499, .3972418979839712000, .36609105937014484000,
			.331522133465107600, .29385787862038116, .2534499544661147000, .21067563806531767000,
			.165934301141063820, .11964368112606854, .7223598079139825e-1, .24153832843869158e-1};

		for( long i=0; i<numPoints; i++ )
		{
			weights[i] = weights_temp[i];
			c[i] = c_temp[i];
		}
		return;
	};
	double sum(double min, double max, Integrand func)
	{
		ASSERT( Method == Gaussian32 );
		double a = 0.5*(max+min),
		b = max-min,
		total = 0.;

		for( long i=0; i< numPoints; i++ )
			total += b * weights[i] * ( func(a+b*c[i]) + func(a-b*c[i]) );

		return total;
	}
};

 /**			 
 32 point gaussian quadrature integration
 \param xl lower limit to integration
 \param xu - upper limit to integration
 \param (*fct) - pointer to routine to be integrated, arg is x val<BR>
 */ 
double qg32( double, double, double(*)(double) );
/* 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
   \param x[] input array to be sorted	   
   \param  n number of values in x 	   
   \param  iperm[]  permutation output array
   \param  kflag flag saying what to do - 1 sorts into increasing order, not changing
   \param  kflag the original routine 
   \param  *ier error condition, should be 0
 */ 
void spsort( realnum x[], long int n, long int iperm[], int kflag, int *ier);

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

/* */
#ifdef _MSC_VER
	/* disable warning that conditional expression is constant, true or false in if */
#	pragma warning( disable : 4127 )
	/* disable strcat warning */
#	pragma warning( disable : 4996 )
	/* disable bogus underflow warning in MS VS*/
#	pragma warning( disable : 4056 )
	/* disable "inline function removed since not used", MS VS*/
#	pragma warning( disable : 4514 )
	/* disable "assignment operator could not be generated", cddefines.h
	 * line 126 */
#	pragma warning( disable : 4512 )
#endif
#ifdef __INTEL_COMPILER
#	pragma warning( disable : 1572 )/**< disable warning that floating-point comparisons are unreliable */
#endif
/* */

/*lint +e129 these resolve several issues pclint has with my system headers */
/*lint +e78 */
/*lint +e830 */
/*lint +e38 */
/*lint +e148 */
/*lint +e114 */
/*lint +e18 */
/*lint +e49 */

#endif /* CDDEFINES_H_ */