Cloudy
Spectral Synthesis Code for Astrophysics
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cddefines.h
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2023 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 
4 #ifndef CDDEFINES_H_
5 #define CDDEFINES_H_
6 
7 #include "cdstd.h"
8 
9 #ifdef _MSC_VER
10 /* we are not using MS foundation class */
11 # ifndef WIN32_LEAN_AND_MEAN
12 # define WIN32_LEAN_AND_MEAN
13 # endif
14 #endif
15 
16 #ifdef __clang__
17 // this would generate lots of warnings about mismatched tags in the STL valarray definition
18 #pragma clang diagnostic ignored "-Wmismatched-tags"
19 #endif
20 
21 /* these headers are needed by all files */
22 /*lint -e129 these resolve several issues pclint has with my system headers */
23 /*lint -e78 */
24 /*lint -e830 */
25 /*lint -e38 */
26 /*lint -e148 */
27 /*lint -e114 */
28 /*lint -e18 */
29 /*lint -e49 */
30 // C++ versions of C headers
31 #include <cstdio>
32 #include <cstdlib>
33 #include <cctype>
34 #ifdef _MSC_VER
35 // MSVC needs this before cmath in order provide numeric constants
36 // (M_PI etc.) defined by C99 but not C++ standards to date.
37 #define _USE_MATH_DEFINES
38 #endif
39 #include <cmath>
40 #include <cassert>
41 #include <cstring>
42 #include <cfloat>
43 #include <climits>
44 #include <ctime>
45 #if defined(__sun) && defined(__SUNPRO_CC)
46 // with Solaris Studio 12.2 under Sparc Solaris, csignal doesn't define sigaction...
47 #include <signal.h>
48 #else
49 #include <csignal>
50 #endif
51 // C++ headers
52 #include <limits>
53 #include <string>
54 #include <sstream>
55 #include <iomanip>
56 #include <array>
57 #include <vector>
58 #include <valarray>
59 #include <complex>
60 #include <map>
61 #include <tuple>
62 #include <regex>
63 #include <memory>
64 #include <stdexcept>
65 #include <algorithm>
66 #include <fstream>
67 #include <bitset>
68 #include <unordered_map>
69 #include <numeric>
70 
71 // Workaround for Windows...
72 #if defined(_MSC_VER) && !defined(SYS_CONFIG)
73 #define SYS_CONFIG "cloudyconfig_vs.h"
74 #endif
75 
76 // platform specific configuration; generated by configure.sh
77 #ifdef SYS_CONFIG
78 #include SYS_CONFIG
79 #else
80 #include "cloudyconfig.h"
81 #endif
82 
83 // prevent problems on platforms with broken support (such as Mac homebrew g++)
84 #ifndef HAVE_AVX_INTRIN
85 #undef __AVX__
86 #endif
87 
88 #ifndef HAVE_FMA_INTRIN
89 #undef __FMA__
90 #endif
91 
92 #ifndef HAVE_AVX2_INTRIN
93 #undef __AVX2__
94 #endif
95 
96 #ifndef HAVE_AVX512F_INTRIN
97 #undef __AVX512F__
98 #endif
99 
100 #ifdef __AVX__
101 #include <immintrin.h>
102 #endif
103 
104 /*lint +e18 */
105 /*lint +e49 */
106 /*lint +e38 */
107 /*lint +e148 */
108 /*lint +e830 */
109 /*lint +e78 */
110 /*lint -e129 */
111 
112 using namespace std;
113 
114 #undef NULL
115 #define NULL nullptr
116 
117 #undef STATIC
118 #ifdef USE_GPROF
119 #define STATIC
120 #else
121 #define STATIC static
122 #endif
123 
124 #ifdef FLT_IS_DBL
125 typedef double realnum;
126 #else
127 typedef float realnum;
128 #endif
129 
130 typedef float sys_float;
131 // prevent explicit float's from creeping back into the code
132 #define float PLEASE_USE_REALNUM_NOT_FLOAT
133 
134 // define realnum literals, use as 12_r or 12.4_r
135 inline realnum operator "" _r( unsigned long long l )
136 {
137  return realnum(l);
138 }
139 
140 inline realnum operator "" _r( long double l )
141 {
142  return realnum(l);
143 }
144 
145 inline FILE *sys_fopen(const char *path, const char *mode)
146 {
147  return fopen( path, mode );
148 }
149 #define fopen PLEASE_USE_open_data_NOT_fopen
150 
151 typedef enum {
152  ES_SUCCESS=0, // everything went fine...
153  ES_FAILURE=1, // general failure exit
154  ES_WARNINGS, // warnings were present
155  ES_BOTCHES, // botched monitors were present
156  ES_CLOUDY_ABORT, // Cloudy aborted
157  ES_BAD_ASSERT, // an assert in the code failed
158  ES_BAD_ALLOC, // a memory allocation failed
159  ES_OUT_OF_RANGE, // an out-of-range exception was thrown
160  ES_DOMAIN_ERROR, // a vectorized math routine threw a domain error
161  ES_ILLEGAL_INSTRUCTION, // the CPU encountered an illegal instruction
162  ES_FP_EXCEPTION, // a floating point exception was caught
163  ES_SEGFAULT, // a segmentation fault occurred
164  ES_BUS_ERROR, // a bus error occurred
165  ES_UNKNOWN_SIGNAL, // an unknown signal was caught
166  ES_UNKNOWN_EXCEPTION, // an unknown exception was caught
167  ES_TOP // NB NB -- this should always be the last entry
168 } exit_type;
169 
170 // make sure the system definitions are on par with ours
171 // especially EXIT_FAILURE does not have a guaranteed value!
172 #undef EXIT_SUCCESS
173 #define EXIT_SUCCESS ES_SUCCESS
174 #undef EXIT_FAILURE
175 #define EXIT_FAILURE ES_FAILURE
176 
177 #ifdef BOUNDS_CHECK
178 #define lgBOUNDSCHECKVAL true
179 #else
180 #define lgBOUNDSCHECKVAL false
181 #endif
182 
183 /* make sure this is globally visible as well! */
184 /* This must be done at the start of every file, to ensure that policy
185  for FPE handling, etc., is guaranteed to be set up before the
186  construction of file-statics and globals. */
187 #include "cpu.h"
188 
189 //*************************************************************************
190 //
208 //
209 // This implementation has been obtained from Wikipedia
210 //
211 //*************************************************************************
212 
213 template<typename T> class Singleton
214 {
215 public:
216  static T& Inst()
217  {
218  static T instance; // assumes T has a protected default constructor
219  return instance;
220  }
221 };
222 
223 /**************************************************************************
224  *
225  * these are variables and pointers for output from the code, used everywhere
226  * declared extern here, and definition is in cddefines.cpp
227  *
228  **************************************************************************/
229 
234 class Output
235 {
236 public:
237  FILE *m_fp;
238  // Implicit conversion
239  Output(FILE* fp) : m_fp(fp) {}
240  FILE *fptr() const
241  {
242  return m_fp;
243  }
244 };
245 
246 extern FILE *ioQQQ;
247 
248 extern FILE *ioStdin;
249 
250 extern FILE *ioMAP;
251 
254 extern FILE* ioPrnErr;
255 
258 extern bool lgTestCodeCalled;
259 
262 extern bool lgTestCodeEnabled;
263 
266 extern bool lgPrnErr;
267 
270 extern long int nzone;
271 
273 extern double fnzone;
274 
277 extern long int iteration;
278 
285 extern const double ZeroNum;
286 
291 extern const void* ZeroPtr;
292 
293 /**************************************************************************
294  *
295  * these are constants used to dimension several vectors and index arrays
296  *
297  **************************************************************************/
298 
303 const int FILENAME_PATH_LENGTH = 200;
304 
307 
311 const int INPUT_LINE_LENGTH = 2000;
312 
314 const int NCHLAB = 20;
315 
318 const int LIMELM = 30;
319 
321 const int NISO = 2;
322 
326 const int NHYDRO_MAX_LEVEL = 401;
327 
329 const double MAX_DENSITY = 1.e24;
330 
332 const double DEPTH_OFFSET = 1.e-30;
333 
334 enum {CHARS_SPECIES=10};
335 enum {CHARS_ISOTOPE_SYM = 6};
336 
337 /* indices within recombination coefficient array */
338 /* ipRecEsc is state specific escape probability*/
339 const int ipRecEsc = 2;
340 /* the net escaping, including destruction by background and optical deepth*/
341 const int ipRecNetEsc = 1;
342 /* ipRecRad is state specific radiative recombination rate*/
343 const int ipRecRad = 0;
348 /* these specify the form of the line redistribution function */
349 /* partial redistribution with wings */
350 const int ipPRD = 1;
351 /* complete redistribution, core only, no wings, Hummer's K2 function */
352 const int ipCRD = -1;
353 /* complete redistribution with wings */
354 const int ipCRDW = 2;
355 /* redistribution function for Lya, calls Hummer routine for H-like series only */
356 const int ipLY_A = -2;
357 
359 const int ipHYDROGEN = 0;
360 const int ipHELIUM = 1;
361 const int ipLITHIUM = 2;
362 const int ipBERYLLIUM = 3;
363 const int ipBORON = 4;
364 const int ipCARBON = 5;
365 const int ipNITROGEN = 6;
366 const int ipOXYGEN = 7;
367 const int ipFLUORINE = 8;
368 const int ipNEON = 9;
369 const int ipSODIUM = 10;
370 const int ipMAGNESIUM = 11;
371 const int ipALUMINIUM = 12;
372 const int ipSILICON = 13;
373 const int ipPHOSPHORUS = 14;
374 const int ipSULPHUR = 15;
375 const int ipCHLORINE = 16;
376 const int ipARGON = 17;
377 const int ipPOTASSIUM = 18;
378 const int ipCALCIUM = 19;
379 const int ipSCANDIUM = 20;
380 const int ipTITANIUM = 21;
381 const int ipVANADIUM = 22;
382 const int ipCHROMIUM = 23;
383 const int ipMANGANESE = 24;
384 const int ipIRON = 25;
385 const int ipCOBALT = 26;
386 const int ipNICKEL = 27;
387 const int ipCOPPER = 28;
388 const int ipZINC = 29;
389 const int ipKRYPTON = 35;
390 
391 /***************************************************************************
392  * the following are prototypes for some routines that are part of the
393  * debugging process - they come and go in any particular sub.
394  * it is not necessary to declare them when used since they are defined here
395  **************************************************************************/
396 
405 double fudge(long int ipnt);
406 
410 void broken(void);
411 
414 void fixit_base(const char* func, const char* file, int line,
415  const char *reason);
416 
417 class Fixit
418 {
419 public:
420  Fixit(const char* func, const char* file, int line,
421  const char *reason)
422  {
423  fixit_base(func,file,line,reason);
424  }
425 };
426 
427 #define fixit(a) \
428  do { \
429  static Fixit fixit_s(__func__,__FILE__,__LINE__, (a)); \
430  } while (0)
431 
433 void CodeReview(void);
434 
436 void TestCode(void);
437 
442 void MyAssert(const char *file, int line, const char *comment);
443 
446 
448 {
449  const char* p_routine;
450  const char* p_file;
451  long p_line;
453 public:
454  cloudy_exit(const char* routine, const char* file, long line, exit_type exit_code)
455  {
456  p_routine = routine;
457  p_file = file;
458  p_line = line;
459  p_exit = exit_code;
460  }
461  const char* routine() const throw()
462  {
463  return p_routine;
464  }
465  const char* file() const throw()
466  {
467  return p_file;
468  }
469  long line() const
470  {
471  return p_line;
472  }
474  {
475  return p_exit;
476  }
477 };
478 
479 // workarounds for __func__ are defined in cpu.h
480 #define cdEXIT( FAIL ) throw cloudy_exit( __func__, __FILE__, __LINE__, FAIL )
481 
482 // calls like puts( "[Stop in MyRoutine]" ) have been integrated in cdEXIT above
483 #define puts( STR ) Using_puts_before_cdEXIT_is_no_longer_needed
484 
486 void ShowMe(void);
487 
489 NORETURN void TotalInsanity(void);
490 
491 /* TotalInsanityAsStub always calls TotalInsanity(), but in such a way that
492  * it can be used as a stub for another routine without generating warnings
493  * about unreachable code after the stub. Hence this should NOT be NORETURN */
494 template<class T>
496 {
497  // this is always true...
498  if( ZeroNum == 0. )
499  TotalInsanity();
500  else
501  return T();
502 }
503 
505 NORETURN void BadRead(void);
506 
510 int dbg_printf(int debug, const char *fmt, ...);
511 
513 int dprintf(FILE *fp, const char *format, ...);
514 
516 inline void cdBacktrace()
517 {
519  cpu.i().PrintBacktrace("DEBUG ", false);
520 }
521 
523 int fprintf (const Output& stream, const char *format, ...);
524 
525 int dprintf(const Output & stream, const char *format, ...);
526 
532 bool read_whole_line( string& chLine, FILE *ioIN );
533 
534 #ifdef HAVE_LIBCPP_BUG
535 
536 #include "service.h"
537 
538 // workaround for the bug described in https://bugs.llvm.org/show_bug.cgi?id=17782
539 inline istringstream& operator>> ( istringstream& s, double& x )
540 {
541  FPRead(s, s.str(), x);
542  return s;
543 }
544 
545 inline istringstream& operator>> ( istringstream& s, sys_float& x )
546 {
547  double y;
548  FPRead(s, s.str(), y);
549  x = sys_float(y);
550  return s;
551 }
552 
553 #endif
554 
555 /**************************************************************************
556  *
557  * various macros used by the code
558  *
559  **************************************************************************/
560 
564 #ifndef NDEBUG
565 # define DEBUG
566 #else
567 # undef DEBUG
568 #endif
569 
571 {
572  int p_sig;
573 public:
574  bad_signal(int sig, void* ptr);
575  bad_signal(const bad_signal&) = default;
576  bad_signal& operator= (const bad_signal&) = default;
577  virtual ~bad_signal() throw() {}
578  int sig() const throw()
579  {
580  return p_sig;
581  }
582 };
583 
585 {
586  const char* p_file;
587  long p_line;
588  const char* p_comment;
589 public:
590  bad_assert(const char* file, long line, const char* comment);
591  bad_assert(const bad_assert&) = default;
592  bad_assert& operator= (const bad_assert&) = default;
593  virtual ~bad_assert() throw() {}
594  void print(void) const
595  {
596  fprintf(ioQQQ,"DISASTER Assertion failure at %s:%ld\n%s\n",
597  p_file, p_line, p_comment);
598  }
599  const char* file() const throw()
600  {
601  return p_file;
602  }
603  long line() const throw()
604  {
605  return p_line;
606  }
607  const char *comment() const throw()
608  {
609  return p_comment;
610  }
611 };
612 
614 {
615  const char* p_comment;
616 public:
617  explicit cloudy_abort(const char* comment);
618  cloudy_abort(const cloudy_abort&) = default;
619  cloudy_abort& operator= (const cloudy_abort&) = default;
620  virtual ~cloudy_abort() throw() {}
621  const char *comment() const throw()
622  {
623  return p_comment;
624  }
625 };
626 
627 /* the do { ... } while ( 0 ) construct prevents bugs in code like this:
628  * if( test )
629  * ASSERT( n == 10 );
630  * else
631  * do something else...
632  */
633 #undef ASSERT
634 #if NDEBUG
635 # define ASSERT(exp) ((void)0)
636 #else
637 # define ASSERT(exp) \
638  do { \
639  if (UNLIKELY(!(exp))) \
640  throw bad_assert(__FILE__,__LINE__,"Failed: " #exp); \
641  } while( 0 )
642 #endif
643 
644 #define MESSAGE_ASSERT(msg, exp) ASSERT( (msg) ? (exp) : false )
645 
646 inline NORETURN void OUT_OF_RANGE(const char* str)
647 {
649  throw out_of_range( str );
650 }
651 
652 inline NORETURN void DOMAIN_ERROR(const string& str)
653 {
655  throw domain_error( str );
656 }
657 
658 #ifdef __SUNPRO_CC
659 #pragma does_not_return(OUT_OF_RANGE)
660 #endif
661 
662 /* Windows does not define isnan */
663 /* use our version on all platforms since the isnanf
664  * function does not exist under Solaris 9 either */
665 #undef isnan
666 #define isnan MyIsnan
667 
670 class t_debug : public Singleton<t_debug>
671 {
672  friend class Singleton<t_debug>;
673  FILE *p_fp;
675 protected:
676  t_debug() : p_fp(stderr)
677  {
678  p_callLevel = 0;
679  }
680 public:
681  void enter(const char *name)
682  {
683  ++p_callLevel;
684  fprintf(p_fp,"%*c%s\n",p_callLevel,'>',name);
685  }
686  void leave(const char *name)
687  {
688  fprintf(p_fp,"%*c%s\n",p_callLevel,'<',name);
689  --p_callLevel;
690  }
691 };
692 
695 class t_nodebug : public Singleton<t_nodebug>
696 {
697  friend class Singleton<t_nodebug>;
698 protected:
700 public:
701  void enter(const char *) const {}
702  void leave(const char *) const {}
703 };
704 
705 template<class Trace>
707 {
708  const char *p_name;
709 public:
710  explicit debugtrace(const char *funcname)
711  {
712  p_name = funcname;
713  Trace::Inst().enter(p_name);
714  }
715  debugtrace(const debugtrace&) = default;
716  debugtrace& operator= (const debugtrace&) = default;
718  {
719  Trace::Inst().leave(p_name);
720  }
721  const char* name() const
722  {
723  return p_name;
724  }
725 };
726 
727 #ifdef DEBUG_FUN
728 #define DEBUG_ENTRY( funcname ) debugtrace<t_debug> DEBUG_ENTRY( funcname )
729 #else
730 #define DEBUG_ENTRY( funcname ) ((void)0)
731 #endif
732 
733 // overload the character manipulation routines
734 inline char tolower(char c)
735 {
736  return static_cast<char>( tolower( static_cast<int>(c) ) );
737 }
738 inline unsigned char tolower(unsigned char c)
739 {
740  return static_cast<unsigned char>( tolower( static_cast<int>(c) ) );
741 }
742 
743 inline char toupper(char c)
744 {
745  return static_cast<char>( toupper( static_cast<int>(c) ) );
746 }
747 inline unsigned char toupper(unsigned char c)
748 {
749  return static_cast<unsigned char>( toupper( static_cast<int>(c) ) );
750 }
751 
752 /* TorF(l) returns a 'T' or 'F' depending on the 'logical' expr 'l' */
753 inline char TorF( bool l ) { return l ? 'T' : 'F'; }
754 /* */
755 
757 inline bool is_odd( int j ) { return (j&1) == 1; }
758 inline bool is_odd( long j ) { return (j&1L) == 1L; }
759 /* */
760 
762 inline long nint( double x ) { return static_cast<long>( (x < 0.) ? x-0.5 : x+0.5 ); }
763 /* */
764 
765 /* define min for mixed arguments, the rest already exists */
766 inline long min( int a, long b ) { long c = a; return ( (c < b) ? c : b ); }
767 inline long min( long a, int b ) { long c = b; return ( (a < c) ? a : c ); }
768 inline double min( sys_float a, double b ) { double c = a; return ( (c < b) ? c : b ); }
769 inline double min( double a, sys_float b ) { double c = b; return ( (a < c) ? a : c ); }
770 
771 /* want to define this only if no native os support exists */
772 #ifndef HAVE_POWI
773 
774 double powi( double , long int );
775 #endif
776 
778 double powpq(double x, int p, int q);
779 
780 /* avoid ambiguous overloads */
781 #ifndef HAVE_POW_DOUBLE_INT
782 inline double pow( double x, int i ) { return powi( x, long(i) ); }
783 #endif
784 
785 #ifndef HAVE_POW_DOUBLE_LONG
786 inline double pow( double x, long i ) { return powi( x, i ); }
787 #endif
788 
789 #ifndef HAVE_POW_FLOAT_INT
790 inline sys_float pow( sys_float x, int i ) { return sys_float( powi( double(x), long(i) ) ); }
791 #endif
792 
793 #ifndef HAVE_POW_FLOAT_LONG
794 inline sys_float pow( sys_float x, long i ) { return sys_float( powi( double(x), i ) ); }
795 #endif
796 
797 #ifndef HAVE_POW_FLOAT_DOUBLE
798 inline double pow( sys_float x, double y ) { return pow( double(x), y ); }
799 #endif
800 
801 #ifndef HAVE_POW_DOUBLE_FLOAT
802 inline double pow( double x, sys_float y ) { return pow( x, double(y) ); }
803 #endif
804 
805 #undef MIN2
806 
807 #define MIN2(a,b) min(a,b)
808 /* */
809 
810 #undef MIN3
811 
812 #define MIN3(a,b,c) (min(min(a,b),c))
813 /* */
814 
815 #undef MIN4
816 
817 #define MIN4(a,b,c,d) (min(min(a,b),min(c,d)))
818 /* */
819 
820 /* define max for mixed arguments, the rest already exists */
821 inline long max( int a, long b ) { long c = a; return ( (c > b) ? c : b ); }
822 inline long max( long a, int b ) { long c = b; return ( (a > c) ? a : c ); }
823 inline double max( sys_float a, double b ) { double c = a; return ( (c > b) ? c : b ); }
824 inline double max( double a, sys_float b ) { double c = b; return ( (a > c) ? a : c ); }
825 
826 #undef MAX2
827 
828 #define MAX2(a,b) max(a,b)
829 /* */
830 
831 #undef MAX3
832 
833 #define MAX3(a,b,c) (max(max(a,b),c))
834 /* */
835 
836 #undef MAX4
837 
838 #define MAX4(a,b,c,d) (max(max(a,b),max(c,d)))
839 /* */
840 
841 template<class T>
842 inline void vzero(vector<T>& vec)
843 {
844  memset(vec.data(), 0, vec.size()*sizeof(T));
845 }
846 
851 template<class T>
852 inline T sign( T x, T y )
853 {
854  return ( y < T() ) ? -abs(x) : abs(x);
855 }
856 /* */
857 
859 template<class T>
860 inline int sign3( T x ) { return ( x < T() ) ? -1 : ( ( x > T() ) ? 1 : 0 ); }
861 /* */
862 
864 inline bool fp_equal( sys_float x, sys_float y, int n=3 )
865 {
866  ASSERT( n >= 1 );
867  // mimic IEEE behavior
868  if( isnan(x) || isnan(y) )
869  return false;
870  int sx = sign3(x);
871  int sy = sign3(y);
872  // treat zero cases first to avoid division by zero below
873  if( sx == 0 && sy == 0 )
874  return true;
875  // either x or y is zero (but not both), or x and y have different sign
876  if( sx*sy != 1 )
877  return false;
878  x = abs(x);
879  y = abs(y);
880  // adjust epsilon value for denormalized numbers
881  int exp;
882  (void)frexp(max(x,y), &exp);
883  sys_float epsilon = ldexp(FLT_EPSILON, max(0,FLT_MIN_EXP-exp));
884  return ( 1.f - min(x,y)/max(x,y) < ((sys_float)n+0.1f)*epsilon );
885 }
886 
887 inline bool fp_equal( double x, double y, int n=3 )
888 {
889  ASSERT( n >= 1 );
890  // mimic IEEE behavior
891  if( isnan(x) || isnan(y) )
892  return false;
893  int sx = sign3(x);
894  int sy = sign3(y);
895  // treat zero cases first to avoid division by zero below
896  if( sx == 0 && sy == 0 )
897  return true;
898  // either x or y is zero (but not both), or x and y have different sign
899  if( sx*sy != 1 )
900  return false;
901  x = abs(x);
902  y = abs(y);
903  // adjust epsilon value for denormalized numbers
904  int exp;
905  (void)frexp(max(x,y), &exp);
906  double epsilon = ldexp(DBL_EPSILON, max(0,DBL_MIN_EXP-exp));
907  return ( 1. - min(x,y)/max(x,y) < ((double)n+0.1)*epsilon );
908 }
909 
910 inline bool fp_equal_tol( sys_float x, sys_float y, sys_float tol )
911 {
912  ASSERT( tol > 0.f );
913  // mimic IEEE behavior
914  if( isnan(x) || isnan(y) )
915  return false;
916  // make sure the tolerance is not too stringent
917  ASSERT( tol >= FLT_EPSILON*max(abs(x),abs(y)) );
918  return ( abs( x-y ) <= tol );
919 }
920 
921 inline bool fp_equal_tol( double x, double y, double tol )
922 {
923  ASSERT( tol > 0. );
924  // mimic IEEE behavior
925  if( isnan(x) || isnan(y) )
926  return false;
927  // make sure the tolerance is not too stringent
928  ASSERT( tol >= DBL_EPSILON*max(abs(x),abs(y)) );
929  return ( abs( x-y ) <= tol );
930 }
931 
933 inline bool fp_bound( sys_float lo, sys_float x, sys_float hi, int n=3 )
934 {
935  ASSERT( n >= 1 );
936  // mimic IEEE behavior
937  if( isnan(x) || isnan(lo) || isnan(hi) )
938  return false;
939  if( fp_equal(lo,hi,n) )
940  return fp_equal(0.5f*(lo+hi),x,n);
941  if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -((sys_float)n+0.1f)*FLT_EPSILON )
942  return false;
943  return true;
944 }
945 inline bool fp_bound( double lo, double x, double hi, int n=3 )
946 {
947  ASSERT( n >= 1 );
948  // mimic IEEE behavior
949  if( isnan(x) || isnan(lo) || isnan(hi) )
950  return false;
951  if( fp_equal(lo,hi,n) )
952  return fp_equal(0.5*(lo+hi),x,n);
953  if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -((double)n+0.1)*DBL_EPSILON )
954  return false;
955  return true;
956 }
957 inline bool fp_bound_tol( sys_float lo, sys_float x, sys_float hi, sys_float tol )
958 {
959  ASSERT( tol > 0.f );
960  // mimic IEEE behavior
961  if( isnan(x) || isnan(lo) || isnan(hi) )
962  return false;
963  if( fp_equal_tol(lo,hi,tol) )
964  return fp_equal_tol(0.5f*(lo+hi),x,tol);
965  if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -tol )
966  return false;
967  return true;
968 }
969 inline bool fp_bound_tol( double lo, double x, double hi, double tol )
970 {
971  ASSERT( tol > 0. );
972  // mimic IEEE behavior
973  if( isnan(x) || isnan(lo) || isnan(hi) )
974  return false;
975  if( fp_equal_tol(lo,hi,tol) )
976  return fp_equal_tol(0.5*(lo+hi),x,tol);
977  if( ((hi-x)/(hi-lo))*((x-lo)/(hi-lo)) < -tol )
978  return false;
979  return true;
980 }
981 
982 
983 #undef POW2
984 
985 #define POW2 pow2
986 template<class T>
987 inline T pow2(T a) { return a*a; }
988 /* */
989 
990 #undef POW3
991 
992 #define POW3 pow3
993 template<class T>
994 inline T pow3(T a) { return a*a*a; }
995 /* */
996 
997 #undef POW4
998 
999 #define POW4 pow4
1000 template<class T>
1001 inline T pow4(T a) { T b = a*a; return b*b; }
1002 /* */
1003 
1004 #undef SDIV
1005 
1008 inline sys_float SDIV( sys_float x ) { return ( fabs((double)x) < (double)SMALLFLOAT ) ? (sys_float)SMALLFLOAT : x; }
1009 /* \todo should we use SMALLDOUBLE here ? it produces overflows now... PvH */
1010 inline double SDIV( double x ) { return ( fabs(x) < (double)SMALLFLOAT ) ? (double)SMALLFLOAT : x; }
1011 // inline double SDIV( double x ) { return ( fabs(x) < SMALLDOUBLE ) ? SMALLDOUBLE : x; }
1012 /* */
1013 
1018 {
1019  // this should crash...
1020  if( isnan(x) || isnan(y) )
1021  return x/y;
1022  int sx = sign3(x);
1023  int sy = sign3(y);
1024  // 0/0 -> NaN, this should crash as well...
1025  if( sx == 0 && sy == 0 )
1026  {
1027  if( isnan(res_0by0) )
1028  return x/y;
1029  else
1030  return res_0by0;
1031  }
1032  if( sx == 0 )
1033  return 0.;
1034  if( sy == 0 )
1035  return ( sx < 0 ) ? -FLT_MAX : FLT_MAX;
1036  // at this stage x != 0. and y != 0.
1037  sys_float ay = abs(y);
1038  if( ay >= 1.f )
1039  return x/y;
1040  else
1041  {
1042  // multiplication is safe since ay < 1.
1043  if( abs(x) < ay*FLT_MAX )
1044  return x/y;
1045  else
1046  return ( sx*sy < 0 ) ? -FLT_MAX : FLT_MAX;
1047  }
1048 }
1049 
1051 {
1052  return safe_div( x, y, numeric_limits<sys_float>::quiet_NaN() );
1053 }
1054 
1058 inline double safe_div(double x, double y, double res_0by0)
1059 {
1060  // this should crash...
1061  if( isnan(x) || isnan(y) )
1062  return x/y;
1063  int sx = sign3(x);
1064  int sy = sign3(y);
1065  // 0/0 -> NaN, this should crash as well...
1066  if( sx == 0 && sy == 0 )
1067  {
1068  if( isnan(res_0by0) )
1069  return x/y;
1070  else
1071  return res_0by0;
1072  }
1073  if( sx == 0 )
1074  return 0.;
1075  if( sy == 0 )
1076  return ( sx < 0 ) ? -DBL_MAX : DBL_MAX;
1077  // at this stage x != 0. and y != 0.
1078  double ay = abs(y);
1079  if( ay >= 1. )
1080  return x/y;
1081  else
1082  {
1083  // multiplication is safe since ay < 1.
1084  if( abs(x) < ay*DBL_MAX )
1085  return x/y;
1086  else
1087  return ( sx*sy < 0 ) ? -DBL_MAX : DBL_MAX;
1088  }
1089 }
1090 
1091 inline double safe_div(double x, double y)
1092 {
1093  return safe_div( x, y, numeric_limits<double>::quiet_NaN() );
1094 }
1095 
1096 template<class T>
1097 inline void invalidate_array(T* p, size_t size)
1098 {
1099  if( size > 0 )
1100  memset( p, -1, size );
1101 }
1102 
1103 inline void invalidate_array(double* p, size_t size)
1104 {
1105  set_NaN( p, (long)(size/sizeof(double)) );
1106 }
1107 
1108 inline void invalidate_array(sys_float* p, size_t size)
1109 {
1110  set_NaN( p, (long)(size/sizeof(sys_float)) );
1111 }
1112 
1115 template<class T> inline T* get_ptr(T *v)
1116 {
1117  return v;
1118 }
1119 template<class T> inline T* get_ptr(valarray<T> &v)
1120 {
1121  return &v[0];
1122 }
1123 template<class T, class U> inline T* get_ptr(vector<T,U> &v)
1124 {
1125  return &v[0];
1126 }
1127 template<class T> inline const T* get_ptr(const valarray<T> &v)
1128 {
1129  return const_cast<const T*>(&const_cast<valarray<T>&>(v)[0]);
1130 }
1131 template<class T, class U> inline const T* get_ptr(const vector<T,U> &v)
1132 {
1133  return const_cast<const T*>(&const_cast<vector<T,U>&>(v)[0]);
1134 }
1135 
1141 double csphot(long int inu, long int ithr, long int iofset);
1142 
1144 double AnuUnit(realnum energy);
1145 
1150 void cap4(char *chCAP , const char *chLab);
1151 
1154 void uncaps(char *chCard);
1155 void uncaps(string& chCard);
1156 
1159 void caps(char *chCard);
1160 void caps(string& chCard);
1161 
1168 double FFmtRead(const char *chCard,
1169  long int *ipnt,
1170  long int last,
1171  bool *lgEOL);
1172 
1178 long nMatch(const char *chKey,
1179  const char *chCard);
1180 
1181 // these are safe versions of strstr, strchr, etc to work around a deficiency in glibc
1182 inline const char *strstr_s(const char *haystack, const char *needle)
1183 {
1184  return const_cast<const char *>(strstr(haystack, needle));
1185 }
1186 
1187 inline char *strstr_s(char *haystack, const char *needle)
1188 {
1189  return const_cast<char *>(strstr(haystack, needle));
1190 }
1191 
1192 inline const char *strchr_s(const char *s, int c)
1193 {
1194  return const_cast<const char *>(strchr(s, c));
1195 }
1196 
1197 inline char *strchr_s(char *s, int c)
1198 {
1199  return const_cast<char *>(strchr(s, c));
1200 }
1201 
1204 long int ipow( long, long );
1205 
1206 // this routine is fully equivalent to snprintf() apart from the fact that it
1207 // concatenates the output to an existing string rather than replacing it.
1208 size_t sncatf( char* buf, size_t bufSize, const char* fmt, ... );
1209 
1210 size_t sncatf( ostringstream& buf, const char* fmt, ... );
1211 
1214 void PrintE82( FILE*, double );
1215 
1217 void PrintE71( FILE*, double );
1218 
1220 void PrintE93( FILE*, double );
1221 
1227 // prevent compiler warnings on non-MS systems
1228 #ifdef _MSC_VER
1229 char *PrintEfmt(const char *fmt, double val );
1230 #else
1231 #define PrintEfmt( F, V ) F, V
1232 #endif
1233 
1235 const double SEXP_LIMIT = 84.;
1237 const double DSEXP_LIMIT = 680.;
1238 
1241 double sexp(double x);
1242 
1247 double dsexp(double x);
1248 
1250 inline double exp10(double x)
1251 {
1252  if( x < -330. )
1253  return 0.;
1254  else if( x > 310. )
1255  return pow2(BIGDOUBLE); // +Inf
1256  else
1257  {
1258  double y = trunc(x/3.);
1259  double z = 3.*y;
1260  x -= z;
1261  x *= 2.302585092994045684; // constant is ln(10)
1262  x -= 7.90550887243868070612e-3*z; // constant is ln(10)-10/3*ln(2)
1263  return ldexp(exp(x),10*int(y));
1264  }
1265 }
1266 
1269 {
1270  if( x < -50.f )
1271  return 0.f;
1272  else if( x > 40.f )
1273  return pow2(BIGFLOAT); // +Inf
1274  else
1275  {
1276  sys_float y = truncf(x/3.f);
1277  sys_float z = 3.f*y;
1278  x -= z;
1279  x *= 2.302585093f; // constant is ln(10)
1280  x -= 7.9055088724e-3f*z; // constant is ln(10)-10/3*ln(2)
1281  return ldexpf(expf(x),10*int(y));
1282  }
1283 }
1284 
1286 {
1287  return exp10f(x);
1288 }
1289 
1294 double plankf(long int ip);
1295 
1296 // safe version of getline() that correctly handles all types of EOL lf, crlf and cr...
1297 istream& SafeGetline(istream& is, string& t);
1298 
1308 void spsort( realnum x[], long int n, long int iperm[], int kflag, int *ier);
1309 
1310 /**************************************************************************
1311  *
1312  * disable some bogus errors in the ms c compiler
1313  *
1314  **************************************************************************/
1315 
1316 /* */
1317 #ifdef _MSC_VER
1318  /* disable strcat warning */
1319 # pragma warning( disable : 4996 )
1320  /* disable bogus underflow warning in MS VS*/
1321 # pragma warning( disable : 4056 )
1322  /* disable "inline function removed since not used", MS VS*/
1323 # pragma warning( disable : 4514 )
1324  /* disable "assignment operator could not be generated", cddefines.h
1325  * line 126 */
1326 # pragma warning( disable : 4512 )
1327 #endif
1328 #ifdef __INTEL_COMPILER
1329 # pragma warning( disable : 1572 )
1330 #endif
1331 /* */
1332 
1333 /*lint +e129 these resolve several issues pclint has with my system headers */
1334 /*lint +e78 */
1335 /*lint +e830 */
1336 /*lint +e38 */
1337 /*lint +e148 */
1338 /*lint +e114 */
1339 /*lint +e18 */
1340 /*lint +e49 */
1341 
1342 #endif /* CDDEFINES_H_ */
1343 
const char * p_name
Definition: cddefines.h:708
T pow4(T a)
Definition: cddefines.h:1001
bool lgTestCodeCalled
Definition: cddefines.cpp:10
FILE * ioMAP
Definition: cdinit.cpp:9
void leave(const char *name)
Definition: cddefines.h:686
const char * file() const
Definition: cddefines.h:599
const int ipBERYLLIUM
Definition: cddefines.h:362
void leave(const char *) const
Definition: cddefines.h:702
FILE * p_fp
Definition: cddefines.h:673
t_debug()
Definition: cddefines.h:676
void PrintE93(FILE *, double)
Definition: service.cpp:1126
const int ipMAGNESIUM
Definition: cddefines.h:370
int p_callLevel
Definition: cddefines.h:674
bool is_odd(int j)
Definition: cddefines.h:757
Definition: cddefines.h:447
Definition: cddefines.h:417
virtual ~bad_signal()
Definition: cddefines.h:577
exit_type exit_status() const
Definition: cddefines.h:473
double exp10(double x)
Definition: cddefines.h:1250
const int FILENAME_PATH_LENGTH_2
Definition: cddefines.h:306
#define NORETURN
Definition: cpu.h:461
Definition: cddefines.h:213
istream & SafeGetline(istream &is, string &t)
Definition: service.cpp:1849
FILE * ioStdin
Definition: cddefines.cpp:8
T * get_ptr(T *v)
Definition: cddefines.h:1115
void cdBacktrace()
Definition: cddefines.h:516
NORETURN void TotalInsanity(void)
Definition: service.cpp:1174
const double BIGDOUBLE
Definition: cpu.h:238
bool read_whole_line(string &chLine, FILE *ioIN)
Definition: service.cpp:68
void set_NaN(sys_float &x)
Definition: cpu.cpp:871
const int ipARGON
Definition: cddefines.h:376
const int ipTITANIUM
Definition: cddefines.h:380
const char * file() const
Definition: cddefines.h:465
const realnum SMALLFLOAT
Definition: cpu.h:235
~debugtrace()
Definition: cddefines.h:717
t_cpu_i & i()
Definition: cpu.h:414
FILE * sys_fopen(const char *path, const char *mode)
Definition: cddefines.h:145
const int NISO
Definition: cddefines.h:321
size_t sncatf(char *buf, size_t bufSize, const char *fmt,...)
Definition: service.cpp:919
const char * routine() const
Definition: cddefines.h:461
Definition: cddefines.h:160
FILE * ioPrnErr
Definition: cddefines.cpp:9
char TorF(bool l)
Definition: cddefines.h:753
const int ipOXYGEN
Definition: cddefines.h:366
const int ipCHLORINE
Definition: cddefines.h:375
#define PrintEfmt(F, V)
Definition: cddefines.h:1231
void enter(const char *) const
Definition: cddefines.h:701
long int nzone
Definition: cddefines.cpp:13
void invalidate_array(T *p, size_t size)
Definition: cddefines.h:1097
long nMatch(const char *chKey, const char *chCard)
Definition: service.cpp:728
Definition: cddefines.h:706
const void * ZeroPtr
Definition: cdinit.cpp:15
T sign(T x, T y)
Definition: cddefines.h:852
int sig() const
Definition: cddefines.h:578
const char * name() const
Definition: cddefines.h:721
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
Definition: cddefines.h:910
T TotalInsanityAsStub()
Definition: cddefines.h:495
const char * strstr_s(const char *haystack, const char *needle)
Definition: cddefines.h:1182
void GenerateBacktrace(void *ptr)
Definition: cpu.cpp:1178
sys_float sexp(sys_float x)
Definition: service.cpp:1203
const int ipRecNetEsc
Definition: cddefines.h:341
long line() const
Definition: cddefines.h:603
int dbg_printf(int debug, const char *fmt,...)
Definition: service.cpp:1357
T pow3(T a)
Definition: cddefines.h:994
double fudge(long int ipnt)
Definition: service.cpp:758
const int ipCOBALT
Definition: cddefines.h:385
int p_sig
Definition: cddefines.h:572
NORETURN void OUT_OF_RANGE(const char *str)
Definition: cddefines.h:646
bool lgTestCodeEnabled
Definition: cddefines.cpp:11
const int ipNICKEL
Definition: cddefines.h:386
void FPRead(istringstream &iss, const string &s, double &value)
Definition: service.cpp:548
const int ipZINC
Definition: cddefines.h:388
exit_type
Definition: cddefines.h:151
void cap4(char *chCAP, const char *chLab)
Definition: service.cpp:257
Output(FILE *fp)
Definition: cddefines.h:239
const char * p_file
Definition: cddefines.h:586
double dsexp(double x)
Definition: service.cpp:1242
void CodeReview(void)
Definition: service.cpp:1297
const int ipSULPHUR
Definition: cddefines.h:374
Definition: cddefines.h:155
Definition: cddefines.h:164
bool fp_bound_tol(sys_float lo, sys_float x, sys_float hi, sys_float tol)
Definition: cddefines.h:957
void fixit_base(const char *func, const char *file, int line, const char *reason)
Definition: service.cpp:1280
const double MAX_DENSITY
Definition: cddefines.h:329
Definition: cddefines.h:335
Definition: cddefines.h:234
static T & Inst()
Definition: cddefines.h:216
bool lgPrnErr
Definition: cddefines.cpp:12
Definition: cddefines.h:158
void uncaps(char *chCard)
Definition: service.cpp:277
char toupper(char c)
Definition: cddefines.h:743
void print(void) const
Definition: cddefines.h:594
#define fopen
Definition: cddefines.h:149
const double DSEXP_LIMIT
Definition: cddefines.h:1237
void PrintBacktrace(const string &s, bool lgPrintSeed=true)
Definition: cpu.cpp:1185
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:864
void broken(void)
Definition: service.cpp:1271
debugtrace(const char *funcname)
Definition: cddefines.h:710
int dprintf(FILE *fp, const char *format,...)
Definition: service.cpp:1306
const int ipVANADIUM
Definition: cddefines.h:381
const int ipIRON
Definition: cddefines.h:384
void PrintE71(FILE *, double)
Definition: service.cpp:1076
Definition: cddefines.h:156
sys_float exp10f(sys_float x)
Definition: cddefines.h:1268
char tolower(char c)
Definition: cddefines.h:734
double energy(const genericState &gs)
Definition: generic_state.cpp:51
Definition: cddefines.h:159
const char * p_routine
Definition: cddefines.h:449
Definition: cddefines.h:152
const int ipSCANDIUM
Definition: cddefines.h:379
Definition: cddefines.h:695
Fixit(const char *func, const char *file, int line, const char *reason)
Definition: cddefines.h:420
const char * p_comment
Definition: cddefines.h:615
long int iteration
Definition: cddefines.cpp:15
const double ZeroNum
Definition: cdinit.cpp:13
const int ipCRD
Definition: cddefines.h:352
void enter(const char *name)
Definition: cddefines.h:681
FILE * m_fp
Definition: cddefines.h:237
virtual ~bad_assert()
Definition: cddefines.h:593
long p_line
Definition: cddefines.h:451
float realnum
Definition: cddefines.h:127
float sys_float
Definition: cddefines.h:130
const int ipLY_A
Definition: cddefines.h:356
const realnum BIGFLOAT
Definition: cpu.h:233
const int ipPHOSPHORUS
Definition: cddefines.h:373
const int INPUT_LINE_LENGTH
Definition: cddefines.h:311
long max(int a, long b)
Definition: cddefines.h:821
NORETURN void DOMAIN_ERROR(const string &str)
Definition: cddefines.h:652
int sign3(T x)
Definition: cddefines.h:860
const char * p_file
Definition: cddefines.h:450
const int ipNEON
Definition: cddefines.h:368
const int ipCHROMIUM
Definition: cddefines.h:382
double powi(double, long int)
Definition: service.cpp:797
const int ipRecRad
Definition: cddefines.h:343
long min(int a, long b)
Definition: cddefines.h:766
bool fp_bound(sys_float lo, sys_float x, sys_float hi, int n=3)
Definition: cddefines.h:933
long p_line
Definition: cddefines.h:587
const char * strchr_s(const char *s, int c)
Definition: cddefines.h:1192
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition: cddefines.h:1017
const int ipRecEsc
Definition: cddefines.h:339
virtual ~cloudy_abort()
Definition: cddefines.h:620
#define NULL
Definition: cddefines.h:115
Definition: cddefines.h:154
Definition: cddefines.h:166
exit_type p_exit
Definition: cddefines.h:452
const int ipKRYPTON
Definition: cddefines.h:389
double AnuUnit(realnum energy)
Definition: service.cpp:189
void MyAssert(const char *file, int line, const char *comment)
Definition: service.cpp:172
Definition: cddefines.h:334
FILE * ioQQQ
Definition: cddefines.cpp:7
const int ipSILICON
Definition: cddefines.h:372
const char * comment() const
Definition: cddefines.h:607
long int ipow(long, long)
Definition: service.cpp:861
#define ASSERT(exp)
Definition: cddefines.h:637
FILE * fptr() const
Definition: cddefines.h:240
const int ipALUMINIUM
Definition: cddefines.h:371
const char * p_comment
Definition: cddefines.h:588
const int ipCALCIUM
Definition: cddefines.h:378
const int ipNITROGEN
Definition: cddefines.h:365
double csphot(long int inu, long int ithr, long int iofset)
Definition: service.cpp:1768
const int LIMELM
Definition: cddefines.h:318
T pow2(T a)
Definition: cddefines.h:987
const int ipFLUORINE
Definition: cddefines.h:367
double powpq(double x, int p, int q)
Definition: service.cpp:833
#define isnan
Definition: cddefines.h:666
const int ipHELIUM
Definition: cddefines.h:360
Definition: cddefines.h:163
const int ipMANGANESE
Definition: cddefines.h:383
const char * comment() const
Definition: cddefines.h:621
double fnzone
Definition: cddefines.cpp:14
long line() const
Definition: cddefines.h:469
Definition: cddefines.h:584
void vzero(vector< T > &vec)
Definition: cddefines.h:842
Definition: cddefines.h:157
const int NCHLAB
Definition: cddefines.h:314
const int FILENAME_PATH_LENGTH
Definition: cddefines.h:303
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1325
Definition: cddefines.h:570
Definition: cddefines.h:613
Definition: cddefines.h:161
Definition: cddefines.h:162
const int ipPRD
Definition: cddefines.h:350
sys_float SDIV(sys_float x)
Definition: cddefines.h:1008
double pow(double x, int i)
Definition: cddefines.h:782
const int ipCARBON
Definition: cddefines.h:364
cloudy_exit(const char *routine, const char *file, long line, exit_type exit_code)
Definition: cddefines.h:454
long nint(double x)
Definition: cddefines.h:762
void caps(char *chCard)
Definition: service.cpp:299
static t_cpu cpu
Definition: cpu.h:422
void ShowMe(void)
Definition: service.cpp:197
void PrintE82(FILE *, double)
Definition: service.cpp:1027
Definition: cddefines.h:670
const double SEXP_LIMIT
Definition: cddefines.h:1235
Definition: cddefines.h:165
double plankf(long int ip)
Definition: service.cpp:1786
const int ipHYDROGEN
Definition: cddefines.h:359
const int ipLITHIUM
Definition: cddefines.h:361
void cdPrepareExit(exit_type)
Definition: cdinit.cpp:131
const double DEPTH_OFFSET
Definition: cddefines.h:332
void TestCode(void)
Definition: service.cpp:1261
Definition: cddefines.h:167
const int ipPOTASSIUM
Definition: cddefines.h:377
const int NHYDRO_MAX_LEVEL
Definition: cddefines.h:326
NORETURN void BadRead(void)
Definition: service.cpp:1190
const int ipBORON
Definition: cddefines.h:363
const int ipCRDW
Definition: cddefines.h:354
const int ipSODIUM
Definition: cddefines.h:369
Definition: cddefines.h:153
void spsort(realnum x[], long int n, long int iperm[], int kflag, int *ier)
Definition: service.cpp:1426
t_nodebug()
Definition: cddefines.h:699
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:393
const int ipCOPPER
Definition: cddefines.h:387