00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include <stdio.h>
00013 #include <stdlib.h>
00014 #include <math.h>
00015 #include <assert.h>
00016 #include <string.h>
00017 #include <float.h>
00018 #include <limits.h>
00019 #include <time.h>
00020 #include <string>
00021 #include <vector>
00022 #include <valarray>
00023 #include <complex>
00024 #ifdef DMALLOC
00025 #include <dmalloc.h>
00026 #endif
00027
00028
00029
00030
00031
00032
00033
00034
00035 using namespace std;
00036
00038 #ifndef EXTERN
00039 #define EXTERN extern
00040 #endif
00041
00042 #undef STATIC
00043 #ifdef USE_GPROF
00044 #define STATIC
00045 #else
00046 #define STATIC static
00047 #endif
00048
00049
00050 #include "cpu.h"
00051
00052
00053
00054
00072
00073
00074
00075
00076
00077 template<typename T> class Singleton
00078 {
00079 public:
00080 static T& Inst()
00081 {
00082 static T instance__;
00083 return instance__;
00084 }
00085 };
00086
00087
00088
00089
00090
00091
00092
00093
00094
00098 EXTERN FILE *ioQQQ;
00099
00100 EXTERN FILE *ioStdin;
00101
00102 extern FILE *ioMAP;
00103
00106 EXTERN FILE* ioPrnErr;
00107
00109 EXTERN bool lgAbort;
00110
00113 EXTERN bool lgTestCodeCalled;
00114
00117 EXTERN bool lgTestCodeEnabled;
00118
00121 EXTERN bool lgPrnErr;
00122
00125 EXTERN bool lgAssertFPE;
00126
00129 EXTERN long int nzone;
00130
00132 EXTERN double fnzone;
00133
00136 EXTERN long int iteration;
00137
00144 extern const double ZeroNum;
00145
00146
00147
00148
00149
00150
00151
00152
00157 const int FILENAME_PATH_LENGTH = 200;
00158
00160 const int FILENAME_PATH_LENGTH_2 = FILENAME_PATH_LENGTH*2;
00161
00165 const int INPUT_LINE_LENGTH = 200;
00166
00169 const int LIMELM = 30;
00170
00172 const int NISO = 2;
00173
00177 const int NHYDRO_MAX_LEVEL = 401;
00178
00181
00182
00184 const int N_H_MOLEC = 8;
00185
00187 const double TEND = 4000.;
00188
00190 const double DEPTH_OFFSET = 1.e-30;
00191
00192
00193 const int ipELECTRON = 0;
00194 const int ipPROTON = 1;
00195 const int ipHE_PLUS = 2;
00196
00197
00198
00199 const int ipRecEsc = 2;
00200
00201 const int ipRecNetEsc = 1;
00202
00203 const int ipRecRad = 0;
00208
00209
00210 const int ipPRD = 1;
00211
00212 const int ipCRD = -1;
00213
00214 const int ipCRDW = 2;
00215
00216 const int ipLY_A = -2;
00217
00218 const int ipDEST_K2 = 1;
00219
00220 const int ipDEST_INCOM = 2;
00221
00222 const int ipDEST_SIMPL = 3;
00223
00226 const int ipH1s = 0;
00227 const int ipH2s = 1;
00228 const int ipH2p = 2;
00229
00232
00233 const int ipHe1s1S = 0;
00234
00235
00236 const int ipHe2s3S = 1;
00237 const int ipHe2s1S = 2;
00238 const int ipHe2p3P0 = 3;
00239 const int ipHe2p3P1 = 4;
00240 const int ipHe2p3P2 = 5;
00241 const int ipHe2p1P = 6;
00242
00243
00244 const int ipHe3s3S = 7;
00245 const int ipHe3s1S = 8;
00246 const int ipHe3p3P = 9;
00247 const int ipHe3d3D = 10;
00248 const int ipHe3d1D = 11;
00249 const int ipHe3p1P = 12;
00250
00251
00252 const int ipHe4s3S = 13;
00253 const int ipHe4s1S = 14;
00254 const int ipHe4p3P = 15;
00255 const int ipHe4d3D = 16;
00256 const int ipHe4d1D = 17;
00257 const int ipHe4f3F = 18;
00258 const int ipHe4f1F = 19;
00259 const int ipHe4p1P = 20;
00260
00261
00262 const int ipHe5s3S = 21;
00263 const int ipHe5s1S = 22;
00264 const int ipHe5p3P = 23;
00265 const int ipHe5d3D = 24;
00266 const int ipHe5d1D = 25;
00267 const int ipHe5f3F = 26;
00268 const int ipHe5f1F = 27;
00269 const int ipHe5g3G = 28;
00270 const int ipHe5g1G = 29;
00271 const int ipHe5p1P = 30;
00272
00276 const int ipH_LIKE = 0;
00277 const int ipHE_LIKE = 1;
00278 const int ipLI_LIKE = 2;
00279 const int ipBE_LIKE = 3;
00280 const int ipB_LIKE = 4;
00281 const int ipC_LIKE = 5;
00282 const int ipN_LIKE = 6;
00283 const int ipO_LIKE = 7;
00284 const int ipF_LIKE = 8;
00285
00287 const int ipHYDROGEN = 0;
00288 const int ipHELIUM = 1;
00289 const int ipLITHIUM = 2;
00290 const int ipBERYLLIUM = 3;
00291 const int ipBORON = 4;
00292 const int ipCARBON = 5;
00293 const int ipNITROGEN = 6;
00294 const int ipOXYGEN = 7;
00295 const int ipFLUORINE = 8;
00296 const int ipNEON = 9;
00297 const int ipSODIUM = 10;
00298 const int ipMAGNESIUM = 11;
00299 const int ipALUMINIUM = 12;
00300 const int ipSILICON = 13;
00301 const int ipPHOSPHORUS = 14;
00302 const int ipSULPHUR = 15;
00303 const int ipCHLORINE = 16;
00304 const int ipARGON = 17;
00305 const int ipPOTASSIUM = 18;
00306 const int ipCALCIUM = 19;
00307 const int ipSCANDIUM = 20;
00308 const int ipTITANIUM = 21;
00309 const int ipVANADIUM = 22;
00310 const int ipCHROMIUM = 23;
00311 const int ipMANGANESE = 24;
00312 const int ipIRON = 25;
00313 const int ipCOBALT = 26;
00314 const int ipNICKEL = 27;
00315 const int ipCOPPER = 28;
00316 const int ipZINC = 29;
00317
00318
00319
00320
00321
00322
00323
00324
00328 #ifndef NDEBUG
00329 # define DEBUG
00330 #else
00331 # undef DEBUG
00332 #endif
00333
00341 #ifndef PARALLEL_MODE
00342 #define PARALLEL_MODE false
00343 #endif
00344
00346 #if defined(malloc)
00347
00348
00349 # define MALLOC(exp) (malloc(exp))
00350 #else
00351
00352 # define MALLOC(exp) (MyMalloc(exp,__FILE__, __LINE__))
00353 #endif
00354
00356 #if defined(calloc)
00357
00358 # define CALLOC calloc
00359 #else
00360
00361 # define CALLOC MyCalloc
00362 #endif
00363
00365 #if defined(realloc)
00366
00367 # define REALLOC realloc
00368 #else
00369
00370 # define REALLOC MyRealloc
00371 #endif
00372
00374 #undef ASSERT
00375 #ifndef STD_ASSERT
00376 # ifdef NDEBUG
00377 # define ASSERT(exp) ((void)0)
00378 # else
00379
00380 # define ASSERT(exp) if (!(exp)) \
00381 MyAssert(__FILE__, __LINE__)
00382 # endif
00383 #else
00384 # define ASSERT(exp) (assert(exp))
00385 #endif
00386
00387
00388
00389
00390 #undef isnan
00391 #define isnan MyIsnan
00392
00395 class t_debugprt : public Singleton<t_debugprt>
00396 {
00397 friend class Singleton<t_debugprt>;
00398 protected:
00399 t_debugprt()
00400 {
00401 indent = 0;
00402 debug_fp = stderr;
00403 }
00404 private:
00405 int indent;
00406 FILE *debug_fp;
00407 void print_str( const char *str, char dir ) const
00408 {
00409 for( int i=0; i < indent; i++ )
00410 fputc( ' ', debug_fp );
00411 fputc( dir, debug_fp );
00412 fputs( str, debug_fp );
00413 fputc( '\n', debug_fp );
00414 }
00415 public:
00416 void printin( const char *str )
00417 {
00418 print_str( str, '>' );
00419 indent++;
00420 }
00421 void printout( const char *str )
00422 {
00423 --indent;
00424 print_str( str, '<' );
00425 }
00426 };
00427
00428 #ifdef DEBUG_FUN
00429 #define DEBUG_ENTRY( STR ) t_debugprt::Inst().printin( STR )
00430 #define DEBUG_EXIT( STR ) t_debugprt::Inst().printout( STR )
00431 #else
00432 #define DEBUG_ENTRY( STR ) ((void)0)
00433 #define DEBUG_EXIT( STR ) ((void)0)
00434 #endif
00435
00436
00437 inline char TorF( bool l ) { return l ? 'T' : 'F'; }
00438
00439
00441 inline bool is_odd( int j ) { return (j&1) == 1; }
00442 inline bool is_odd( long j ) { return (j&1L) == 1L; }
00443
00444
00446 inline long nint( double x ) { return static_cast<long>( (x < 0.) ? x-0.5 : x+0.5 ); }
00447
00448
00449
00450 inline long min( int a, long b ) { long c = a; return ( (c < b) ? c : b ); }
00451 inline long min( long a, int b ) { long c = b; return ( (a < c) ? a : c ); }
00452 inline double min( float a, double b ) { double c = a; return ( (c < b) ? c : b ); }
00453 inline double min( double a, float b ) { double c = b; return ( (a < c) ? a : c ); }
00454
00455 #undef MIN2
00456
00457
00458 #define MIN2 min
00459
00460
00461 #undef MIN3
00462
00463 #define MIN3(a,b,c) (MIN2(MIN2(a,b),c))
00464
00465
00466 #undef MIN4
00467
00468 #define MIN4(a,b,c,d) (MIN2(MIN2(a,b),MIN2(c,d)))
00469
00470
00471
00472 inline long max( int a, long b ) { long c = a; return ( (c > b) ? c : b ); }
00473 inline long max( long a, int b ) { long c = b; return ( (a > c) ? a : c ); }
00474 inline double max( float a, double b ) { double c = a; return ( (c > b) ? c : b ); }
00475 inline double max( double a, float b ) { double c = b; return ( (a > c) ? a : c ); }
00476
00477 #undef MAX2
00478
00479
00480 #define MAX2 max
00481
00482
00483 #undef MAX3
00484
00485 #define MAX3(a,b,c) (MAX2(MAX2(a,b),c))
00486
00487
00488 #undef MAX4
00489
00490 #define MAX4(a,b,c,d) (MAX2(MAX2(a,b),MAX2(c,d)))
00491
00492
00497 template<class T>
00498 inline T sign( T x, T y )
00499 {
00500 return ( y < T() ) ? -abs(x) : abs(x);
00501 }
00502
00503
00505 template<class T>
00506 inline int sign3( T x ) { return ( x < T() ) ? -1 : ( ( x > T() ) ? 1 : 0 ); }
00507
00508
00510 inline bool fp_equal( float x, float y )
00511 {
00512
00513 if( isnan(x) || isnan(y) )
00514 return false;
00515 int sx = sign3(x);
00516 int sy = sign3(y);
00517
00518 if( sx == 0 && sy == 0 )
00519 return true;
00520
00521 if( sx*sy != 1 )
00522 return false;
00523 x = abs(x);
00524 y = abs(y);
00525 return ( 1.f - min(x,y)/max(x,y) < 3.1f*FLT_EPSILON );
00526 }
00527
00528 inline bool fp_equal( double x, double y )
00529 {
00530
00531 if( isnan(x) || isnan(y) )
00532 return false;
00533 int sx = sign3(x);
00534 int sy = sign3(y);
00535
00536 if( sx == 0 && sy == 0 )
00537 return true;
00538
00539 if( sx*sy != 1 )
00540 return false;
00541 x = abs(x);
00542 y = abs(y);
00543 return ( 1. - min(x,y)/max(x,y) < 3.1*DBL_EPSILON );
00544 }
00545
00546 #undef POW2
00547
00548
00549 #define POW2 pow2
00550 inline int pow2(int a) { return a*a; }
00551 inline long pow2(long a) { return a*a; }
00552 inline float pow2(float a) { return a*a; }
00553 inline double pow2(double a) { return a*a; }
00554
00555
00556 #undef POW3
00557
00558
00559 #define POW3 pow3
00560 inline int pow3(int a) { return a*a*a; }
00561 inline long pow3(long a) { return a*a*a; }
00562 inline float pow3(float a) { return a*a*a; }
00563 inline double pow3(double a) { return a*a*a; }
00564
00565
00566 #undef POW4
00567
00568 #define POW4(X) (POW2(X)*POW2(X))
00569
00570
00571 #undef SDIV
00572
00575
00576 inline float SDIV( float x ) { return ( fabs((double)x) < (double)SMALLFLOAT ) ? SMALLFLOAT : x; }
00577
00578 inline double SDIV( double x ) { return ( fabs(x) < (double)SMALLFLOAT ) ? (double)SMALLFLOAT : x; }
00579
00580
00581
00590
00591
00592
00593
00594
00595 #ifndef _ERR_
00596 #define _ERR_ -1
00597 #endif
00598
00599
00600 #undef HMRATE
00601
00602
00603
00604
00605
00606 #define HMRATE(a,b,c) hmrate4(a,b,c,phycon.te)
00607
00608 inline double hmrate4( double a, double b, double c, double te )
00609 {
00610 if( b == 0. && c == 0. )
00611 return a;
00612 else if( c == 0. )
00613 return a*pow(te/300.,b);
00614 else if( b == 0. )
00615 return ( c/te <= 50. ) ? a*exp(-c/te) : 0.;
00616 else
00617 return ( c/te <= 50. ) ? a*pow(te/300.,b)*exp(-c/te) : 0.;
00618 }
00619
00620
00621
00622
00623
00624
00625
00626 typedef struct
00627 {
00636 int iRedisFun;
00637
00642 int ipCont;
00643
00645 long int ipFine;
00646
00648 int IonStg;
00649
00651 int nelem;
00652
00654 float ColUL;
00655
00669 float TauIn;
00670
00681 float TauTot;
00682
00688 float TauCon;
00689
00692 float FracInwd;
00693
00696 double pump;
00697
00699 double xIntensity;
00700
00702 double phots;
00703
00705 float gf;
00706
00709 float Pesc;
00710
00713 float Pelec_esc;
00714
00717 float Pdest;
00718
00722 float dampXvel;
00723
00725 float damp;
00726
00728 double cool , heat;
00729
00731 float ColOvTot;
00732
00735 float cs;
00736
00738 float WLAng;
00739
00741 float EnergyK;
00742
00744 float EnergyErg;
00745
00747 float EnergyWN;
00748
00754 float opacity;
00755
00757 float gLo , gHi;
00758
00761 double PopLo , PopHi;
00762
00764 double PopOpc;
00765
00768 float Aul;
00769
00771 double ots;
00772
00773 } EmLine;
00774
00775
00776
00777
00778
00779
00780
00781
00786 double fudge(long int ipnt);
00787
00791 void broken(void);
00792
00795 void fixit(void);
00796
00798 void CodeReview(void);
00799
00801 void TestCode(void);
00802
00808 void *MyMalloc(size_t size, const char *file, int line);
00809
00814 void *MyCalloc(size_t num, size_t size);
00815
00820 void *MyRealloc(void *p, size_t size);
00821
00826 NORETURN void MyAssert(const char *file, int line);
00827
00829 NORETURN void cdEXIT(int);
00830
00832 void path_not_set( void );
00833
00835 void ShowMe(void);
00836
00838 NORETURN void TotalInsanity(void);
00839
00841 NORETURN void BadRead(void);
00842
00844 NORETURN void BadOpen(void);
00845
00850 NORETURN void NoNumb(char *chCard);
00851
00855 void dbg_printf(int debug, const char *fmt, ...);
00856
00857
00858
00859
00860
00861
00862
00863
00864 void cdDefines(void);
00865
00871 double csphot(long int inu, long int ithr, long int iofset);
00872
00877 double RandGauss(double xMean, double s );
00878
00882 double MyGaussRand( double PctUncertainty );
00883
00885 double AnuUnit(float energy);
00886
00891 void cap4(char *chCAP , char *chLab);
00892
00895 void caps(char *chCard );
00896
00899 double e2(
00900 double t );
00901
00904 double ee1(double x);
00905
00909 double ee1_safe(double x);
00910
00917 double FFmtRead(const char *chCard,
00918 long int *ipnt,
00919 long int last,
00920 bool *lgEOL);
00921
00927 long nMatch(const char *chKey,
00928 const char *chCard);
00929
00935 long int GetElem( char *chCARD_CAPS );
00936
00945 int GetQuote(char *chLabel, char *chCard, bool lgABORT );
00946
00947
00948 #if !HAVE_POWI
00949
00950 double powi( double , long int );
00951 #endif
00952
00955 long int ipow( long, long );
00956
00959 void PrintE82( FILE*, double );
00960
00962 void PrintE71( FILE*, double );
00963
00965 void PrintE93( FILE*, double );
00966
00972 char *PrintEfmt(const char *fmt, double val );
00973
00975 const double SEXP_LIMIT = 84.;
00977 double sexp(double x);
00978
00983 double dsexp(double x);
00984
00989 double plankf(long int ip);
00990
00997 double qg32( double, double, double(*)(double) );
00998
00999
01000
01001
01011 void spsort( float x[], long int n, long int iperm[], int kflag, int *ier);
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021 #ifdef _MSC_VER
01022
01023 # pragma warning( disable : 4127 )
01024
01025 # pragma warning( disable : 4996 )
01026
01027 # pragma warning( disable : 4056 )
01028
01029 # pragma warning( disable : 4514 )
01030 #endif
01031 #ifdef __INTEL_COMPILER
01032 # pragma warning( disable : 1572 )
01033 #endif
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043