cddefines.h

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2022 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

// prevent problems on platforms with broken support (such as Mac homebrew g++)
#ifndef HAVE_AVX_INTRIN
#undef __AVX__
#endif

#ifndef HAVE_FMA_INTRIN
#undef __FMA__
#endif

#ifndef HAVE_AVX2_INTRIN
#undef __AVX2__
#endif

#ifndef HAVE_AVX512F_INTRIN
#undef __AVX512F__
#endif

#ifdef __AVX__
#include <immintrin.h>
#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

inline FILE *sys_fopen(const char *path, const char *mode)
{
        return fopen( path, mode );
}
#define fopen PLEASE_USE_open_data_NOT_fopen

//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_DOMAIN_ERROR,         // a vectorized math routine threw a domain error
        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

#ifdef BOUNDS_CHECK
#define lgBOUNDSCHECKVAL true
#else
#define lgBOUNDSCHECKVAL false
#endif

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

class Output
{
public:
        FILE *m_fp;
        // Implicit conversion
   Output(FILE* fp) : m_fp(fp) {}
        FILE *fptr() const
        {
                return m_fp;
        }
};

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 NCHLAB = 20;

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;

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_base(const char* func, const char* file, int line,
                                         const char *reason);

class Fixit
{
public:
        Fixit(const char* func, const char* file, int line,
                                         const char *reason)
        {
                fixit_base(func,file,line,reason);
        }
};

#define fixit(a) \
        do { \
                static Fixit fixit_s(__func__,__FILE__,__LINE__, (a));  \
        } while (0)

void CodeReview(void);

void TestCode(void);

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

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

int fprintf (const Output& stream, const char *format, ...);

int dprintf(const Output & stream, 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

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

#ifdef __SUNPRO_CC
#pragma does_not_return(OUT_OF_RANGE)
#endif

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

double powpq(double x, int p, int q);

/* 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(a,b) min(a,b)
/* */

#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(a,b) max(a,b)
/* */

#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);
        // adjust epsilon value for denormalized numbers
        int exp;
        (void)frexp(max(x,y), &exp);
        sys_float epsilon = ldexp(FLT_EPSILON, max(0,FLT_MIN_EXP-exp));
        return ( 1.f - min(x,y)/max(x,y) < ((sys_float)n+0.1f)*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);
        // adjust epsilon value for denormalized numbers
        int exp;
        (void)frexp(max(x,y), &exp);
        double epsilon = ldexp(DBL_EPSILON, max(0,DBL_MIN_EXP-exp));
        return ( 1. - min(x,y)/max(x,y) < ((double)n+0.1)*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() );
}

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, class U> inline T* get_ptr(vector<T,U> &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, class U> inline const T* get_ptr(const vector<T,U> &v)
{       
        return const_cast<const T*>(&const_cast<vector<T,U>&>(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() );
        }
};

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 FFmtRead(const char *chCard, 
                long int *ipnt, 
                long int last, 
                bool *lgEOL);

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

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

// this routine is fully equivalent to snprintf() apart from the fact that it
// concatenates the output to an existing string rather than replacing it.
size_t sncatf( char* buf, size_t bufSize, const char* fmt, ... );

size_t sncatf( ostringstream& buf, const char* fmt, ... );

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

inline double exp10(double x)
{
        if( x < -330. )
                return 0.;
        else if( x > 310. )
                return pow2(BIGDOUBLE); // +Inf
        else
        {
                double y = trunc(x/3.);
                double z = 3.*y;
                x -= z;
                x *= 2.302585092994045684;       // constant is ln(10)
                x -= 7.90550887243868070612e-3*z; // constant is ln(10)-10/3*ln(2)
                return ldexp(exp(x),10*int(y));
        }
}

inline sys_float exp10f(sys_float x)
{
        if( x < -50.f )
                return 0.f;
        else if( x > 40.f )
                return pow2(BIGFLOAT); // +Inf
        else
        {
                sys_float y = truncf(x/3.f);
                sys_float z = 3.f*y;
                x -= z;
                x *= 2.302585093f;       // constant is ln(10)
                x -= 7.9055088724e-3f*z; // constant is ln(10)-10/3*ln(2)
                return ldexpf(expf(x),10*int(y));
        }
}

inline sys_float exp10(sys_float x)
{
        return exp10f(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);

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