00001
00002
00003
00004 #ifndef CDDEFINES_H_
00005 #define CDDEFINES_H_
00006
00007 #ifdef _MSC_VER
00008
00009 # pragma warning( disable : 4127 )
00010
00011 # ifndef WIN32_LEAN_AND_MEAN
00012 # define WIN32_LEAN_AND_MEAN
00013 # endif
00014 #endif
00015
00016 #ifdef __clang__
00017
00018 #pragma clang diagnostic ignored "-Wmismatched-tags"
00019 #endif
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include <cstdio>
00032 #include <cstdlib>
00033 #include <cctype>
00034
00035
00036 #define _USE_MATH_DEFINES
00037 #include <cmath>
00038 #include <cassert>
00039 #include <cstring>
00040 #include <cfloat>
00041 #include <climits>
00042 #include <ctime>
00043 #if defined(__sun) && defined(__SUNPRO_CC)
00044
00045 #include <signal.h>
00046 #else
00047 #include <csignal>
00048 #endif
00049
00050 #include <limits>
00051 #include <string>
00052 #include <sstream>
00053 #include <iomanip>
00054 #include <vector>
00055 #include <valarray>
00056 #include <complex>
00057 #include <map>
00058 #include <memory>
00059 #include <stdexcept>
00060 #include <algorithm>
00061 #include <fstream>
00062 #include <bitset>
00063 #ifdef DMALLOC
00064 #include <dmalloc.h>
00065 #endif
00066
00067
00068 #if defined(_MSC_VER) && !defined(SYS_CONFIG)
00069 #define SYS_CONFIG "cloudyconfig_vs.h"
00070 #endif
00071
00072
00073 #ifdef SYS_CONFIG
00074 #include SYS_CONFIG
00075 #else
00076 #include "cloudyconfig.h"
00077 #endif
00078
00079 #ifdef MPI_GRID_RUN
00080 #define MPI_ENABLED
00081 #endif
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091 using namespace std;
00092
00094 #ifndef EXTERN
00095 #define EXTERN extern
00096 #endif
00097
00098 #undef STATIC
00099 #ifdef USE_GPROF
00100 #define STATIC
00101 #else
00102 #define STATIC static
00103 #endif
00104
00105 #ifdef FLT_IS_DBL
00106 typedef double realnum;
00107 #else
00108 typedef float realnum;
00109 #endif
00110
00111 typedef float sys_float;
00112
00113 #define float PLEASE_USE_REALNUM_NOT_FLOAT
00114
00115
00116 template<bool> struct StaticAssertFailed;
00117 template<> struct StaticAssertFailed<true> {};
00118 #define STATIC_ASSERT(x) ((void)StaticAssertFailed< (x) == true >())
00119
00120
00121 #include "cpu.h"
00122
00123
00124
00142
00143
00144
00145
00146
00147 template<typename T> class Singleton
00148 {
00149 public:
00150 static T& Inst()
00151 {
00152 static T instance;
00153 return instance;
00154 }
00155 };
00156
00157
00158
00159
00160
00161
00162
00163
00167 EXTERN FILE *ioQQQ;
00168
00169 EXTERN FILE *ioStdin;
00170
00171 extern FILE *ioMAP;
00172
00175 EXTERN FILE* ioPrnErr;
00176
00178 EXTERN bool lgAbort;
00179
00182 EXTERN bool lgTestCodeCalled;
00183
00186 EXTERN bool lgTestCodeEnabled;
00187
00190 EXTERN bool lgPrnErr;
00191
00194 EXTERN long int nzone;
00195
00197 EXTERN double fnzone;
00198
00201 EXTERN long int iteration;
00202
00209 extern const double ZeroNum;
00210
00211
00212
00213
00214
00215
00216
00221 const int FILENAME_PATH_LENGTH = 200;
00222
00224 const int FILENAME_PATH_LENGTH_2 = FILENAME_PATH_LENGTH*2;
00225
00229 const int INPUT_LINE_LENGTH = 2000;
00230
00233 const int LIMELM = 30;
00234
00236 const int NISO = 2;
00237
00241 const int NHYDRO_MAX_LEVEL = 401;
00242
00244 const int N_H_MOLEC = 8;
00245
00247 const double MAX_DENSITY = 1.e24;
00248
00250 const double DEPTH_OFFSET = 1.e-30;
00251
00252
00253 enum collider {
00254 ipELECTRON,
00255 ipPROTON,
00256 ipHE_PLUS,
00257 ipALPHA,
00258
00259
00260
00261
00262 ipNCOLLIDER
00263 };
00264
00265
00266
00267 const int ipRecEsc = 2;
00268
00269 const int ipRecNetEsc = 1;
00270
00271 const int ipRecRad = 0;
00276
00277
00278 const int ipPRD = 1;
00279
00280 const int ipCRD = -1;
00281
00282 const int ipCRDW = 2;
00283
00284 const int ipLY_A = -2;
00285
00286 const int ipDEST_K2 = 1;
00287
00288 const int ipDEST_INCOM = 2;
00289
00290 const int ipDEST_SIMPL = 3;
00291
00294 const int ipH1s = 0;
00295 const int ipH2s = 1;
00296 const int ipH2p = 2;
00297 const int ipH3s = 3;
00298 const int ipH3p = 4;
00299 const int ipH3d = 5;
00300 const int ipH4s = 6;
00301 const int ipH4p = 7;
00302 const int ipH4d = 8;
00303 const int ipH4f = 9;
00304
00307
00308 const int ipHe1s1S = 0;
00309
00310
00311 const int ipHe2s3S = 1;
00312 const int ipHe2s1S = 2;
00313 const int ipHe2p3P0 = 3;
00314 const int ipHe2p3P1 = 4;
00315 const int ipHe2p3P2 = 5;
00316 const int ipHe2p1P = 6;
00317
00318
00319 const int ipHe3s3S = 7;
00320 const int ipHe3s1S = 8;
00321 const int ipHe3p3P = 9;
00322 const int ipHe3d3D = 10;
00323 const int ipHe3d1D = 11;
00324 const int ipHe3p1P = 12;
00325
00326
00327 const int ipHe4s3S = 13;
00328 const int ipHe4s1S = 14;
00329 const int ipHe4p3P = 15;
00330 const int ipHe4d3D = 16;
00331 const int ipHe4d1D = 17;
00332 const int ipHe4f3F = 18;
00333 const int ipHe4f1F = 19;
00334 const int ipHe4p1P = 20;
00335
00339 const int ipH_LIKE = 0;
00340 const int ipHE_LIKE = 1;
00341 const int ipLI_LIKE = 2;
00342 const int ipBE_LIKE = 3;
00343 const int ipB_LIKE = 4;
00344 const int ipC_LIKE = 5;
00345 const int ipN_LIKE = 6;
00346 const int ipO_LIKE = 7;
00347 const int ipF_LIKE = 8;
00348 const int ipNE_LIKE = 9;
00349 const int ipNA_LIKE = 10;
00350 const int ipMG_LIKE = 11;
00351 const int ipAL_LIKE = 12;
00352 const int ipSI_LIKE = 13;
00353 const int ipP_LIKE = 14;
00354 const int ipS_LIKE = 15;
00355 const int ipCL_LIKE = 16;
00356 const int ipAR_LIKE = 17;
00357
00359 const int ipHYDROGEN = 0;
00360 const int ipHELIUM = 1;
00361 const int ipLITHIUM = 2;
00362 const int ipBERYLLIUM = 3;
00363 const int ipBORON = 4;
00364 const int ipCARBON = 5;
00365 const int ipNITROGEN = 6;
00366 const int ipOXYGEN = 7;
00367 const int ipFLUORINE = 8;
00368 const int ipNEON = 9;
00369 const int ipSODIUM = 10;
00370 const int ipMAGNESIUM = 11;
00371 const int ipALUMINIUM = 12;
00372 const int ipSILICON = 13;
00373 const int ipPHOSPHORUS = 14;
00374 const int ipSULPHUR = 15;
00375 const int ipCHLORINE = 16;
00376 const int ipARGON = 17;
00377 const int ipPOTASSIUM = 18;
00378 const int ipCALCIUM = 19;
00379 const int ipSCANDIUM = 20;
00380 const int ipTITANIUM = 21;
00381 const int ipVANADIUM = 22;
00382 const int ipCHROMIUM = 23;
00383 const int ipMANGANESE = 24;
00384 const int ipIRON = 25;
00385 const int ipCOBALT = 26;
00386 const int ipNICKEL = 27;
00387 const int ipCOPPER = 28;
00388 const int ipZINC = 29;
00389
00390
00391
00392
00393
00394
00395
00404 double fudge(long int ipnt);
00405
00409 void broken(void);
00410
00413 void fixit(void);
00414
00416 void CodeReview(void);
00417
00419 void TestCode(void);
00420
00426 void *MyMalloc(size_t size, const char *file, int line);
00427
00432 void *MyCalloc(size_t num, size_t size);
00433
00438 void *MyRealloc(void *p, size_t size);
00439
00444 void MyAssert(const char *file, int line, const char *comment);
00445
00447 void cdPrepareExit();
00448
00449 class cloudy_exit
00450 {
00451 const char* p_routine;
00452 const char* p_file;
00453 long p_line;
00454 int p_exit;
00455 public:
00456 cloudy_exit(const char* routine, const char* file, long line, int exit_code)
00457 {
00458 p_routine = routine;
00459 p_file = file;
00460 p_line = line;
00461 p_exit = exit_code;
00462 }
00463 virtual ~cloudy_exit() throw()
00464 {
00465 p_routine = NULL;
00466 p_file = NULL;
00467 }
00468 const char* routine() const throw()
00469 {
00470 return p_routine;
00471 }
00472 const char* file() const throw()
00473 {
00474 return p_file;
00475 }
00476 long line() const
00477 {
00478 return p_line;
00479 }
00480 int exit_status() const
00481 {
00482 return p_exit;
00483 }
00484 };
00485
00486
00487 #define cdEXIT( FAIL ) throw cloudy_exit( __func__, __FILE__, __LINE__, FAIL )
00488
00489
00490 #define puts( STR ) Using_puts_before_cdEXIT_is_no_longer_needed
00491
00493 void ShowMe(void);
00494
00496 NORETURN void TotalInsanity(void);
00497
00498
00499
00500
00501 template<class T>
00502 T TotalInsanityAsStub()
00503 {
00504
00505 if( ZeroNum == 0. )
00506 TotalInsanity();
00507 else
00508 return T();
00509 }
00510
00512 NORETURN void BadRead(void);
00513
00517 int dbg_printf(int debug, const char *fmt, ...);
00518
00520 int dprintf(FILE *fp, const char *format, ...);
00521
00530 char *read_whole_line( char *chLine , int nChar , FILE *ioIN );
00531
00532
00533
00534
00535
00536
00537
00541 #ifndef NDEBUG
00542 # define DEBUG
00543 #else
00544 # undef DEBUG
00545 #endif
00546
00548 #if defined(malloc)
00549
00550
00551 # define MALLOC(exp) (malloc(exp))
00552 #else
00553
00554 # define MALLOC(exp) (MyMalloc(exp,__FILE__, __LINE__))
00555 #endif
00556
00558 #if defined(calloc)
00559
00560 # define CALLOC calloc
00561 #else
00562
00563 # define CALLOC MyCalloc
00564 #endif
00565
00567 #if defined(realloc)
00568
00569 # define REALLOC realloc
00570 #else
00571
00572 # define REALLOC MyRealloc
00573 #endif
00574
00575 class bad_signal
00576 {
00577 int p_sig;
00578 public:
00579 explicit bad_signal(int sig)
00580 {
00581 p_sig = sig;
00582 }
00583 virtual ~bad_signal() throw() {}
00584 int sig() const throw()
00585 {
00586 return p_sig;
00587 }
00588 };
00589
00590 class bad_assert
00591 {
00592 const char* p_file;
00593 long p_line;
00594 const char* p_comment;
00595 public:
00596 bad_assert(const char* file, long line, const char* comment);
00597 void print(void) const
00598 {
00599 fprintf(ioQQQ,"DISASTER Assertion failure at %s:%ld\n%s\n",
00600 p_file, p_line, p_comment);
00601 }
00602 virtual ~bad_assert() throw()
00603 {
00604 p_file = NULL;
00605 }
00606 const char* file() const throw()
00607 {
00608 return p_file;
00609 }
00610 long line() const throw()
00611 {
00612 return p_line;
00613 }
00614 const char *comment() const throw()
00615 {
00616 return p_comment;
00617 }
00618 };
00619
00620
00621
00622
00623
00624
00625
00626 #undef ASSERT
00627 #ifndef OLD_ASSERT
00628 # if NDEBUG
00629 # define ASSERT(exp) ((void)0)
00630 # else
00631 # define ASSERT(exp) \
00632 do { \
00633 if (!(exp)) \
00634 { \
00635 bad_assert a(__FILE__,__LINE__,"Failed: " #exp); \
00636 if( cpu.lgAssertAbort() ) \
00637 { \
00638 a.print(); \
00639 abort(); \
00640 } \
00641 else \
00642 throw a; \
00643 } \
00644 } while( 0 )
00645 # endif
00646 #else
00647
00648 # ifdef NDEBUG
00649 # define ASSERT(exp) ((void)0)
00650 # else
00651 # define ASSERT(exp) \
00652 do { \
00653 if (!(exp)) \
00654 MyAssert(__FILE__, __LINE__, "Failed: " #exp); \
00655 } while( 0 )
00656 # endif
00657 #endif
00658
00659 #define MESSAGE_ASSERT(msg, exp) ASSERT( (msg) ? (exp) : false )
00660
00661 inline NORETURN void OUT_OF_RANGE(const char* str)
00662 {
00663 if( cpu.lgAssertAbort() )
00664 abort();
00665 else
00666 throw out_of_range( str );
00667 }
00668
00669
00670
00671
00672 #undef isnan
00673 #define isnan MyIsnan
00674
00677 class t_debug : public Singleton<t_debug>
00678 {
00679 friend class Singleton<t_debug>;
00680 FILE *p_fp;
00681 int p_callLevel;
00682 protected:
00683 t_debug()
00684 {
00685 p_fp = stderr;
00686 p_callLevel = 0;
00687 }
00688 public:
00689 void enter(const char *name)
00690 {
00691 ++p_callLevel;
00692 fprintf(p_fp,"%*c%s\n",p_callLevel,'>',name);
00693 }
00694 void leave(const char *name)
00695 {
00696 fprintf(p_fp,"%*c%s\n",p_callLevel,'<',name);
00697 --p_callLevel;
00698 }
00699 };
00700
00701 template<bool lgTrace>
00702 class debugtrace
00703 {
00704 const char *p_name;
00705 public:
00706 explicit debugtrace(const char *funcname)
00707 {
00708 p_name = funcname;
00709 # ifdef _MSC_VER
00710
00711 # pragma warning( disable : 4127 )
00712 # endif
00713 if( lgTrace )
00714 t_debug::Inst().enter(p_name);
00715 }
00716 ~debugtrace()
00717 {
00718 # ifdef _MSC_VER
00719
00720 # pragma warning( disable : 4127 )
00721 # endif
00722 if( lgTrace )
00723 t_debug::Inst().leave(p_name);
00724 p_name = NULL;
00725 }
00726 const char* name() const
00727 {
00728 return p_name;
00729 }
00730 };
00731
00732 #ifdef DEBUG_FUN
00733 #define DEBUG_ENTRY( funcname ) debugtrace<true> DEBUG_ENTRY( funcname )
00734 #else
00735 #ifdef HAVE_FUNC
00736 #define DEBUG_ENTRY( funcname ) ((void)0)
00737 #else
00738 #define DEBUG_ENTRY( funcname ) debugtrace<false> DEBUG_ENTRY( funcname )
00739 #endif
00740 #endif
00741
00742
00743 inline char tolower(char c)
00744 {
00745 return static_cast<char>( tolower( static_cast<int>(c) ) );
00746 }
00747 inline unsigned char tolower(unsigned char c)
00748 {
00749 return static_cast<unsigned char>( tolower( static_cast<int>(c) ) );
00750 }
00751
00752 inline char toupper(char c)
00753 {
00754 return static_cast<char>( toupper( static_cast<int>(c) ) );
00755 }
00756 inline unsigned char toupper(unsigned char c)
00757 {
00758 return static_cast<unsigned char>( toupper( static_cast<int>(c) ) );
00759 }
00760
00761
00762 inline char TorF( bool l ) { return l ? 'T' : 'F'; }
00763
00764
00766 inline bool is_odd( int j ) { return (j&1) == 1; }
00767 inline bool is_odd( long j ) { return (j&1L) == 1L; }
00768
00769
00771 inline long nint( double x ) { return static_cast<long>( (x < 0.) ? x-0.5 : x+0.5 ); }
00772
00773
00774
00775 inline long min( int a, long b ) { long c = a; return ( (c < b) ? c : b ); }
00776 inline long min( long a, int b ) { long c = b; return ( (a < c) ? a : c ); }
00777 inline double min( sys_float a, double b ) { double c = a; return ( (c < b) ? c : b ); }
00778 inline double min( double a, sys_float b ) { double c = b; return ( (a < c) ? a : c ); }
00779
00780
00781 #ifndef HAVE_POWI
00782
00783 double powi( double , long int );
00784 #endif
00785
00786
00787 #ifndef HAVE_POW_DOUBLE_INT
00788 inline double pow( double x, int i ) { return powi( x, long(i) ); }
00789 #endif
00790
00791 #ifndef HAVE_POW_DOUBLE_LONG
00792 inline double pow( double x, long i ) { return powi( x, i ); }
00793 #endif
00794
00795 #ifndef HAVE_POW_FLOAT_INT
00796 inline sys_float pow( sys_float x, int i ) { return sys_float( powi( double(x), long(i) ) ); }
00797 #endif
00798
00799 #ifndef HAVE_POW_FLOAT_LONG
00800 inline sys_float pow( sys_float x, long i ) { return sys_float( powi( double(x), i ) ); }
00801 #endif
00802
00803 #ifndef HAVE_POW_FLOAT_DOUBLE
00804 inline double pow( sys_float x, double y ) { return pow( double(x), y ); }
00805 #endif
00806
00807 #ifndef HAVE_POW_DOUBLE_FLOAT
00808 inline double pow( double x, sys_float y ) { return pow( x, double(y) ); }
00809 #endif
00810
00811 #undef MIN2
00812
00813 #define MIN2 min
00814
00815
00816 #undef MIN3
00817
00818 #define MIN3(a,b,c) (min(min(a,b),c))
00819
00820
00821 #undef MIN4
00822
00823 #define MIN4(a,b,c,d) (min(min(a,b),min(c,d)))
00824
00825
00826
00827 inline long max( int a, long b ) { long c = a; return ( (c > b) ? c : b ); }
00828 inline long max( long a, int b ) { long c = b; return ( (a > c) ? a : c ); }
00829 inline double max( sys_float a, double b ) { double c = a; return ( (c > b) ? c : b ); }
00830 inline double max( double a, sys_float b ) { double c = b; return ( (a > c) ? a : c ); }
00831
00832 #undef MAX2
00833
00834 #define MAX2 max
00835
00836
00837 #undef MAX3
00838
00839 #define MAX3(a,b,c) (max(max(a,b),c))
00840
00841
00842 #undef MAX4
00843
00844 #define MAX4(a,b,c,d) (max(max(a,b),max(c,d)))
00845
00846
00851 template<class T>
00852 inline T sign( T x, T y )
00853 {
00854 return ( y < T() ) ? -abs(x) : abs(x);
00855 }
00856
00857
00859 template<class T>
00860 inline int sign3( T x ) { return ( x < T() ) ? -1 : ( ( x > T() ) ? 1 : 0 ); }
00861
00862
00864 inline bool fp_equal( sys_float x, sys_float y, int n=3 )
00865 {
00866 #ifdef _MSC_VER
00867
00868 # pragma warning( disable : 4127 )
00869 #endif
00870 ASSERT( n >= 1 );
00871
00872 if( isnan(x) || isnan(y) )
00873 return false;
00874 int sx = sign3(x);
00875 int sy = sign3(y);
00876
00877 if( sx == 0 && sy == 0 )
00878 return true;
00879
00880 if( sx*sy != 1 )
00881 return false;
00882 x = abs(x);
00883 y = abs(y);
00884 return ( 1.f - min(x,y)/max(x,y) < ((sys_float)n+0.1f)*FLT_EPSILON );
00885 }
00886
00887 inline bool fp_equal( double x, double y, int n=3 )
00888 {
00889 ASSERT( n >= 1 );
00890
00891 if( isnan(x) || isnan(y) )
00892 return false;
00893 int sx = sign3(x);
00894 int sy = sign3(y);
00895
00896 if( sx == 0 && sy == 0 )
00897 return true;
00898
00899 if( sx*sy != 1 )
00900 return false;
00901 x = abs(x);
00902 y = abs(y);
00903 return ( 1. - min(x,y)/max(x,y) < ((double)n+0.1)*DBL_EPSILON );
00904 }
00905
00906 inline bool fp_equal_tol( sys_float x, sys_float y, sys_float tol )
00907 {
00908 ASSERT( tol > 0.f );
00909
00910 if( isnan(x) || isnan(y) )
00911 return false;
00912
00913 ASSERT( tol >= FLT_EPSILON*max(abs(x),abs(y)) );
00914 return ( abs( x-y ) <= tol );
00915 }
00916
00917 inline bool fp_equal_tol( double x, double y, double tol )
00918 {
00919 ASSERT( tol > 0. );
00920
00921 if( isnan(x) || isnan(y) )
00922 return false;
00923
00924 ASSERT( tol >= DBL_EPSILON*max(abs(x),abs(y)) );
00925 return ( abs( x-y ) <= tol );
00926 }
00927
00929 inline bool fp_bound( sys_float lo, sys_float x, sys_float hi, int n=3 )
00930 {
00931 ASSERT( n >= 1 );
00932
00933 if( isnan(x) || isnan(lo) || isnan(hi) )
00934 return false;
00935 if( fp_equal(lo,hi,n) )
00936 return fp_equal(0.5f*(lo+hi),x,n);
00937 if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -((sys_float)n+0.1f)*FLT_EPSILON )
00938 return false;
00939 return true;
00940 }
00941 inline bool fp_bound( double lo, double x, double hi, int n=3 )
00942 {
00943 ASSERT( n >= 1 );
00944
00945 if( isnan(x) || isnan(lo) || isnan(hi) )
00946 return false;
00947 if( fp_equal(lo,hi,n) )
00948 return fp_equal(0.5*(lo+hi),x,n);
00949 if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -((double)n+0.1)*DBL_EPSILON )
00950 return false;
00951 return true;
00952 }
00953 inline bool fp_bound_tol( sys_float lo, sys_float x, sys_float hi, sys_float tol )
00954 {
00955 ASSERT( tol > 0.f );
00956
00957 if( isnan(x) || isnan(lo) || isnan(hi) )
00958 return false;
00959 if( fp_equal_tol(lo,hi,tol) )
00960 return fp_equal_tol(0.5f*(lo+hi),x,tol);
00961 if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -tol )
00962 return false;
00963 return true;
00964 }
00965 inline bool fp_bound_tol( double lo, double x, double hi, double tol )
00966 {
00967 ASSERT( tol > 0. );
00968
00969 if( isnan(x) || isnan(lo) || isnan(hi) )
00970 return false;
00971 if( fp_equal_tol(lo,hi,tol) )
00972 return fp_equal_tol(0.5*(lo+hi),x,tol);
00973 if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -tol )
00974 return false;
00975 return true;
00976 }
00977
00978
00979 #undef POW2
00980
00981 #define POW2 pow2
00982 template<class T>
00983 inline T pow2(T a) { return a*a; }
00984
00985
00986 #undef POW3
00987
00988 #define POW3 pow3
00989 template<class T>
00990 inline T pow3(T a) { return a*a*a; }
00991
00992
00993 #undef POW4
00994
00995 #define POW4 pow4
00996 template<class T>
00997 inline T pow4(T a) { T b = a*a; return b*b; }
00998
00999
01000 #undef SDIV
01001
01004 inline sys_float SDIV( sys_float x ) { return ( fabs((double)x) < (double)SMALLFLOAT ) ? (sys_float)SMALLFLOAT : x; }
01005
01006 inline double SDIV( double x ) { return ( fabs(x) < (double)SMALLFLOAT ) ? (double)SMALLFLOAT : x; }
01007
01008
01009
01013 inline sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
01014 {
01015
01016 if( isnan(x) || isnan(y) )
01017 return x/y;
01018 int sx = sign3(x);
01019 int sy = sign3(y);
01020
01021 if( sx == 0 && sy == 0 )
01022 {
01023 if( isnan(res_0by0) )
01024 return x/y;
01025 else
01026 return res_0by0;
01027 }
01028 if( sx == 0 )
01029 return 0.;
01030 if( sy == 0 )
01031 return ( sx < 0 ) ? -FLT_MAX : FLT_MAX;
01032
01033 sys_float ay = abs(y);
01034 if( ay >= 1.f )
01035 return x/y;
01036 else
01037 {
01038
01039 if( abs(x) < ay*FLT_MAX )
01040 return x/y;
01041 else
01042 return ( sx*sy < 0 ) ? -FLT_MAX : FLT_MAX;
01043 }
01044 }
01045
01046 inline sys_float safe_div(sys_float x, sys_float y)
01047 {
01048 return safe_div( x, y, numeric_limits<sys_float>::quiet_NaN() );
01049 }
01050
01054 inline double safe_div(double x, double y, double res_0by0)
01055 {
01056
01057 if( isnan(x) || isnan(y) )
01058 return x/y;
01059 int sx = sign3(x);
01060 int sy = sign3(y);
01061
01062 if( sx == 0 && sy == 0 )
01063 {
01064 if( isnan(res_0by0) )
01065 return x/y;
01066 else
01067 return res_0by0;
01068 }
01069 if( sx == 0 )
01070 return 0.;
01071 if( sy == 0 )
01072 return ( sx < 0 ) ? -DBL_MAX : DBL_MAX;
01073
01074 double ay = abs(y);
01075 if( ay >= 1. )
01076 return x/y;
01077 else
01078 {
01079
01080 if( abs(x) < ay*DBL_MAX )
01081 return x/y;
01082 else
01083 return ( sx*sy < 0 ) ? -DBL_MAX : DBL_MAX;
01084 }
01085 }
01086
01087 inline double safe_div(double x, double y)
01088 {
01089 return safe_div( x, y, numeric_limits<double>::quiet_NaN() );
01090 }
01091
01092 #undef HMRATE
01093
01094
01095
01096
01097
01098 #define HMRATE(a,b,c) hmrate4(a,b,c,phycon.te)
01099
01100 inline double hmrate4( double a, double b, double c, double te )
01101 {
01102 if( b == 0. && c == 0. )
01103 return a;
01104 else if( c == 0. )
01105 return a*pow(te/300.,b);
01106 else if( b == 0. )
01107 return ( c/te <= 50. ) ? a*exp(-c/te) : 0.;
01108 else
01109 return ( c/te <= 50. ) ? a*pow(te/300.,b)*exp(-c/te) : 0.;
01110 }
01111
01112 template<class T>
01113 inline void invalidate_array(T* p, size_t size)
01114 {
01115 if( size > 0 )
01116 memset( p, -1, size );
01117 }
01118
01119 inline void invalidate_array(double* p, size_t size)
01120 {
01121 set_NaN( p, (long)(size/sizeof(double)) );
01122 }
01123
01124 inline void invalidate_array(sys_float* p, size_t size)
01125 {
01126 set_NaN( p, (long)(size/sizeof(sys_float)) );
01127 }
01128
01131 template<class T> inline T* get_ptr(valarray<T> &v)
01132 {
01133 return &v[0];
01134 }
01135 template<class T> inline T* get_ptr(vector<T> &v)
01136 {
01137 return &v[0];
01138 }
01139 template<class T> inline const T* get_ptr(const valarray<T> &v)
01140 {
01141 return const_cast<const T*>(&const_cast<valarray<T>&>(v)[0]);
01142 }
01143 template<class T> inline const T* get_ptr(const vector<T> &v)
01144 {
01145 return const_cast<const T*>(&const_cast<vector<T>&>(v)[0]);
01146 }
01147
01173 template<class T>
01174 class auto_vec
01175 {
01176 T* ptr;
01177
01178 template<class U>
01179 struct auto_vec_ref
01180 {
01181 U* ptr;
01182
01183 explicit auto_vec_ref( U* p )
01184 {
01185 ptr = p;
01186 }
01187 };
01188
01189 public:
01190 typedef T element_type;
01191
01192
01193
01194 explicit auto_vec( element_type* p = NULL ) throw()
01195 {
01196 ptr = p;
01197 }
01198 auto_vec( auto_vec& p ) throw()
01199 {
01200 ptr = p.release();
01201 }
01202 auto_vec& operator= ( auto_vec& p ) throw()
01203 {
01204 reset( p.release() );
01205 return *this;
01206 }
01207 ~auto_vec() throw()
01208 {
01209 delete[] ptr;
01210 }
01211
01212
01213
01214 element_type& operator[] ( ptrdiff_t n ) const throw()
01215 {
01216 return *(ptr+n);
01217 }
01218 element_type* get() const throw()
01219 {
01220 return ptr;
01221 }
01222
01223 element_type* data() const throw()
01224 {
01225 return ptr;
01226 }
01227 element_type* release() throw()
01228 {
01229 element_type* p = ptr;
01230 ptr = NULL;
01231 return p;
01232 }
01233 void reset( element_type* p = NULL ) throw()
01234 {
01235 if( p != ptr )
01236 {
01237 delete[] ptr;
01238 ptr = p;
01239 }
01240 }
01241
01242
01243
01244 auto_vec( auto_vec_ref<element_type> r ) throw()
01245 {
01246 ptr = r.ptr;
01247 }
01248 auto_vec& operator= ( auto_vec_ref<element_type> r ) throw()
01249 {
01250 if( r.ptr != ptr )
01251 {
01252 delete[] ptr;
01253 ptr = r.ptr;
01254 }
01255 return *this;
01256 }
01257 operator auto_vec_ref<element_type>() throw()
01258 {
01259 return auto_vec_ref<element_type>( this->release() );
01260 }
01261 };
01262
01263 #include "container_classes.h"
01264 #include "iter_track.h"
01265
01266
01267
01268
01269
01270
01271
01272 class transition;
01273 typedef struct t_quantumState quantumState;
01274 typedef struct t_emission emission;
01275 typedef struct t_collision collision;
01276 typedef struct t_species species;
01277
01278 struct t_emission
01279 {
01288 int iRedisFun;
01289
01291 long int ipFine;
01292
01306 realnum TauIn;
01307
01318 realnum TauTot;
01319
01324 iter_track_basic<realnum> TauTrack;
01325
01331 realnum TauCon;
01332
01334 realnum FracInwd;
01335
01338 double pump;
01339
01341 double xIntensity;
01342
01344 double phots;
01345
01347 realnum gf;
01348
01350 realnum Pesc;
01351
01353 realnum Pelec_esc;
01354
01356 realnum Pdest;
01357
01361 realnum dampXvel;
01362
01364 realnum damp;
01365
01367 double ColOvTot;
01368
01370 realnum AutoIonizFrac;
01371
01377 realnum opacity;
01378
01380 double PopOpc;
01381
01383 realnum Aul;
01384
01386 double ots;
01387
01388
01389 transition *tran;
01390
01391
01392 emission *next;
01393 };
01394
01395
01396
01397
01398
01399 struct t_species
01400 {
01401
01402 char *chLabel;
01403
01404 long numLevels_max;
01405
01406 long numLevels_local;
01407
01408 realnum fmolweight;
01409
01410 bool lgMolecular;
01411
01412 realnum fracType;
01413
01414 realnum fracIsotopologue;
01416 double CoolTotal;
01418 bool lgActive;
01419 };
01420
01421
01422
01423 typedef struct t_CollRatesArray
01424 {
01425
01426 long ntemps;
01427
01428 double *temps;
01429
01430 double ***collrates;
01431
01432 } CollRateCoeffArray ;
01433
01434
01435
01436 typedef struct t_CollSplinesArray
01437 {
01438
01439
01440
01441
01442
01443 double *collspline;
01444 double *SplineSecDer;
01445
01446 long nSplinePts;
01447 long intTranType;
01448 double gf;
01449 double EnergyDiff;
01450 double ScalingParam;
01451
01452 } CollSplinesArray ;
01453
01454 struct t_collision
01455 {
01456
01457
01458
01459
01461 realnum ColUL;
01462
01464 realnum col_str;
01465
01467 realnum col_stri[ipNCOLLIDER];
01468
01470 double cool , heat;
01471 };
01472
01473 #include "energy.h"
01474
01475 struct t_quantumState
01476 {
01477 char chLabel[5];
01478 char chConfig[11];
01479
01480 species *sp;
01481
01483 Energy energy;
01484
01486 realnum g;
01487
01489 double Pop;
01490
01492 int IonStg;
01494 int nelem;
01495
01497 double ConBoltz;
01498
01499
01500 long n, l, S, j;
01501
01503 double lifetime;
01504
01505
01506 quantumState *next;
01507 };
01508
01509
01510
01511 INSTANTIATE_MULTI_ARR( quantumState, lgBOUNDSCHECKVAL );
01512
01513
01514
01515
01516
01517
01518
01519
01526 typedef enum { SPM_RELAX, SPM_KEEP_EMPTY, SPM_STRICT } split_mode;
01527
01529 void Split(const string& str,
01530 const string& sep,
01531 vector<string>& lst,
01532 split_mode mode);
01533
01536 inline bool FindAndReplace(string& str,
01537 const string& substr,
01538 const string& newstr)
01539 {
01540 string::size_type ptr = str.find( substr );
01541 if( ptr != string::npos )
01542 str.replace( ptr, substr.length(), newstr );
01543 return ptr != string::npos;
01544 }
01545
01548 inline bool FindAndErase(string& str,
01549 const string& substr)
01550 {
01551 return FindAndReplace( str, substr, "" );
01552 }
01553
01559 double csphot(long int inu, long int ithr, long int iofset);
01560
01565 double RandGauss(double xMean, double s );
01566
01570 double MyGaussRand( double PctUncertainty );
01571
01573 double AnuUnit(realnum energy);
01574
01579 void cap4(char *chCAP , const char *chLab);
01580
01583 void uncaps(char *chCard );
01584
01587 void caps(char *chCard );
01588
01591 double e2(
01592 double t );
01593
01596 double ee1(double x);
01597
01601 double ee1_safe(double x);
01602
01609 double FFmtRead(const char *chCard,
01610 long int *ipnt,
01611 long int last,
01612 bool *lgEOL);
01613
01619 long nMatch(const char *chKey,
01620 const char *chCard);
01621
01631 int GetQuote( char *chLabel, char *chCard, char *chCardRaw, bool lgABORT );
01632
01633
01634 inline const char *strstr_s(const char *haystack, const char *needle)
01635 {
01636 return const_cast<const char *>(strstr(haystack, needle));
01637 }
01638
01639 inline char *strstr_s(char *haystack, const char *needle)
01640 {
01641 return const_cast<char *>(strstr(haystack, needle));
01642 }
01643
01644 inline const char *strchr_s(const char *s, int c)
01645 {
01646 return const_cast<const char *>(strchr(s, c));
01647 }
01648
01649 inline char *strchr_s(char *s, int c)
01650 {
01651 return const_cast<char *>(strchr(s, c));
01652 }
01653
01656 long int ipow( long, long );
01657
01660 void PrintE82( FILE*, double );
01661
01663 void PrintE71( FILE*, double );
01664
01666 void PrintE93( FILE*, double );
01667
01673
01674 #ifdef _MSC_VER
01675 char *PrintEfmt(const char *fmt, double val );
01676 #else
01677 #define PrintEfmt( F, V ) F, V
01678 #endif
01679
01681 const double SEXP_LIMIT = 84.;
01683 const double DSEXP_LIMIT = 680.;
01684
01686 sys_float sexp(sys_float x);
01687 double sexp(double x);
01688
01693 double dsexp(double x);
01694
01699 double plankf(long int ip);
01700
01701
01702 typedef enum { Gaussian32, Legendre } methods;
01703
01704
01705 template<typename Integrand, methods Method>
01706 class Integrator
01707 {
01708 public:
01709 double numPoints, weights[16], c[16];
01710
01711 Integrator( void )
01712 {
01713 numPoints = 16;
01714 double weights_temp[16] = {
01715 .35093050047350483e-2, .81371973654528350e-2, .12696032654631030e-1, .17136931456510717e-1,
01716 .21417949011113340e-1, .25499029631188088e-1, .29342046739267774e-1, .32911111388180923e-1,
01717 .36172897054424253e-1, .39096947893535153e-1, .41655962113473378e-1, .43826046502201906e-1,
01718 .45586939347881942e-1, .46922199540402283e-1, .47819360039637430e-1, .48270044257363900e-1};
01719
01720 double c_temp[16] = {
01721 .498631930924740780, .49280575577263417, .4823811277937532200, .46745303796886984000,
01722 .448160577883026060, .42468380686628499, .3972418979839712000, .36609105937014484000,
01723 .331522133465107600, .29385787862038116, .2534499544661147000, .21067563806531767000,
01724 .165934301141063820, .11964368112606854, .7223598079139825e-1, .24153832843869158e-1};
01725
01726 for( long i=0; i<numPoints; i++ )
01727 {
01728 weights[i] = weights_temp[i];
01729 c[i] = c_temp[i];
01730 }
01731 return;
01732 };
01733 double sum(double min, double max, Integrand func)
01734 {
01735 ASSERT( Method == Gaussian32 );
01736 double a = 0.5*(max+min),
01737 b = max-min,
01738 total = 0.;
01739
01740 for( long i=0; i< numPoints; i++ )
01741 total += b * weights[i] * ( func(a+b*c[i]) + func(a-b*c[i]) );
01742
01743 return total;
01744 }
01745 };
01746
01753 double qg32( double, double, double(*)(double) );
01754
01755
01756
01757
01767 void spsort( realnum x[], long int n, long int iperm[], int kflag, int *ier);
01768
01769
01770 inline double vfun(double damp, double x)
01771 {
01772
01773 return sexp(x*x) + damp/1.772453850905516027298167/(1. + x*x);
01774 }
01775
01776
01777
01778
01779
01780
01781
01782
01783
01784 #ifdef _MSC_VER
01785
01786 # pragma warning( disable : 4127 )
01787
01788 # pragma warning( disable : 4996 )
01789
01790 # pragma warning( disable : 4056 )
01791
01792 # pragma warning( disable : 4514 )
01793
01794
01795 # pragma warning( disable : 4512 )
01796 #endif
01797 #ifdef __INTEL_COMPILER
01798 # pragma warning( disable : 1572 )
01799 #endif
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811 #endif
01812