cddefines.h

Go to the documentation of this file.

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

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

extern FILE *ioQQQ;

extern FILE *ioStdin;

extern FILE *ioMAP;

extern FILE* ioPrnErr;

extern bool lgAbort;

extern bool lgTestCodeCalled; 

extern bool lgTestCodeEnabled;

extern bool lgPrnErr;

extern long int nzone;

extern double fnzone;

extern long int iteration;

extern const double ZeroNum;

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

const int FILENAME_PATH_LENGTH = 200; 

const int FILENAME_PATH_LENGTH_2  = FILENAME_PATH_LENGTH*2; 

const int INPUT_LINE_LENGTH = 2000;

const int LIMELM = 30;

const int NISO = 2;

const int NHYDRO_MAX_LEVEL = 401;

const double MAX_DENSITY = 1.e24;

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

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

double fudge(long int ipnt);

void broken(void);

void fixit(void);

void CodeReview(void);

void TestCode(void);

void *MyMalloc(size_t size, const char *file, int line);

void *MyCalloc(size_t num, size_t size);

void *MyRealloc(void *p, size_t size);

void MyAssert(const char *file, int line, const char *comment);

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

void ShowMe(void);

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();
}

NORETURN void BadRead(void);

int dbg_printf(int debug, const char *fmt, ...);

int dprintf(FILE *fp, const char *format, ...);

char *read_whole_line( char *chLine , int nChar , FILE *ioIN );

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

#ifndef NDEBUG
#       define DEBUG
#else
#       undef DEBUG
#endif

#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

#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

#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

#       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

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;
        }               
};

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'; }
/* */

inline bool is_odd( int j ) { return (j&1) == 1; }
inline bool is_odd( long j ) { return (j&1L) == 1L; }
/* */

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

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

#define MIN2 min
/* */

#undef MIN3

#define MIN3(a,b,c) (min(min(a,b),c))
/* */

#undef MIN4

#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

#define MAX2 max
/* */

#undef MAX3

#define MAX3(a,b,c) (max(max(a,b),c))
/* */

#undef MAX4

#define MAX4(a,b,c,d) (max(max(a,b),max(c,d)))
/* */

template<class T>
inline T sign( T x, T y )
{
        return ( y < T() ) ? -abs(x) : abs(x);
}
/* */

template<class T>
inline int sign3( T x ) { return ( x < T() ) ? -1 : ( ( x > T() ) ? 1 : 0 ); }
/* */

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 );
}

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

#define POW2 pow2
template<class T>
inline T pow2(T a) { return a*a; }
/* */

#undef POW3

#define POW3 pow3
template<class T>
inline T pow3(T a) { return a*a*a; }
/* */

#undef POW4

#define POW4 pow4
template<class T>
inline T pow4(T a) { T b = a*a; return b*b; }
/* */

#undef SDIV

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

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() );
}

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)) );
}

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]);
}

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;
        double CoolTotal;
        bool lgActive;
        /* maximum wavenumber in chianti */
        double maxWN;
        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,
 *
 **************************************************************************/

typedef enum { SPM_RELAX, SPM_KEEP_EMPTY, SPM_STRICT } split_mode;

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

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;
}

inline bool FindAndErase(string& str,
                         const string& substr)
{
        return FindAndReplace( str, substr, "" );
}

double csphot(long int inu, long int ithr, long int iofset);

double RandGauss(double xMean, double s );

double MyGaussRand( double PctUncertainty );

double AnuUnit(realnum energy);

void cap4(char *chCAP , const char *chLab);

void uncaps(char *chCard );

void caps(char *chCard );

double e2(
          double t );

double ee1(double x);

double ee1_safe(double x);

double FFmtRead(const char *chCard, 
                long int *ipnt, 
                long int last, 
                bool *lgEOL);

long nMatch(const char *chKey, 
            const char *chCard);

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));
}

long int ipow( long, long );

void PrintE82( FILE*, double );

void PrintE71( FILE*, double );

void PrintE93( FILE*, double );

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

const double SEXP_LIMIT = 84.;
const double DSEXP_LIMIT = 680.;

sys_float sexp(sys_float x);
double sexp(double x);

double dsexp(double x);

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;
        }
};

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


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 )
#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_ */