cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 00004 #ifndef _CDDEFINES_H_ 00005 #define _CDDEFINES_H_ 00006 00007 #ifdef _MSC_VER 00008 /* disable warning that conditional expression is constant, true or false in if */ 00009 # pragma warning( disable : 4127 ) 00010 /* we are not using MS foundation class */ 00011 # ifndef WIN32_LEAN_AND_MEAN 00012 # define WIN32_LEAN_AND_MEAN 00013 # endif 00014 #endif 00015 /* these headers are needed by all files */ 00016 /*lint -e129 these resolve several issues pclint has with my system headers */ 00017 /*lint -e78 */ 00018 /*lint -e830 */ 00019 /*lint -e38 */ 00020 /*lint -e148 */ 00021 /*lint -e114 */ 00022 /*lint -e18 */ 00023 /*lint -e49 */ 00024 #include <stdio.h> 00025 #include <stdlib.h> 00026 #include <math.h> 00027 #include <assert.h> 00028 #include <string.h> 00029 #include <float.h> 00030 #include <limits.h> 00031 #include <time.h> 00032 #include <signal.h> 00033 #include <string> 00034 #include <sstream> 00035 #include <vector> 00036 #include <valarray> 00037 #include <complex> 00038 #include <map> 00039 #include <memory> 00040 #include <stdexcept> 00041 #ifdef DMALLOC 00042 #include <dmalloc.h> 00043 #endif 00044 /*lint +e18 */ 00045 /*lint +e49 */ 00046 /*lint +e38 */ 00047 /*lint +e148 */ 00048 /*lint +e830 */ 00049 /*lint +e78 */ 00050 /*lint -e129 */ 00051 00052 using namespace std; 00053 00055 #ifndef EXTERN 00056 #define EXTERN extern 00057 #endif 00058 00059 #undef STATIC 00060 #ifdef USE_GPROF 00061 #define STATIC 00062 #else 00063 #define STATIC static 00064 #endif 00065 00066 #ifdef FLT_IS_DBL 00067 typedef double realnum; 00068 #else 00069 typedef float realnum; 00070 #endif 00071 00072 typedef float sys_float; 00073 // prevent explicit float's from creeping back into the code 00074 #define float PLEASE_USE_REALNUM_NOT_FLOAT 00075 00076 /* make sure this is globally visible as well! */ 00077 #include "cpu.h" 00078 00079 //************************************************************************* 00080 // 00098 // 00099 // This implementation has been obtained from Wikipedia 00100 // 00101 //************************************************************************* 00102 00103 template<typename T> class Singleton 00104 { 00105 public: 00106 static T& Inst() 00107 { 00108 static T instance; // assumes T has a protected default constructor 00109 return instance; 00110 } 00111 }; 00112 00113 /************************************************************************** 00114 * 00115 * these are variables and pointers for output from the code, used everywhere 00116 * declared extern here, and definition is in cddefines.cpp 00117 * 00118 **************************************************************************/ 00119 00123 EXTERN FILE *ioQQQ; 00124 00125 EXTERN FILE *ioStdin; 00126 00127 extern FILE *ioMAP; 00128 00131 EXTERN FILE* ioPrnErr; 00132 00134 EXTERN bool lgAbort; 00135 00138 EXTERN bool lgTestCodeCalled; 00139 00142 EXTERN bool lgTestCodeEnabled; 00143 00146 EXTERN bool lgPrnErr; 00147 00150 EXTERN long int nzone; 00151 00153 EXTERN double fnzone; 00154 00157 EXTERN long int iteration; 00158 00160 EXTERN long int nSimThisCoreload; 00161 00168 extern const double ZeroNum; 00169 00170 /************************************************************************** 00171 * 00172 * these are constants used to dimension several vectors and index arrays 00173 * 00174 **************************************************************************/ 00175 00180 const int FILENAME_PATH_LENGTH = 200; 00181 00183 const int FILENAME_PATH_LENGTH_2 = FILENAME_PATH_LENGTH*2; 00184 00188 const int INPUT_LINE_LENGTH = 200; 00189 00192 const int LIMELM = 30; 00193 00195 const int NISO = 2; 00196 00200 const int NHYDRO_MAX_LEVEL = 401; 00201 00203 const int N_H_MOLEC = 8; 00204 00205 #if 0 00206 00207 const double TEMP_STOP_DEFAULT = 4000.; 00209 const double TEMP_LIMIT_HIGH = 1e10; 00211 const double TEMP_LIMIT_LOW = 2.8; 00212 #endif 00213 00215 const double DEPTH_OFFSET = 1.e-30; 00216 00217 /* these are flags for various colliders that are used across the code */ 00218 enum collider { 00219 ipELECTRON, 00220 ipPROTON, 00221 ipHE_PLUS, 00222 ipALPHA, 00223 //ipATOMH, 00224 //ipHE_2PLUS, 00225 //ipH2_ORTHO, 00226 //ipH2_PARA, 00227 ipNCOLLIDER 00228 }; 00229 00230 /* indices within recombination coefficient array */ 00231 /* ipRecEsc is state specific escape probability*/ 00232 const int ipRecEsc = 2; 00233 /* the net escaping, including destruction by background and optical deepth*/ 00234 const int ipRecNetEsc = 1; 00235 /* ipRecRad is state specific radiative recombination rate*/ 00236 const int ipRecRad = 0; 00241 /* these specify the form of the line redistribution function */ 00242 /* partial redistribution with wings */ 00243 const int ipPRD = 1; 00244 /* complete redistribution, core only, no wings, Hummer's K2 function */ 00245 const int ipCRD = -1; 00246 /* complete redistribution with wings */ 00247 const int ipCRDW = 2; 00248 /* redistribution function for Lya, calls Hummer routine for H-like series only */ 00249 const int ipLY_A = -2; 00250 /* core function for K2 destruction */ 00251 const int ipDEST_K2 = 1; 00252 /* core function for complete redist destruction */ 00253 const int ipDEST_INCOM = 2; 00254 /* core function for simple destruction */ 00255 const int ipDEST_SIMPL = 3; 00256 00259 const int ipH1s = 0; 00260 const int ipH2s = 1; 00261 const int ipH2p = 2; 00262 const int ipH3s = 3; 00263 const int ipH3p = 4; 00264 const int ipH3d = 5; 00265 const int ipH4s = 6; 00266 const int ipH4p = 7; 00267 const int ipH4d = 8; 00268 const int ipH4f = 9; 00269 const int ipH5s = 10; 00270 const int ipH5p = 11; 00271 const int ipH5d = 12; 00272 const int ipH5f = 13; 00273 const int ipH5g = 14; 00274 const int ipH6s = 15; 00275 00278 /* level 1 */ 00279 const int ipHe1s1S = 0; 00280 00281 /* level 2 */ 00282 const int ipHe2s3S = 1; 00283 const int ipHe2s1S = 2; 00284 const int ipHe2p3P0 = 3; 00285 const int ipHe2p3P1 = 4; 00286 const int ipHe2p3P2 = 5; 00287 const int ipHe2p1P = 6; 00288 00289 /* level 3 */ 00290 const int ipHe3s3S = 7; 00291 const int ipHe3s1S = 8; 00292 const int ipHe3p3P = 9; 00293 const int ipHe3d3D = 10; 00294 const int ipHe3d1D = 11; 00295 const int ipHe3p1P = 12; 00296 00297 /* level 4 */ 00298 const int ipHe4s3S = 13; 00299 const int ipHe4s1S = 14; 00300 const int ipHe4p3P = 15; 00301 const int ipHe4d3D = 16; 00302 const int ipHe4d1D = 17; 00303 const int ipHe4f3F = 18; 00304 const int ipHe4f1F = 19; 00305 const int ipHe4p1P = 20; 00306 00310 const int ipH_LIKE = 0; 00311 const int ipHE_LIKE = 1; 00312 const int ipLI_LIKE = 2; 00313 const int ipBE_LIKE = 3; 00314 const int ipB_LIKE = 4; 00315 const int ipC_LIKE = 5; 00316 const int ipN_LIKE = 6; 00317 const int ipO_LIKE = 7; 00318 const int ipF_LIKE = 8; 00319 00321 const int ipHYDROGEN = 0; 00322 const int ipHELIUM = 1; 00323 const int ipLITHIUM = 2; 00324 const int ipBERYLLIUM = 3; 00325 const int ipBORON = 4; 00326 const int ipCARBON = 5; 00327 const int ipNITROGEN = 6; 00328 const int ipOXYGEN = 7; 00329 const int ipFLUORINE = 8; 00330 const int ipNEON = 9; 00331 const int ipSODIUM = 10; 00332 const int ipMAGNESIUM = 11; 00333 const int ipALUMINIUM = 12; 00334 const int ipSILICON = 13; 00335 const int ipPHOSPHORUS = 14; 00336 const int ipSULPHUR = 15; 00337 const int ipCHLORINE = 16; 00338 const int ipARGON = 17; 00339 const int ipPOTASSIUM = 18; 00340 const int ipCALCIUM = 19; 00341 const int ipSCANDIUM = 20; 00342 const int ipTITANIUM = 21; 00343 const int ipVANADIUM = 22; 00344 const int ipCHROMIUM = 23; 00345 const int ipMANGANESE = 24; 00346 const int ipIRON = 25; 00347 const int ipCOBALT = 26; 00348 const int ipNICKEL = 27; 00349 const int ipCOPPER = 28; 00350 const int ipZINC = 29; 00351 00352 /*************************************************************************** 00353 * the following are prototypes for some routines that are part of the 00354 * debugging process - they come and go in any particular sub. 00355 * it is not necessary to declare them when used since they are defined here 00356 **************************************************************************/ 00357 00366 double fudge(long int ipnt); 00367 00371 void broken(void); 00372 00375 void fixit(void); 00376 00378 void CodeReview(void); 00379 00381 void TestCode(void); 00382 00388 void *MyMalloc(size_t size, const char *file, int line); 00389 00394 void *MyCalloc(size_t num, size_t size); 00395 00400 void *MyRealloc(void *p, size_t size); 00401 00406 void MyAssert(const char *file, int line); 00407 00411 NORETURN void cdExit(int iexit); 00412 00413 class cloudy_exit 00414 { 00415 const char* p_routine; 00416 const char* p_file; 00417 long p_line; 00418 int p_exit; 00419 public: 00420 cloudy_exit(const char* routine, const char* file, long line, int exit_code) 00421 { 00422 p_routine = routine; 00423 p_file = file; 00424 p_line = line; 00425 p_exit = exit_code; 00426 } 00427 virtual ~cloudy_exit() throw() 00428 { 00429 p_routine = NULL; 00430 p_file = NULL; 00431 } 00432 const char* routine() const throw() 00433 { 00434 return p_routine; 00435 } 00436 const char* file() const throw() 00437 { 00438 return p_file; 00439 } 00440 long line() const 00441 { 00442 return p_line; 00443 } 00444 int exit_status() const 00445 { 00446 return p_exit; 00447 } 00448 }; 00449 00450 // workarounds for __func__ are defined in cpu.h 00451 #define cdEXIT( FAIL ) throw cloudy_exit( __func__, __FILE__, __LINE__, FAIL ) 00452 00453 // calls like puts( "[Stop in MyRoutine]" ) have been integrated in cdEXIT above 00454 #define puts( STR ) Using_puts_before_cdEXIT_is_no_longer_needed 00455 00457 void ShowMe(void); 00458 00460 NORETURN void TotalInsanity(void); 00461 00463 NORETURN void BadRead(void); 00464 00466 NORETURN void BadOpen(void); 00467 00472 NORETURN void NoNumb(char *chCard); 00473 00477 int dbg_printf(int debug, const char *fmt, ...); 00478 00480 int dprintf(FILE *fp, const char *format, ...); 00481 00490 char *read_whole_line( char *chLine , int nChar , FILE *ioIN ); 00491 00492 /************************************************************************** 00493 * 00494 * various macros used by the code 00495 * 00496 **************************************************************************/ 00497 00501 #ifndef NDEBUG 00502 # define DEBUG 00503 #else 00504 # undef DEBUG 00505 #endif 00506 00514 #ifndef PARALLEL_MODE 00515 #define PARALLEL_MODE false 00516 #endif 00517 00519 #if defined(malloc) 00520 /* ...but if malloc is a macro, assume it is instrumented by a memory debugging tool 00521 * (e.g. dmalloc) */ 00522 # define MALLOC(exp) (malloc(exp)) 00523 #else 00524 /* Otherwise instrument and protect it ourselves */ 00525 # define MALLOC(exp) (MyMalloc(exp,__FILE__, __LINE__)) 00526 #endif 00527 00529 #if defined(calloc) 00530 /* ...but if calloc is a macro, assume it is instrumented by a memory debugging tool */ 00531 # define CALLOC calloc 00532 #else 00533 /* Otherwise instrument and protect it ourselves */ 00534 # define CALLOC MyCalloc 00535 #endif 00536 00538 #if defined(realloc) 00539 /* ...but if calloc is a macro, assume it is instrumented by a memory debugging tool */ 00540 # define REALLOC realloc 00541 #else 00542 /* Otherwise instrument and protect it ourselves */ 00543 # define REALLOC MyRealloc 00544 #endif 00545 00546 class bad_signal 00547 { 00548 int p_sig; 00549 public: 00550 explicit bad_signal(int sig) 00551 { 00552 p_sig = sig; 00553 } 00554 virtual ~bad_signal() throw() {} 00555 int sig() const throw() 00556 { 00557 return p_sig; 00558 } 00559 }; 00560 00561 class bad_assert 00562 { 00563 const char* p_file; 00564 long p_line; 00565 public: 00566 bad_assert(const char* file, long line) 00567 { 00568 p_file = file; 00569 p_line = line; 00570 } 00571 virtual ~bad_assert() throw() 00572 { 00573 p_file = NULL; 00574 } 00575 const char* file() const throw() 00576 { 00577 return p_file; 00578 } 00579 long line() const throw() 00580 { 00581 return p_line; 00582 } 00583 }; 00584 00585 class bad_mpi { 00586 long p_failMode; 00587 public: 00588 bad_mpi(long failMode) : p_failMode(failMode) {} 00589 long failMode() const { return p_failMode; } 00590 }; 00591 00593 #undef ASSERT 00594 #ifndef STD_ASSERT 00595 # ifdef NDEBUG 00596 # define ASSERT(exp) ((void)0) 00597 # elif defined ASSERTDEBUG 00598 # define ASSERT(exp) if (!(exp)) MyAssert(__FILE__, __LINE__) 00599 # else 00600 /* the do { ... } while ( 0 ) prevents bugs in code like this: 00601 * if( test ) 00602 * ASSERT( n == 10 ); 00603 * else 00604 * do something else... 00605 */ 00606 # define ASSERT(exp) \ 00607 do { \ 00608 if (!(exp)) \ 00609 { \ 00610 if( cpu.lgAssertAbort() ) \ 00611 abort(); \ 00612 else \ 00613 throw bad_assert(__FILE__,__LINE__); \ 00614 } \ 00615 } while( 0 ) 00616 # endif 00617 #else 00618 # define ASSERT(exp) (assert(exp)) 00619 #endif 00620 00621 NORETURN inline void OUT_OF_RANGE(const char* str) 00622 { 00623 if( cpu.lgAssertAbort() ) 00624 abort(); 00625 else 00626 throw out_of_range( str ); 00627 } 00628 00629 /* Windows does not define isnan */ 00630 /* use our version on all platforms since the isnanf 00631 * function does not exist under Solaris 9 either */ 00632 #undef isnan 00633 #define isnan MyIsnan 00634 00637 class t_debug : public Singleton<t_debug> 00638 { 00639 friend class Singleton<t_debug>; 00640 FILE *p_fp; 00641 int p_callLevel; 00642 protected: 00643 t_debug() 00644 { 00645 p_fp = stderr; 00646 p_callLevel = 0; 00647 } 00648 public: 00649 void enter(const char *name) 00650 { 00651 ++p_callLevel; 00652 fprintf(p_fp,"%*c%s\n",p_callLevel,'>',name); 00653 } 00654 void leave(const char *name) 00655 { 00656 fprintf(p_fp,"%*c%s\n",p_callLevel,'<',name); 00657 --p_callLevel; 00658 } 00659 }; 00660 00661 template<bool lgTrace> 00662 class debugtrace 00663 { 00664 const char *p_name; 00665 public: 00666 explicit debugtrace(const char *funcname) 00667 { 00668 p_name = funcname; 00669 # ifdef _MSC_VER 00670 /* disable warning that conditional expression is constant, true or false in if */ 00671 # pragma warning( disable : 4127 ) 00672 # endif 00673 if( lgTrace ) 00674 t_debug::Inst().enter(p_name); 00675 } 00676 ~debugtrace() 00677 { 00678 # ifdef _MSC_VER 00679 /* disable warning that conditional expression is constant, true or false in if */ 00680 # pragma warning( disable : 4127 ) 00681 # endif 00682 if( lgTrace ) 00683 t_debug::Inst().leave(p_name); 00684 p_name = NULL; 00685 } 00686 const char* name() const 00687 { 00688 return p_name; 00689 } 00690 }; 00691 00692 #ifdef DEBUG_FUN 00693 #define DEBUG_ENTRY( funcname ) debugtrace<true> DEBUG_ENTRY( funcname ) 00694 #else 00695 #ifdef __GNUC__ 00696 #define DEBUG_ENTRY( funcname ) ((void)0) 00697 #else 00698 #define DEBUG_ENTRY( funcname ) debugtrace<false> DEBUG_ENTRY( funcname ) 00699 #endif 00700 #endif 00701 00702 /* TorF(l) returns a 'T' or 'F' depending on the 'logical' expr 'l' */ 00703 inline char TorF( bool l ) { return l ? 'T' : 'F'; } 00704 /* */ 00705 00707 inline bool is_odd( int j ) { return (j&1) == 1; } 00708 inline bool is_odd( long j ) { return (j&1L) == 1L; } 00709 /* */ 00710 00712 inline long nint( double x ) { return static_cast<long>( (x < 0.) ? x-0.5 : x+0.5 ); } 00713 /* */ 00714 00715 /* define min for mixed arguments, the rest already exists */ 00716 inline long min( int a, long b ) { long c = a; return ( (c < b) ? c : b ); } 00717 inline long min( long a, int b ) { long c = b; return ( (a < c) ? a : c ); } 00718 inline double min( sys_float a, double b ) { double c = a; return ( (c < b) ? c : b ); } 00719 inline double min( double a, sys_float b ) { double c = b; return ( (a < c) ? a : c ); } 00720 00721 #undef MIN2 00722 00723 #define MIN2 min 00724 /* */ 00725 00726 #undef MIN3 00727 00728 #define MIN3(a,b,c) (min(min(a,b),c)) 00729 /* */ 00730 00731 #undef MIN4 00732 00733 #define MIN4(a,b,c,d) (min(min(a,b),min(c,d))) 00734 /* */ 00735 00736 /* define max for mixed arguments, the rest already exists */ 00737 inline long max( int a, long b ) { long c = a; return ( (c > b) ? c : b ); } 00738 inline long max( long a, int b ) { long c = b; return ( (a > c) ? a : c ); } 00739 inline double max( sys_float a, double b ) { double c = a; return ( (c > b) ? c : b ); } 00740 inline double max( double a, sys_float b ) { double c = b; return ( (a > c) ? a : c ); } 00741 00742 #undef MAX2 00743 00744 #define MAX2 max 00745 /* */ 00746 00747 #undef MAX3 00748 00749 #define MAX3(a,b,c) (max(max(a,b),c)) 00750 /* */ 00751 00752 #undef MAX4 00753 00754 #define MAX4(a,b,c,d) (max(max(a,b),max(c,d))) 00755 /* */ 00756 00761 template<class T> 00762 inline T sign( T x, T y ) 00763 { 00764 return ( y < T() ) ? -abs(x) : abs(x); 00765 } 00766 /* */ 00767 00769 template<class T> 00770 inline int sign3( T x ) { return ( x < T() ) ? -1 : ( ( x > T() ) ? 1 : 0 ); } 00771 /* */ 00772 00774 inline bool fp_equal( sys_float x, sys_float y, int n=3 ) 00775 { 00776 # ifdef _MSC_VER 00777 /* disable warning that conditional expression is constant, true or false in if */ 00778 # pragma warning( disable : 4127 ) 00779 # endif 00780 ASSERT( n >= 1 ); 00781 // mimic IEEE behavior 00782 if( isnan(x) || isnan(y) ) 00783 return false; 00784 int sx = sign3(x); 00785 int sy = sign3(y); 00786 // treat zero cases first to avoid division by zero below 00787 if( sx == 0 && sy == 0 ) 00788 return true; 00789 // either x or y is zero (but not both), or x and y have different sign 00790 if( sx*sy != 1 ) 00791 return false; 00792 x = abs(x); 00793 y = abs(y); 00794 return ( 1.f - min(x,y)/max(x,y) < ((sys_float)n+0.1f)*FLT_EPSILON ); 00795 } 00796 00797 inline bool fp_equal( double x, double y, int n=3 ) 00798 { 00799 ASSERT( n >= 1 ); 00800 // mimic IEEE behavior 00801 if( isnan(x) || isnan(y) ) 00802 return false; 00803 int sx = sign3(x); 00804 int sy = sign3(y); 00805 // treat zero cases first to avoid division by zero below 00806 if( sx == 0 && sy == 0 ) 00807 return true; 00808 // either x or y is zero (but not both), or x and y have different sign 00809 if( sx*sy != 1 ) 00810 return false; 00811 x = abs(x); 00812 y = abs(y); 00813 return ( 1. - min(x,y)/max(x,y) < ((double)n+0.1)*DBL_EPSILON ); 00814 } 00815 00816 inline bool fp_equal_tol( sys_float x, sys_float y, sys_float tol ) 00817 { 00818 ASSERT( tol > 0.f ); 00819 // mimic IEEE behavior 00820 if( isnan(x) || isnan(y) ) 00821 return false; 00822 // make sure the tolerance is not too stringent 00823 ASSERT( tol >= FLT_EPSILON*max(abs(x),abs(y)) ); 00824 return ( abs( x-y ) <= tol ); 00825 } 00826 00827 inline bool fp_equal_tol( double x, double y, double tol ) 00828 { 00829 ASSERT( tol > 0. ); 00830 // mimic IEEE behavior 00831 if( isnan(x) || isnan(y) ) 00832 return false; 00833 // make sure the tolerance is not too stringent 00834 ASSERT( tol >= DBL_EPSILON*max(abs(x),abs(y)) ); 00835 return ( abs( x-y ) <= tol ); 00836 } 00837 00838 #undef POW2 00839 00840 #define POW2 pow2 00841 template<class T> 00842 inline T pow2(T a) { return a*a; } 00843 /* */ 00844 00845 #undef POW3 00846 00847 #define POW3 pow3 00848 template<class T> 00849 inline T pow3(T a) { return a*a*a; } 00850 /* */ 00851 00852 #undef POW4 00853 00854 #define POW4(X) (pow2(X)*pow2(X)) 00855 /* */ 00856 00857 #undef SDIV 00858 00861 inline sys_float SDIV( sys_float x ) { return ( fabs((double)x) < (double)SMALLFLOAT ) ? (sys_float)SMALLFLOAT : x; } 00862 /* \todo should we use SMALLDOUBLE here ? it produces overflows now... PvH */ 00863 inline double SDIV( double x ) { return ( fabs(x) < (double)SMALLFLOAT ) ? (double)SMALLFLOAT : x; } 00864 // inline double SDIV( double x ) { return ( fabs(x) < SMALLDOUBLE ) ? SMALLDOUBLE : x; } 00865 /* */ 00866 00870 inline sys_float safe_div(sys_float x, sys_float y) 00871 { 00872 // this should crash... 00873 if( isnan(x) || isnan(y) ) 00874 return x/y; 00875 int sx = sign3(x); 00876 int sy = sign3(y); 00877 // 0/0 -> NaN, this should crash as well... 00878 if( sx == 0 && sy == 0 ) 00879 return x/y; 00880 if( sx == 0 ) 00881 return 0.; 00882 if( sy == 0 ) 00883 return ( sx < 0 ) ? -FLT_MAX : FLT_MAX; 00884 // at this stage x != 0. and y != 0. 00885 sys_float ay = abs(y); 00886 if( ay >= 1.f ) 00887 return x/y; 00888 else 00889 { 00890 // multiplication is safe since ay < 1. 00891 if( abs(x) < ay*FLT_MAX ) 00892 return x/y; 00893 else 00894 return ( sx*sy < 0 ) ? -FLT_MAX : FLT_MAX; 00895 } 00896 } 00897 00901 inline double safe_div(double x, double y) 00902 { 00903 // this should crash... 00904 if( isnan(x) || isnan(y) ) 00905 return x/y; 00906 int sx = sign3(x); 00907 int sy = sign3(y); 00908 // 0/0 -> NaN, this should crash as well... 00909 if( sx == 0 && sy == 0 ) 00910 return x/y; 00911 if( sx == 0 ) 00912 return 0.; 00913 if( sy == 0 ) 00914 return ( sx < 0 ) ? -DBL_MAX : DBL_MAX; 00915 // at this stage x != 0. and y != 0. 00916 double ay = abs(y); 00917 if( ay >= 1. ) 00918 return x/y; 00919 else 00920 { 00921 // multiplication is safe since ay < 1. 00922 if( abs(x) < ay*DBL_MAX ) 00923 return x/y; 00924 else 00925 return ( sx*sy < 0 ) ? -DBL_MAX : DBL_MAX; 00926 } 00927 } 00928 00929 #undef HMRATE 00930 /*HMRATE compile molecular rates using Hollenbach and McKee fits */ 00931 /* #define HMRATE(a,b,c) ( ((b) == 0 && (c) == 0) ? (a) : \ 00932 * ( ((c) == 0) ? (a)*pow(phycon.te/300.,(b)) : \ 00933 * ( ((c)/phycon.te > 50.) ? 0. : ( ((b) == 0) ? (a)*exp(-(c)/phycon.te) : \ 00934 * (a)*pow(phycon.te/300.,(b))*exp(-(c)/phycon.te) ) ) ) ) */ 00935 #define HMRATE(a,b,c) hmrate4(a,b,c,phycon.te) 00936 00937 inline double hmrate4( double a, double b, double c, double te ) 00938 { 00939 if( b == 0. && c == 0. ) 00940 return a; 00941 else if( c == 0. ) 00942 return a*pow(te/300.,b); 00943 else if( b == 0. ) 00944 return ( c/te <= 50. ) ? a*exp(-c/te) : 0.; 00945 else 00946 return ( c/te <= 50. ) ? a*pow(te/300.,b)*exp(-c/te) : 0.; 00947 } 00948 00949 template<class T> 00950 inline void invalidate_array(T* p, size_t size) 00951 { 00952 memset( p, -1, size ); 00953 } 00954 00955 inline void invalidate_array(double* p, size_t size) 00956 { 00957 set_NaN( p, (long)(size/sizeof(double)) ); 00958 } 00959 00960 inline void invalidate_array(sys_float* p, size_t size) 00961 { 00962 set_NaN( p, (long)(size/sizeof(sys_float)) ); 00963 } 00964 00990 template<class T> 00991 class auto_vec 00992 { 00993 T* ptr; 00994 00995 template<class U> 00996 struct auto_vec_ref 00997 { 00998 U* ptr; 00999 01000 explicit auto_vec_ref( U* p ) 01001 { 01002 ptr = p; 01003 } 01004 }; 01005 01006 public: 01007 typedef T element_type; 01008 01009 // 20.4.5.1 construct/copy/destroy: 01010 01011 explicit auto_vec( element_type* p = NULL ) throw() 01012 { 01013 ptr = p; 01014 } 01015 auto_vec( auto_vec& p ) throw() 01016 { 01017 ptr = p.release(); 01018 } 01019 auto_vec& operator= ( auto_vec& p ) throw() 01020 { 01021 reset( p.release() ); 01022 return *this; 01023 } 01024 ~auto_vec() throw() 01025 { 01026 delete[] ptr; 01027 } 01028 01029 // 20.4.5.2 members: 01030 01031 element_type& operator[] ( ptrdiff_t n ) const throw() 01032 { 01033 return *(ptr+n); 01034 } 01035 element_type* get() const throw() 01036 { 01037 return ptr; 01038 } 01039 // for consistency with other container classes 01040 element_type* data() const throw() 01041 { 01042 return ptr; 01043 } 01044 element_type* release() throw() 01045 { 01046 element_type* p = ptr; 01047 ptr = NULL; 01048 return p; 01049 } 01050 void reset( element_type* p = NULL ) throw() 01051 { 01052 if( p != ptr ) 01053 { 01054 delete[] ptr; 01055 ptr = p; 01056 } 01057 } 01058 01059 // 20.4.5.3 conversions: 01060 01061 auto_vec( auto_vec_ref<element_type> r ) throw() 01062 { 01063 ptr = r.ptr; 01064 } 01065 auto_vec& operator= ( auto_vec_ref<element_type> r ) throw() 01066 { 01067 if( r.ptr != ptr ) 01068 { 01069 delete[] ptr; 01070 ptr = r.ptr; 01071 } 01072 return *this; 01073 } 01074 operator auto_vec_ref<element_type>() throw() 01075 { 01076 return auto_vec_ref<element_type>( this->release() ); 01077 } 01078 }; 01079 01080 #include "container_classes.h" 01081 01082 /*Many structure were intoduced by Humeshkar B Nemala as a part of his Thesis 01083 *The structures were designed to read in transition,radiative and collisional data 01084 *from two major databases:LEIDEN and CHIANTI 01085 01086 * these structures define the emission, collision, state, and transition classes*/ 01087 01088 typedef struct t_transition transition; 01089 typedef struct t_quantumState quantumState; 01090 typedef struct t_emission emission; 01091 typedef struct t_collision collision; 01092 typedef struct t_species species; 01093 01094 struct t_emission 01095 { 01104 int iRedisFun; 01105 01107 long int ipFine; 01108 01122 realnum TauIn; 01123 01134 realnum TauTot; 01135 01141 realnum TauCon; 01142 01144 realnum FracInwd; 01145 01148 double pump; 01149 01151 double xIntensity; 01152 01154 double phots; 01155 01157 realnum gf; 01158 01160 realnum Pesc; 01161 01163 realnum Pelec_esc; 01164 01166 realnum Pdest; 01167 01171 realnum dampXvel; 01172 01174 realnum damp; 01175 01177 double ColOvTot; 01178 01184 realnum opacity; 01185 01187 double PopOpc; 01188 01190 realnum Aul; 01191 01193 double ots; 01194 01195 /* this points back to the transition that points here. */ 01196 transition *tran; 01197 01198 /* linked-list */ 01199 emission *next; 01200 01201 }; 01202 01203 /*The species structure is used to hold information about a particular atom,ion or molecule 01204 mentioned in the species.ini file.The name of the atom/ion/molecule is used to obtain the density 01205 of molecules in the case of the Leiden Database and along with atomic number and ion stage the 01206 density of atoms/ions in the case of the CHIANTI database */ 01207 struct t_species 01208 { 01209 /*Name of the atom/ion/ molecule*/ 01210 char *chptrSpName; 01211 /*Atomic Number*/ 01212 int intAtNo; 01213 /*Ionization Stage*/ 01214 int intIonStage; 01215 /*Actual Number of energy levels in the data file*/ 01216 long numLevels_max; 01217 /*Number of energy levels used locally*/ 01218 long numLevels_local; 01219 /*Molecular weight*/ 01220 realnum fmolweight; 01221 01222 01223 }; 01224 01225 /*This structure is specifically used to hold the collision data in the format given in the LEIDEN Database 01226 The data is available as collsion rate coefficients(cm3 s-1) over different temperatures*/ 01227 typedef struct t_CollRatesArray 01228 { 01229 /*Number of temps*/ 01230 long ntemps; 01231 /*Array of temps*/ 01232 double *temps; 01233 /*Matrix of collision rates(temp,up,lo)*/ 01234 double ***collrates; 01235 01236 } CollRateCoeffArray ; 01237 01238 /*This structure is specifically used to hold the collision data in the format given in the CHIANTI Database 01239 The data is available as spline fits to the maxwellian averaged collision strengths */ 01240 typedef struct t_CollSplinesArray 01241 { 01242 /*Matrix of spline fits(hi,lo,spline index)* 01243 *The first five columns gives the no of spline pts,transition type,gf value,delta E 01244 *& Scaling parameter ,in the specified order*/ 01245 /*The transition type basically tells how the temperature and collision 01246 strengths have been scaled*/ 01247 double *collspline; 01248 double *SplineSecDer; 01249 01250 long nSplinePts; 01251 long intTranType; 01252 double gf; 01253 double EnergyDiff; 01254 double ScalingParam; 01255 01256 } CollSplinesArray ; 01257 01258 struct t_collision 01259 { 01260 /* species *collider;*/ 01261 /* Linked list for possible colliders */ 01262 /* collision *next; */ 01263 01265 realnum ColUL; 01266 01268 realnum cs; 01269 01271 realnum csi[ipNCOLLIDER]; 01272 01274 double cool , heat; 01275 01276 }; 01277 01278 /*Generalized structure used to hold the energy level information for both atoms/ions and molecules*/ 01279 struct t_quantumState 01280 { 01281 char chLabel[5]; 01282 char chConfig[11]; 01283 01284 species *sp; 01285 01287 realnum energy; 01288 01290 realnum g; 01291 01293 double Pop; 01294 01296 int IonStg; 01298 int nelem; 01299 01301 double ConBoltz; 01302 01303 /* S is multiplicity. */ 01304 long n, l, S, j; 01305 01307 double lifetime; 01308 01309 /* linked-list */ 01310 quantumState *next; 01311 }; 01312 01313 /*Generalized structure used to hold the transition information for both atoms/ions and molecules*/ 01314 struct t_transition 01315 { 01316 quantumState *Lo, *Hi; 01317 emission *Emis; 01318 collision Coll; 01319 01320 /* Possible linked-list optimization */ 01321 /*transition *nextemis, *nextcoll; */ 01322 01324 realnum WLAng; 01325 01327 realnum EnergyK; 01328 01330 realnum EnergyErg; 01331 01333 realnum EnergyWN; 01334 01339 long ipCont; 01340 }; 01341 01342 /* Explicit instantiations for debugging purposes */ 01343 INSTANTIATE_MULTI_ARR( quantumState, MEM_LAYOUT_VAL, lgBOUNDSCHECKVAL ); 01344 INSTANTIATE_MULTI_ARR( transition, MEM_LAYOUT_VAL, lgBOUNDSCHECKVAL ); 01345 01346 01347 /*************************************************************************** 01348 * 01349 * a series of Cloudy service routines, used throughout code, 01350 * 01351 **************************************************************************/ 01352 01359 typedef enum { SPM_RELAX, SPM_KEEP_EMPTY, SPM_STRICT } split_mode; 01360 01362 void Split(const string& str, // input string 01363 const string& sep, // separator, may be multiple characters 01364 vector<string>& lst, // the separated items will be appended here 01365 split_mode mode); // see above 01366 01369 inline bool FindAndReplace(string& str, 01370 const string& substr, 01371 const string& newstr) 01372 { 01373 string::size_type ptr = str.find( substr ); 01374 if( ptr != string::npos ) 01375 str.replace( ptr, substr.length(), newstr ); 01376 return ptr != string::npos; 01377 } 01378 01381 inline bool FindAndErase(string& str, 01382 const string& substr) 01383 { 01384 return FindAndReplace( str, substr, "" ); 01385 } 01386 01392 double csphot(long int inu, long int ithr, long int iofset); 01393 01398 double RandGauss(double xMean, double s ); 01399 01403 double MyGaussRand( double PctUncertainty ); 01404 01406 double AnuUnit(realnum energy); 01407 01412 void cap4(char *chCAP , char *chLab); 01413 01416 void uncaps(char *chCard ); 01417 01420 void caps(char *chCard ); 01421 01424 double e2( 01425 double t ); 01426 01429 double ee1(double x); 01430 01434 double ee1_safe(double x); 01435 01442 double FFmtRead(const char *chCard, 01443 long int *ipnt, 01444 long int last, 01445 bool *lgEOL); 01446 01452 long nMatch(const char *chKey, 01453 const char *chCard); 01454 01460 long int GetElem( char *chCARD_CAPS ); 01461 01470 int GetQuote(char *chLabel, char *chCard, bool lgABORT ); 01471 01472 /* want to define this only if no native os support exists */ 01473 #if !HAVE_POWI 01474 01475 double powi( double , long int ); 01476 #endif 01477 01480 long int ipow( long, long ); 01481 01484 void PrintE82( FILE*, double ); 01485 01487 void PrintE71( FILE*, double ); 01488 01490 void PrintE93( FILE*, double ); 01491 01497 char *PrintEfmt(const char *fmt, double val ); 01498 01500 const double SEXP_LIMIT = 84.; 01502 sys_float sexp(sys_float x); 01503 double sexp(double x); 01504 01509 double dsexp(double x); 01510 01515 double plankf(long int ip); 01516 01523 double qg32( double, double, double(*)(double) ); 01524 /* declar of optimize_func, the last arg, changed from double(*)() to above, 01525 * seemed to fix flags that were raised */ 01526 01527 01537 void spsort( realnum x[], long int n, long int iperm[], int kflag, int *ier); 01538 01539 /*vfun approximate form of Voigt function */ 01540 inline double vfun(double damp, double x) 01541 { 01542 // constant is SQRTPI 01543 return sexp(x*x) + damp/1.772453850905516027298167/(1. + x*x); 01544 } 01545 01546 /************************************************************************** 01547 * 01548 * disable some bogus errors in the ms c compiler 01549 * 01550 **************************************************************************/ 01551 01552 /* */ 01553 #ifdef _MSC_VER 01554 /* disable warning that conditional expression is constant, true or false in if */ 01555 # pragma warning( disable : 4127 ) 01556 /* disable strcat warning */ 01557 # pragma warning( disable : 4996 ) 01558 /* disable bogus underflow warning in MS VS*/ 01559 # pragma warning( disable : 4056 ) 01560 /* disable "inline function removed since not used", MS VS*/ 01561 # pragma warning( disable : 4514 ) 01562 /* disable "assignment operator could not be generated", cddefines.h 01563 * line 126 */ 01564 # pragma warning( disable : 4512 ) 01565 #endif 01566 #ifdef __INTEL_COMPILER 01567 # pragma warning( disable : 1572 ) 01568 #endif 01569 /* */ 01570 01571 /*lint +e129 these resolve several issues pclint has with my system headers */ 01572 /*lint +e78 */ 01573 /*lint +e830 */ 01574 /*lint +e38 */ 01575 /*lint +e148 */ 01576 /*lint +e114 */ 01577 /*lint +e18 */ 01578 /*lint +e49 */ 01579 01580 #endif /* _CDDEFINES_H_ */ 01581