00001
00002
00003
00004 #ifndef CDDEFINES_H_
00005 #define CDDEFINES_H_
00006
00007 #include "cdstd.h"
00008
00009 #ifdef _MSC_VER
00010
00011 # pragma warning( disable : 4127 )
00012
00013 # ifndef WIN32_LEAN_AND_MEAN
00014 # define WIN32_LEAN_AND_MEAN
00015 # endif
00016 #endif
00017
00018 #ifdef __clang__
00019
00020 #pragma clang diagnostic ignored "-Wmismatched-tags"
00021 #endif
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #include <cstdio>
00034 #include <cstdlib>
00035 #include <cctype>
00036
00037
00038 #define _USE_MATH_DEFINES
00039 #include <cmath>
00040 #include <cassert>
00041 #include <cstring>
00042 #include <cfloat>
00043 #include <climits>
00044 #include <ctime>
00045 #if defined(__sun) && defined(__SUNPRO_CC)
00046
00047 #include <signal.h>
00048 #else
00049 #include <csignal>
00050 #endif
00051
00052 #include <limits>
00053 #include <string>
00054 #include <sstream>
00055 #include <iomanip>
00056 #include <vector>
00057 #include <valarray>
00058 #include <complex>
00059 #include <map>
00060 #include <memory>
00061 #include <stdexcept>
00062 #include <algorithm>
00063 #include <fstream>
00064 #include <bitset>
00065 #ifdef DMALLOC
00066 #include <dmalloc.h>
00067 #endif
00068
00069
00070 #if defined(_MSC_VER) && !defined(SYS_CONFIG)
00071 #define SYS_CONFIG "cloudyconfig_vs.h"
00072 #endif
00073
00074
00075 #ifdef SYS_CONFIG
00076 #include SYS_CONFIG
00077 #else
00078 #include "cloudyconfig.h"
00079 #endif
00080
00081 #ifdef MPI_GRID_RUN
00082 #define MPI_ENABLED
00083 #endif
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093 using namespace std;
00094
00095 #undef STATIC
00096 #ifdef USE_GPROF
00097 #define STATIC
00098 #else
00099 #define STATIC static
00100 #endif
00101
00102 #ifdef FLT_IS_DBL
00103 typedef double realnum;
00104 #else
00105 typedef float realnum;
00106 #endif
00107
00108 typedef float sys_float;
00109
00110 #define float PLEASE_USE_REALNUM_NOT_FLOAT
00111
00112
00113 template<bool> struct StaticAssertFailed;
00114 template<> struct StaticAssertFailed<true> {};
00115 #define STATIC_ASSERT(x) ((void)StaticAssertFailed< (x) == true >())
00116
00117 typedef enum {
00118 ES_SUCCESS=0,
00119 ES_FAILURE=1,
00120 ES_WARNINGS,
00121 ES_BOTCHES,
00122 ES_CLOUDY_ABORT,
00123 ES_BAD_ASSERT,
00124 ES_BAD_ALLOC,
00125 ES_OUT_OF_RANGE,
00126 ES_USER_INTERRUPT,
00127 ES_TERMINATION_REQUEST,
00128 ES_ILLEGAL_INSTRUCTION,
00129 ES_FP_EXCEPTION,
00130 ES_SEGFAULT,
00131 ES_BUS_ERROR,
00132 ES_UNKNOWN_SIGNAL,
00133 ES_UNKNOWN_EXCEPTION,
00134 ES_TOP
00135 } exit_type;
00136
00137
00138
00139 #undef EXIT_SUCCESS
00140 #define EXIT_SUCCESS ES_SUCCESS
00141 #undef EXIT_FAILURE
00142 #define EXIT_FAILURE ES_FAILURE
00143
00144
00145
00146
00147
00148 #include "cpu.h"
00149
00150
00151
00169
00170
00171
00172
00173
00174 template<typename T> class Singleton
00175 {
00176 public:
00177 static T& Inst()
00178 {
00179 static T instance;
00180 return instance;
00181 }
00182 };
00183
00184
00185
00186
00187
00188
00189
00190
00194 extern FILE *ioQQQ;
00195
00196 extern FILE *ioStdin;
00197
00198 extern FILE *ioMAP;
00199
00202 extern FILE* ioPrnErr;
00203
00205 extern bool lgAbort;
00206
00209 extern bool lgTestCodeCalled;
00210
00213 extern bool lgTestCodeEnabled;
00214
00217 extern bool lgPrnErr;
00218
00221 extern long int nzone;
00222
00224 extern double fnzone;
00225
00228 extern long int iteration;
00229
00236 extern const double ZeroNum;
00237
00238
00239
00240
00241
00242
00243
00248 const int FILENAME_PATH_LENGTH = 200;
00249
00251 const int FILENAME_PATH_LENGTH_2 = FILENAME_PATH_LENGTH*2;
00252
00256 const int INPUT_LINE_LENGTH = 2000;
00257
00260 const int LIMELM = 30;
00261
00263 const int NISO = 2;
00264
00268 const int NHYDRO_MAX_LEVEL = 401;
00269
00271 const double MAX_DENSITY = 1.e24;
00272
00274 const double DEPTH_OFFSET = 1.e-30;
00275
00276 enum {CHARS_SPECIES=10};
00277 enum {CHARS_ISOTOPE_SYM = 6};
00278
00279
00280
00281 const int ipRecEsc = 2;
00282
00283 const int ipRecNetEsc = 1;
00284
00285 const int ipRecRad = 0;
00290
00291
00292 const int ipPRD = 1;
00293
00294 const int ipCRD = -1;
00295
00296 const int ipCRDW = 2;
00297
00298 const int ipLY_A = -2;
00299
00300 const int ipDEST_K2 = 1;
00301
00302 const int ipDEST_INCOM = 2;
00303
00304 const int ipDEST_SIMPL = 3;
00305
00307 const int ipHYDROGEN = 0;
00308 const int ipHELIUM = 1;
00309 const int ipLITHIUM = 2;
00310 const int ipBERYLLIUM = 3;
00311 const int ipBORON = 4;
00312 const int ipCARBON = 5;
00313 const int ipNITROGEN = 6;
00314 const int ipOXYGEN = 7;
00315 const int ipFLUORINE = 8;
00316 const int ipNEON = 9;
00317 const int ipSODIUM = 10;
00318 const int ipMAGNESIUM = 11;
00319 const int ipALUMINIUM = 12;
00320 const int ipSILICON = 13;
00321 const int ipPHOSPHORUS = 14;
00322 const int ipSULPHUR = 15;
00323 const int ipCHLORINE = 16;
00324 const int ipARGON = 17;
00325 const int ipPOTASSIUM = 18;
00326 const int ipCALCIUM = 19;
00327 const int ipSCANDIUM = 20;
00328 const int ipTITANIUM = 21;
00329 const int ipVANADIUM = 22;
00330 const int ipCHROMIUM = 23;
00331 const int ipMANGANESE = 24;
00332 const int ipIRON = 25;
00333 const int ipCOBALT = 26;
00334 const int ipNICKEL = 27;
00335 const int ipCOPPER = 28;
00336 const int ipZINC = 29;
00337 const int ipKRYPTON = 35;
00338
00339
00340
00341
00342
00343
00344
00353 double fudge(long int ipnt);
00354
00358 void broken(void);
00359
00362 void fixit(void);
00363
00365 void CodeReview(void);
00366
00368 void TestCode(void);
00369
00375 void *MyMalloc(size_t size, const char *file, int line);
00376
00381 void *MyCalloc(size_t num, size_t size);
00382
00387 void *MyRealloc(void *p, size_t size);
00388
00393 void MyAssert(const char *file, int line, const char *comment);
00394
00396 void cdPrepareExit(exit_type);
00397
00398 class cloudy_exit
00399 {
00400 const char* p_routine;
00401 const char* p_file;
00402 long p_line;
00403 exit_type p_exit;
00404 public:
00405 cloudy_exit(const char* routine, const char* file, long line, exit_type exit_code)
00406 {
00407 p_routine = routine;
00408 p_file = file;
00409 p_line = line;
00410 p_exit = exit_code;
00411 }
00412 virtual ~cloudy_exit() throw()
00413 {
00414 p_routine = NULL;
00415 p_file = NULL;
00416 }
00417 const char* routine() const throw()
00418 {
00419 return p_routine;
00420 }
00421 const char* file() const throw()
00422 {
00423 return p_file;
00424 }
00425 long line() const
00426 {
00427 return p_line;
00428 }
00429 exit_type exit_status() const
00430 {
00431 return p_exit;
00432 }
00433 };
00434
00435
00436 #define cdEXIT( FAIL ) throw cloudy_exit( __func__, __FILE__, __LINE__, FAIL )
00437
00438
00439 #define puts( STR ) Using_puts_before_cdEXIT_is_no_longer_needed
00440
00442 void ShowMe(void);
00443
00445 NORETURN void TotalInsanity(void);
00446
00447
00448
00449
00450 template<class T>
00451 T TotalInsanityAsStub()
00452 {
00453
00454 if( ZeroNum == 0. )
00455 TotalInsanity();
00456 else
00457 return T();
00458 }
00459
00461 NORETURN void BadRead(void);
00462
00466 int dbg_printf(int debug, const char *fmt, ...);
00467
00469 int dprintf(FILE *fp, const char *format, ...);
00470
00479 char *read_whole_line( char *chLine , int nChar , FILE *ioIN );
00480
00481
00482
00483
00484
00485
00486
00490 #ifndef NDEBUG
00491 # define DEBUG
00492 #else
00493 # undef DEBUG
00494 #endif
00495
00497 #if defined(malloc)
00498
00499
00500 # define MALLOC(exp) (malloc(exp))
00501 #else
00502
00503 # define MALLOC(exp) (MyMalloc(exp,__FILE__, __LINE__))
00504 #endif
00505
00507 #if defined(calloc)
00508
00509 # define CALLOC calloc
00510 #else
00511
00512 # define CALLOC MyCalloc
00513 #endif
00514
00516 #if defined(realloc)
00517
00518 # define REALLOC realloc
00519 #else
00520
00521 # define REALLOC MyRealloc
00522 #endif
00523
00524 class bad_signal
00525 {
00526 int p_sig;
00527 public:
00528 explicit bad_signal(int sig)
00529 {
00530 p_sig = sig;
00531 }
00532 virtual ~bad_signal() throw() {}
00533 int sig() const throw()
00534 {
00535 return p_sig;
00536 }
00537 };
00538
00539 class bad_assert
00540 {
00541 const char* p_file;
00542 long p_line;
00543 const char* p_comment;
00544 public:
00545 bad_assert(const char* file, long line, const char* comment);
00546 void print(void) const
00547 {
00548 fprintf(ioQQQ,"DISASTER Assertion failure at %s:%ld\n%s\n",
00549 p_file, p_line, p_comment);
00550 }
00551 virtual ~bad_assert() throw()
00552 {
00553 p_file = NULL;
00554 }
00555 const char* file() const throw()
00556 {
00557 return p_file;
00558 }
00559 long line() const throw()
00560 {
00561 return p_line;
00562 }
00563 const char *comment() const throw()
00564 {
00565 return p_comment;
00566 }
00567 };
00568
00569
00570
00571
00572
00573
00574
00575 #undef ASSERT
00576 #ifndef OLD_ASSERT
00577 # if NDEBUG
00578 # define ASSERT(exp) ((void)0)
00579 # else
00580 # define ASSERT(exp) \
00581 do { \
00582 if (UNLIKELY(!(exp))) \
00583 { \
00584 bad_assert aa(__FILE__,__LINE__,"Failed: " #exp); \
00585 if( cpu.i().lgAssertAbort() ) \
00586 { \
00587 aa.print(); \
00588 abort(); \
00589 } \
00590 else \
00591 throw aa; \
00592 } \
00593 } while( 0 )
00594 # endif
00595 #else
00596
00597 # ifdef NDEBUG
00598 # define ASSERT(exp) ((void)0)
00599 # else
00600 # define ASSERT(exp) \
00601 do { \
00602 if (!(exp)) \
00603 MyAssert(__FILE__, __LINE__, "Failed: " #exp); \
00604 } while( 0 )
00605 # endif
00606 #endif
00607
00608 #define MESSAGE_ASSERT(msg, exp) ASSERT( (msg) ? (exp) : false )
00609
00610 inline NORETURN void OUT_OF_RANGE(const char* str)
00611 {
00612 if( cpu.i().lgAssertAbort() )
00613 abort();
00614 else
00615 throw out_of_range( str );
00616 }
00617
00618
00619
00620
00621 #undef isnan
00622 #define isnan MyIsnan
00623
00626 class t_debug : public Singleton<t_debug>
00627 {
00628 friend class Singleton<t_debug>;
00629 FILE *p_fp;
00630 int p_callLevel;
00631 protected:
00632 t_debug() : p_fp(stderr)
00633 {
00634 p_callLevel = 0;
00635 }
00636 public:
00637 void enter(const char *name)
00638 {
00639 ++p_callLevel;
00640 fprintf(p_fp,"%*c%s\n",p_callLevel,'>',name);
00641 }
00642 void leave(const char *name)
00643 {
00644 fprintf(p_fp,"%*c%s\n",p_callLevel,'<',name);
00645 --p_callLevel;
00646 }
00647 };
00648
00651 class t_nodebug : public Singleton<t_nodebug>
00652 {
00653 friend class Singleton<t_nodebug>;
00654 protected:
00655 t_nodebug() {}
00656 public:
00657 void enter(const char *) const {}
00658 void leave(const char *) const {}
00659 };
00660
00661 template<class Trace>
00662 class debugtrace
00663 {
00664 const char *p_name;
00665 public:
00666 explicit debugtrace(const char *funcname)
00667 {
00668 p_name = funcname;
00669 Trace::Inst().enter(p_name);
00670 }
00671 ~debugtrace()
00672 {
00673 Trace::Inst().leave(p_name);
00674 p_name = NULL;
00675 }
00676 const char* name() const
00677 {
00678 return p_name;
00679 }
00680 };
00681
00682 #ifdef DEBUG_FUN
00683 #define DEBUG_ENTRY( funcname ) debugtrace<t_debug> DEBUG_ENTRY( funcname )
00684 #else
00685 #ifdef HAVE_FUNC
00686 #define DEBUG_ENTRY( funcname ) ((void)0)
00687 #else
00688 #define DEBUG_ENTRY( funcname ) debugtrace<t_nodebug> DEBUG_ENTRY( funcname )
00689 #endif
00690 #endif
00691
00692
00693 inline char tolower(char c)
00694 {
00695 return static_cast<char>( tolower( static_cast<int>(c) ) );
00696 }
00697 inline unsigned char tolower(unsigned char c)
00698 {
00699 return static_cast<unsigned char>( tolower( static_cast<int>(c) ) );
00700 }
00701
00702 inline char toupper(char c)
00703 {
00704 return static_cast<char>( toupper( static_cast<int>(c) ) );
00705 }
00706 inline unsigned char toupper(unsigned char c)
00707 {
00708 return static_cast<unsigned char>( toupper( static_cast<int>(c) ) );
00709 }
00710
00711
00712 inline char TorF( bool l ) { return l ? 'T' : 'F'; }
00713
00714
00716 inline bool is_odd( int j ) { return (j&1) == 1; }
00717 inline bool is_odd( long j ) { return (j&1L) == 1L; }
00718
00719
00721 inline long nint( double x ) { return static_cast<long>( (x < 0.) ? x-0.5 : x+0.5 ); }
00722
00723
00724
00725 inline long min( int a, long b ) { long c = a; return ( (c < b) ? c : b ); }
00726 inline long min( long a, int b ) { long c = b; return ( (a < c) ? a : c ); }
00727 inline double min( sys_float a, double b ) { double c = a; return ( (c < b) ? c : b ); }
00728 inline double min( double a, sys_float b ) { double c = b; return ( (a < c) ? a : c ); }
00729
00730
00731 #ifndef HAVE_POWI
00732
00733 double powi( double , long int );
00734 #endif
00735
00736
00737 #ifndef HAVE_POW_DOUBLE_INT
00738 inline double pow( double x, int i ) { return powi( x, long(i) ); }
00739 #endif
00740
00741 #ifndef HAVE_POW_DOUBLE_LONG
00742 inline double pow( double x, long i ) { return powi( x, i ); }
00743 #endif
00744
00745 #ifndef HAVE_POW_FLOAT_INT
00746 inline sys_float pow( sys_float x, int i ) { return sys_float( powi( double(x), long(i) ) ); }
00747 #endif
00748
00749 #ifndef HAVE_POW_FLOAT_LONG
00750 inline sys_float pow( sys_float x, long i ) { return sys_float( powi( double(x), i ) ); }
00751 #endif
00752
00753 #ifndef HAVE_POW_FLOAT_DOUBLE
00754 inline double pow( sys_float x, double y ) { return pow( double(x), y ); }
00755 #endif
00756
00757 #ifndef HAVE_POW_DOUBLE_FLOAT
00758 inline double pow( double x, sys_float y ) { return pow( x, double(y) ); }
00759 #endif
00760
00761 #undef MIN2
00762
00763 #define MIN2 min
00764
00765
00766 #undef MIN3
00767
00768 #define MIN3(a,b,c) (min(min(a,b),c))
00769
00770
00771 #undef MIN4
00772
00773 #define MIN4(a,b,c,d) (min(min(a,b),min(c,d)))
00774
00775
00776
00777 inline long max( int a, long b ) { long c = a; return ( (c > b) ? c : b ); }
00778 inline long max( long a, int b ) { long c = b; return ( (a > c) ? a : c ); }
00779 inline double max( sys_float a, double b ) { double c = a; return ( (c > b) ? c : b ); }
00780 inline double max( double a, sys_float b ) { double c = b; return ( (a > c) ? a : c ); }
00781
00782 #undef MAX2
00783
00784 #define MAX2 max
00785
00786
00787 #undef MAX3
00788
00789 #define MAX3(a,b,c) (max(max(a,b),c))
00790
00791
00792 #undef MAX4
00793
00794 #define MAX4(a,b,c,d) (max(max(a,b),max(c,d)))
00795
00796
00801 template<class T>
00802 inline T sign( T x, T y )
00803 {
00804 return ( y < T() ) ? -abs(x) : abs(x);
00805 }
00806
00807
00809 template<class T>
00810 inline int sign3( T x ) { return ( x < T() ) ? -1 : ( ( x > T() ) ? 1 : 0 ); }
00811
00812
00814 inline bool fp_equal( sys_float x, sys_float y, int n=3 )
00815 {
00816 #ifdef _MSC_VER
00817
00818 # pragma warning( disable : 4127 )
00819 #endif
00820 ASSERT( n >= 1 );
00821
00822 if( isnan(x) || isnan(y) )
00823 return false;
00824 int sx = sign3(x);
00825 int sy = sign3(y);
00826
00827 if( sx == 0 && sy == 0 )
00828 return true;
00829
00830 if( sx*sy != 1 )
00831 return false;
00832 x = abs(x);
00833 y = abs(y);
00834 return ( 1.f - min(x,y)/max(x,y) < ((sys_float)n+0.1f)*FLT_EPSILON );
00835 }
00836
00837 inline bool fp_equal( double x, double y, int n=3 )
00838 {
00839 ASSERT( n >= 1 );
00840
00841 if( isnan(x) || isnan(y) )
00842 return false;
00843 int sx = sign3(x);
00844 int sy = sign3(y);
00845
00846 if( sx == 0 && sy == 0 )
00847 return true;
00848
00849 if( sx*sy != 1 )
00850 return false;
00851 x = abs(x);
00852 y = abs(y);
00853 return ( 1. - min(x,y)/max(x,y) < ((double)n+0.1)*DBL_EPSILON );
00854 }
00855
00856 inline bool fp_equal_tol( sys_float x, sys_float y, sys_float tol )
00857 {
00858 ASSERT( tol > 0.f );
00859
00860 if( isnan(x) || isnan(y) )
00861 return false;
00862
00863 ASSERT( tol >= FLT_EPSILON*max(abs(x),abs(y)) );
00864 return ( abs( x-y ) <= tol );
00865 }
00866
00867 inline bool fp_equal_tol( double x, double y, double tol )
00868 {
00869 ASSERT( tol > 0. );
00870
00871 if( isnan(x) || isnan(y) )
00872 return false;
00873
00874 ASSERT( tol >= DBL_EPSILON*max(abs(x),abs(y)) );
00875 return ( abs( x-y ) <= tol );
00876 }
00877
00879 inline bool fp_bound( sys_float lo, sys_float x, sys_float hi, int n=3 )
00880 {
00881 ASSERT( n >= 1 );
00882
00883 if( isnan(x) || isnan(lo) || isnan(hi) )
00884 return false;
00885 if( fp_equal(lo,hi,n) )
00886 return fp_equal(0.5f*(lo+hi),x,n);
00887 if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -((sys_float)n+0.1f)*FLT_EPSILON )
00888 return false;
00889 return true;
00890 }
00891 inline bool fp_bound( double lo, double x, double hi, int n=3 )
00892 {
00893 ASSERT( n >= 1 );
00894
00895 if( isnan(x) || isnan(lo) || isnan(hi) )
00896 return false;
00897 if( fp_equal(lo,hi,n) )
00898 return fp_equal(0.5*(lo+hi),x,n);
00899 if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -((double)n+0.1)*DBL_EPSILON )
00900 return false;
00901 return true;
00902 }
00903 inline bool fp_bound_tol( sys_float lo, sys_float x, sys_float hi, sys_float tol )
00904 {
00905 ASSERT( tol > 0.f );
00906
00907 if( isnan(x) || isnan(lo) || isnan(hi) )
00908 return false;
00909 if( fp_equal_tol(lo,hi,tol) )
00910 return fp_equal_tol(0.5f*(lo+hi),x,tol);
00911 if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -tol )
00912 return false;
00913 return true;
00914 }
00915 inline bool fp_bound_tol( double lo, double x, double hi, double tol )
00916 {
00917 ASSERT( tol > 0. );
00918
00919 if( isnan(x) || isnan(lo) || isnan(hi) )
00920 return false;
00921 if( fp_equal_tol(lo,hi,tol) )
00922 return fp_equal_tol(0.5*(lo+hi),x,tol);
00923 if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -tol )
00924 return false;
00925 return true;
00926 }
00927
00928
00929 #undef POW2
00930
00931 #define POW2 pow2
00932 template<class T>
00933 inline T pow2(T a) { return a*a; }
00934
00935
00936 #undef POW3
00937
00938 #define POW3 pow3
00939 template<class T>
00940 inline T pow3(T a) { return a*a*a; }
00941
00942
00943 #undef POW4
00944
00945 #define POW4 pow4
00946 template<class T>
00947 inline T pow4(T a) { T b = a*a; return b*b; }
00948
00949
00950 #undef SDIV
00951
00954 inline sys_float SDIV( sys_float x ) { return ( fabs((double)x) < (double)SMALLFLOAT ) ? (sys_float)SMALLFLOAT : x; }
00955
00956 inline double SDIV( double x ) { return ( fabs(x) < (double)SMALLFLOAT ) ? (double)SMALLFLOAT : x; }
00957
00958
00959
00963 inline sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
00964 {
00965
00966 if( isnan(x) || isnan(y) )
00967 return x/y;
00968 int sx = sign3(x);
00969 int sy = sign3(y);
00970
00971 if( sx == 0 && sy == 0 )
00972 {
00973 if( isnan(res_0by0) )
00974 return x/y;
00975 else
00976 return res_0by0;
00977 }
00978 if( sx == 0 )
00979 return 0.;
00980 if( sy == 0 )
00981 return ( sx < 0 ) ? -FLT_MAX : FLT_MAX;
00982
00983 sys_float ay = abs(y);
00984 if( ay >= 1.f )
00985 return x/y;
00986 else
00987 {
00988
00989 if( abs(x) < ay*FLT_MAX )
00990 return x/y;
00991 else
00992 return ( sx*sy < 0 ) ? -FLT_MAX : FLT_MAX;
00993 }
00994 }
00995
00996 inline sys_float safe_div(sys_float x, sys_float y)
00997 {
00998 return safe_div( x, y, numeric_limits<sys_float>::quiet_NaN() );
00999 }
01000
01004 inline double safe_div(double x, double y, double res_0by0)
01005 {
01006
01007 if( isnan(x) || isnan(y) )
01008 return x/y;
01009 int sx = sign3(x);
01010 int sy = sign3(y);
01011
01012 if( sx == 0 && sy == 0 )
01013 {
01014 if( isnan(res_0by0) )
01015 return x/y;
01016 else
01017 return res_0by0;
01018 }
01019 if( sx == 0 )
01020 return 0.;
01021 if( sy == 0 )
01022 return ( sx < 0 ) ? -DBL_MAX : DBL_MAX;
01023
01024 double ay = abs(y);
01025 if( ay >= 1. )
01026 return x/y;
01027 else
01028 {
01029
01030 if( abs(x) < ay*DBL_MAX )
01031 return x/y;
01032 else
01033 return ( sx*sy < 0 ) ? -DBL_MAX : DBL_MAX;
01034 }
01035 }
01036
01037 inline double safe_div(double x, double y)
01038 {
01039 return safe_div( x, y, numeric_limits<double>::quiet_NaN() );
01040 }
01041
01042 #undef HMRATE
01043
01044
01045
01046
01047
01048 #define HMRATE(a,b,c) hmrate4(a,b,c,phycon.te)
01049
01050 inline double hmrate4( double a, double b, double c, double te )
01051 {
01052 if( b == 0. && c == 0. )
01053 return a;
01054 else if( c == 0. )
01055 return a*pow(te/300.,b);
01056 else if( b == 0. )
01057 return ( c/te <= 50. ) ? a*exp(-c/te) : 0.;
01058 else
01059 return ( c/te <= 50. ) ? a*pow(te/300.,b)*exp(-c/te) : 0.;
01060 }
01061
01062 template<class T>
01063 inline void invalidate_array(T* p, size_t size)
01064 {
01065 if( size > 0 )
01066 memset( p, -1, size );
01067 }
01068
01069 inline void invalidate_array(double* p, size_t size)
01070 {
01071 set_NaN( p, (long)(size/sizeof(double)) );
01072 }
01073
01074 inline void invalidate_array(sys_float* p, size_t size)
01075 {
01076 set_NaN( p, (long)(size/sizeof(sys_float)) );
01077 }
01078
01081 template<class T> inline T* get_ptr(T *v)
01082 {
01083 return v;
01084 }
01085 template<class T> inline T* get_ptr(valarray<T> &v)
01086 {
01087 return &v[0];
01088 }
01089 template<class T> inline T* get_ptr(vector<T> &v)
01090 {
01091 return &v[0];
01092 }
01093 template<class T> inline const T* get_ptr(const valarray<T> &v)
01094 {
01095 return const_cast<const T*>(&const_cast<valarray<T>&>(v)[0]);
01096 }
01097 template<class T> inline const T* get_ptr(const vector<T> &v)
01098 {
01099 return const_cast<const T*>(&const_cast<vector<T>&>(v)[0]);
01100 }
01101
01127 template<class T>
01128 class auto_vec
01129 {
01130 T* ptr;
01131
01132 template<class U>
01133 struct auto_vec_ref
01134 {
01135 U* ptr;
01136
01137 explicit auto_vec_ref( U* p )
01138 {
01139 ptr = p;
01140 }
01141 };
01142
01143 public:
01144 typedef T element_type;
01145
01146
01147
01148 explicit auto_vec( element_type* p = NULL ) throw()
01149 {
01150 ptr = p;
01151 }
01152 auto_vec( auto_vec& p ) throw()
01153 {
01154 ptr = p.release();
01155 }
01156 auto_vec& operator= ( auto_vec& p ) throw()
01157 {
01158 reset( p.release() );
01159 return *this;
01160 }
01161 ~auto_vec() throw()
01162 {
01163 delete[] ptr;
01164 }
01165
01166
01167
01168 element_type& operator[] ( ptrdiff_t n ) const throw()
01169 {
01170 return *(ptr+n);
01171 }
01172 element_type* get() const throw()
01173 {
01174 return ptr;
01175 }
01176
01177 element_type* data() const throw()
01178 {
01179 return ptr;
01180 }
01181 element_type* release() throw()
01182 {
01183 element_type* p = ptr;
01184 ptr = NULL;
01185 return p;
01186 }
01187 void reset( element_type* p = NULL ) throw()
01188 {
01189 if( p != ptr )
01190 {
01191 delete[] ptr;
01192 ptr = p;
01193 }
01194 }
01195
01196
01197
01198 auto_vec( auto_vec_ref<element_type> r ) throw()
01199 {
01200 ptr = r.ptr;
01201 }
01202 auto_vec& operator= ( auto_vec_ref<element_type> r ) throw()
01203 {
01204 if( r.ptr != ptr )
01205 {
01206 delete[] ptr;
01207 ptr = r.ptr;
01208 }
01209 return *this;
01210 }
01211 operator auto_vec_ref<element_type>() throw()
01212 {
01213 return auto_vec_ref<element_type>( this->release() );
01214 }
01215 };
01216
01217 #include "container_classes.h"
01218 #include "iter_track.h"
01219
01220
01221
01222
01223
01224
01225
01226 typedef struct t_species species;
01227
01228 #include "lines_service.h"
01229
01230
01231
01232
01233
01234 struct t_species
01235 {
01236
01237 char *chLabel;
01238
01239 long index;
01240
01241 long numLevels_max;
01242
01243 long numLevels_local;
01244
01245 realnum fmolweight;
01246
01247 bool lgMolecular;
01248
01249 bool lgLAMDA;
01250
01251 realnum fracType;
01252
01253 realnum fracIsotopologue;
01255 double CoolTotal;
01257 bool lgActive;
01258
01259 double maxWN;
01261 bool lgLTE;
01262 };
01263
01264
01265
01266 typedef struct t_CollRatesArray
01267 {
01268
01269 vector<double> temps;
01270
01271 multi_arr<double,3> collrates;
01272
01273 } CollRateCoeffArray ;
01274
01275
01276
01277 typedef struct t_CollSplinesArray
01278 {
01279
01280
01281
01282
01283
01284 double *collspline;
01285 double *SplineSecDer;
01286
01287 long nSplinePts;
01288 long intTranType;
01289 double EnergyDiff;
01290 double ScalingParam;
01291
01292 } CollSplinesArray ;
01293
01294
01295
01296 typedef struct t_StoutColls
01297 {
01298
01299 long ntemps;
01300
01301 double *temps;
01302
01303 double *collstrs;
01304
01305 bool lgIsRate;
01306
01307 } StoutColls ;
01308
01309 #include "physconst.h"
01310
01311
01312
01313
01314
01315
01316
01323 typedef enum { SPM_RELAX, SPM_KEEP_EMPTY, SPM_STRICT } split_mode;
01324
01326 void Split(const string& str,
01327 const string& sep,
01328 vector<string>& lst,
01329 split_mode mode);
01330
01333 inline bool FindAndReplace(string& str,
01334 const string& substr,
01335 const string& newstr)
01336 {
01337 string::size_type ptr = str.find( substr );
01338 if( ptr != string::npos )
01339 str.replace( ptr, substr.length(), newstr );
01340 return ptr != string::npos;
01341 }
01342
01345 inline bool FindAndErase(string& str,
01346 const string& substr)
01347 {
01348 return FindAndReplace( str, substr, "" );
01349 }
01350
01356 double csphot(long int inu, long int ithr, long int iofset);
01357
01362 double RandGauss(double xMean, double s );
01363
01367 double MyGaussRand( double PctUncertainty );
01368
01370 double AnuUnit(realnum energy);
01371
01376 void cap4(char *chCAP , const char *chLab);
01377
01380 void uncaps(char *chCard );
01381
01384 void caps(char *chCard );
01385
01388 double e2(
01389 double t );
01390
01393 double ee1(double x);
01394
01398 double ee1_safe(double x);
01399
01406 double FFmtRead(const char *chCard,
01407 long int *ipnt,
01408 long int last,
01409 bool *lgEOL);
01410
01416 long nMatch(const char *chKey,
01417 const char *chCard);
01418
01428 int GetQuote( char *chLabel, char *chCard, char *chCardRaw, bool lgABORT );
01429
01430
01431 inline const char *strstr_s(const char *haystack, const char *needle)
01432 {
01433 return const_cast<const char *>(strstr(haystack, needle));
01434 }
01435
01436 inline char *strstr_s(char *haystack, const char *needle)
01437 {
01438 return const_cast<char *>(strstr(haystack, needle));
01439 }
01440
01441 inline const char *strchr_s(const char *s, int c)
01442 {
01443 return const_cast<const char *>(strchr(s, c));
01444 }
01445
01446 inline char *strchr_s(char *s, int c)
01447 {
01448 return const_cast<char *>(strchr(s, c));
01449 }
01450
01453 long int ipow( long, long );
01454
01457 void PrintE82( FILE*, double );
01458
01460 void PrintE71( FILE*, double );
01461
01463 void PrintE93( FILE*, double );
01464
01470
01471 #ifdef _MSC_VER
01472 char *PrintEfmt(const char *fmt, double val );
01473 #else
01474 #define PrintEfmt( F, V ) F, V
01475 #endif
01476
01478 const double SEXP_LIMIT = 84.;
01480 const double DSEXP_LIMIT = 680.;
01481
01483 sys_float sexp(sys_float x);
01484 double sexp(double x);
01485
01490 double dsexp(double x);
01491
01496 double plankf(long int ip);
01497
01498
01499 typedef enum { Gaussian32, Legendre } methods;
01500
01501
01502 template<typename Integrand, methods Method>
01503 class Integrator
01504 {
01505 public:
01506 double numPoints, weights[16], c[16];
01507
01508 Integrator( void )
01509 {
01510 numPoints = 16;
01511 double weights_temp[16] = {
01512 .35093050047350483e-2, .81371973654528350e-2, .12696032654631030e-1, .17136931456510717e-1,
01513 .21417949011113340e-1, .25499029631188088e-1, .29342046739267774e-1, .32911111388180923e-1,
01514 .36172897054424253e-1, .39096947893535153e-1, .41655962113473378e-1, .43826046502201906e-1,
01515 .45586939347881942e-1, .46922199540402283e-1, .47819360039637430e-1, .48270044257363900e-1};
01516
01517 double c_temp[16] = {
01518 .498631930924740780, .49280575577263417, .4823811277937532200, .46745303796886984000,
01519 .448160577883026060, .42468380686628499, .3972418979839712000, .36609105937014484000,
01520 .331522133465107600, .29385787862038116, .2534499544661147000, .21067563806531767000,
01521 .165934301141063820, .11964368112606854, .7223598079139825e-1, .24153832843869158e-1};
01522
01523 for( long i=0; i<numPoints; i++ )
01524 {
01525 weights[i] = weights_temp[i];
01526 c[i] = c_temp[i];
01527 }
01528 return;
01529 };
01530 double sum(double min, double max, Integrand func)
01531 {
01532 ASSERT( Method == Gaussian32 );
01533 double a = 0.5*(max+min),
01534 b = max-min,
01535 total = 0.;
01536
01537 for( long i=0; i< numPoints; i++ )
01538 total += b * weights[i] * ( func(a+b*c[i]) + func(a-b*c[i]) );
01539
01540 return total;
01541 }
01542 };
01543
01550 double qg32( double, double, double(*)(double) );
01551
01552
01553
01554
01564 void spsort( realnum x[], long int n, long int iperm[], int kflag, int *ier);
01565
01566
01567
01568
01569
01570
01571
01572
01573 #ifdef _MSC_VER
01574
01575 # pragma warning( disable : 4127 )
01576
01577 # pragma warning( disable : 4996 )
01578
01579 # pragma warning( disable : 4056 )
01580
01581 # pragma warning( disable : 4514 )
01582
01583
01584 # pragma warning( disable : 4512 )
01585 #endif
01586 #ifdef __INTEL_COMPILER
01587 # pragma warning( disable : 1572 )
01588 #endif
01589
01590
01591
01592
01593
01594
01595
01596
01597
01598
01599
01600 #endif
01601