00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 #include <cstdarg>
00045 #include "cddefines.h"
00046 #include "physconst.h"
00047 #include "cddrive.h"
00048 #include "called.h"
00049 #include "opacity.h"
00050 #include "rfield.h"
00051 #include "hextra.h"
00052 #include "struc.h"
00053 #include "hmi.h"
00054 #include "fudgec.h"
00055 #include "broke.h"
00056 #include "trace.h"
00057 #include "input.h"
00058 #include "save.h"
00059 #include "version.h"
00060 #include "warnings.h"
00061 #include "conv.h"
00062 #include "thirdparty.h"
00063
00064
00065
00066 char *read_whole_line( char *chLine , int nChar , FILE *ioIN )
00067 {
00068 char *chRet;
00069
00070 DEBUG_ENTRY( "read_whole_line()" );
00071
00072
00073 memset( chLine, 0, (size_t)nChar );
00074
00075
00076
00077 if( (chRet = fgets( chLine, nChar, ioIN )) == NULL )
00078 return NULL;
00079
00080
00081 char *chEOL = (char*)memchr( chLine, '\n', nChar );
00082
00083
00084
00085
00086 if( chEOL == NULL )
00087 {
00088 if( called.lgTalk )
00089 fprintf( ioQQQ, "DISASTER PROBLEM read_whole_line found input"
00090 " with a line too long to be read, limit is %i char. "
00091 "Start of line follows:\n%s\n",
00092 nChar , chLine );
00093
00094 lgAbort = true;
00095 return NULL;
00096 }
00097 return chRet;
00098 }
00099
00101 void Split(const string& str,
00102 const string& sep,
00103 vector<string>& lst,
00104 split_mode mode)
00105 {
00106 DEBUG_ENTRY( "Split()" );
00107
00108 bool lgStrict = ( mode == SPM_STRICT );
00109 bool lgKeep = ( mode == SPM_KEEP_EMPTY );
00110 bool lgFail = false;
00111 string::size_type ptr1 = 0;
00112 string::size_type ptr2 = str.find( sep );
00113 string sstr = str.substr( ptr1, ptr2-ptr1 );
00114 if( sstr.length() > 0 )
00115 lst.push_back( sstr );
00116 else {
00117 if( lgStrict ) lgFail = true;
00118 if( lgKeep ) lst.push_back( sstr );
00119 }
00120 while( ptr2 != string::npos ) {
00121
00122 ptr1 = ptr2 + sep.length();
00123 if( ptr1 < str.length() ) {
00124 ptr2 = str.find( sep, ptr1 );
00125 sstr = str.substr( ptr1, ptr2-ptr1 );
00126 if( sstr.length() > 0 )
00127 lst.push_back( sstr );
00128 else {
00129 if( lgStrict ) lgFail = true;
00130 if( lgKeep ) lst.push_back( sstr );
00131 }
00132 }
00133 else {
00134 ptr2 = string::npos;
00135 if( lgStrict ) lgFail = true;
00136 if( lgKeep ) lst.push_back( "" );
00137 }
00138 }
00139 if( lgFail )
00140 {
00141 fprintf( ioQQQ, " A syntax error occurred while splitting the string: \"%s\"\n", str.c_str() );
00142 fprintf( ioQQQ, " The separator is \"%s\". Empty substrings are not allowed.\n", sep.c_str() );
00143 cdEXIT(EXIT_FAILURE);
00144 }
00145 }
00146
00147
00148 void MyAssert(const char *file, int line, const char *comment)
00149 {
00150 DEBUG_ENTRY( "MyAssert()" );
00151
00152 fprintf(ioQQQ,"\n\n\n PROBLEM DISASTER\n An assert has been thrown, this is bad.\n");
00153 fprintf(ioQQQ," %s\n",comment);
00154 fprintf(ioQQQ," It happened in the file %s at line number %i\n", file, line );
00155 fprintf(ioQQQ," This is iteration %li, nzone %li, fzone %.2f, lgSearch=%c.\n",
00156 iteration ,
00157 nzone ,
00158 fnzone ,
00159 TorF(conv.lgSearch) );
00160
00161 ShowMe();
00162 # ifdef OLD_ASSERT
00163 cdEXIT(EXIT_FAILURE);
00164 # endif
00165 }
00166
00167
00168 double AnuUnit(realnum energy_ryd)
00169
00170 {
00171 DEBUG_ENTRY( "AnuUnit()" );
00172
00173 return Energy((double)energy_ryd).get(save.chConPunEnr[save.ipConPun]);
00174 }
00175
00176
00177 void ShowMe(void)
00178 {
00179
00180 DEBUG_ENTRY( "ShowMe()" );
00181
00182
00183 if( ioQQQ != NULL )
00184 {
00185
00186 if( (hextra.cryden == 0.) && hmi.BiggestH2 > 0.1 )
00187 {
00188 fprintf( ioQQQ, " >>> \n >>> \n >>> Cosmic rays are not included and the gas is molecular. "
00189 "THIS IS KNOWN TO BE UNSTABLE. Add cosmic rays and try again.\n >>> \n >>>\n\n");
00190 }
00191 else
00192 {
00193 fprintf( ioQQQ, "\n\n\n" );
00194 fprintf( ioQQQ, " vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv \n" );
00195 fprintf( ioQQQ, " > PROBLEM DISASTER PROBLEM DISASTER. <\n" );
00196 fprintf( ioQQQ, " > Sorry, something bad has happened. <\n" );
00197 fprintf( ioQQQ, " > Please post this on the Cloudy web site <\n" );
00198 fprintf( ioQQQ, " > discussion board at www.nublado.org <\n" );
00199 fprintf( ioQQQ, " > Please send all following information: <\n" );
00200 fprintf( ioQQQ, " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ \n" );
00201 fprintf( ioQQQ, "\n\n" );
00202
00203
00204 fprintf( ioQQQ, " Cloudy version number is %s\n",
00205 t_version::Inst().chVersion );
00206 fprintf( ioQQQ, " %s\n\n", t_version::Inst().chInfo );
00207
00208 fprintf( ioQQQ, "%5ld warnings,%3ld cautions,%3ld temperature failures. Messages follow.\n",
00209 warnings.nwarn, warnings.ncaun, conv.nTeFail );
00210
00211
00212 cdWarnings(ioQQQ);
00213
00214
00215 cdCautions(ioQQQ);
00216
00217
00218 cdPrintCommands(ioQQQ);
00219
00220
00221
00222 if( input.nSaveIni )
00223 {
00224 fprintf(ioQQQ," This input stream included an init file.\n");
00225 fprintf(ioQQQ," If this init file is not part of the standard Cloudy distribution\n");
00226 fprintf(ioQQQ," then I will need a copy of it too.\n");
00227 }
00228 }
00229 }
00230 return;
00231 }
00232
00233
00234 void cap4(
00235 char *chCAP ,
00236
00237 const char *chLab)
00238 {
00239 long int
00240 i;
00241
00242 DEBUG_ENTRY( "cap4()" );
00243
00244
00245 for( i=0; i < 4; i++ )
00246 {
00247
00248 chCAP[i] = toupper( chLab[i] );
00249 }
00250
00251
00252 chCAP[4] = '\0';
00253 return;
00254 }
00255
00256
00257 void uncaps(char *chCard )
00258 {
00259 long int i;
00260
00261 DEBUG_ENTRY( "caps()" );
00262
00263
00264 i = 0;
00265 while( chCard[i]!= '\0' )
00266 {
00267 chCard[i] = tolower( chCard[i] );
00268 ++i;
00269 }
00270 return;
00271 }
00272
00273
00274 void caps(char *chCard )
00275 {
00276 long int i;
00277
00278 DEBUG_ENTRY( "caps()" );
00279
00280
00281 i = 0;
00282 while( chCard[i]!= '\0' )
00283 {
00284 chCard[i] = toupper( chCard[i] );
00285 ++i;
00286 }
00287 return;
00288 }
00289
00290
00291
00292
00293 double e2(
00294
00295 double t )
00296 {
00297
00298
00299 double hold = sexp(t) - t*ee1(t);
00300 DEBUG_ENTRY( "e2()" );
00301
00302 return max( hold, 0. );
00303 }
00304
00305
00306 double ee1(double x)
00307 {
00308 double ans,
00309 bot,
00310 top;
00311 static double a[6]={-.57721566,.99999193,-.24991055,.05519968,-.00976004,
00312 .00107857};
00313 static double b[4]={8.5733287401,18.0590169730,8.6347608925,.2677737343};
00314 static double c[4]={9.5733223454,25.6329561486,21.0996530827,3.9584969228};
00315
00316 DEBUG_ENTRY( "ee1()" );
00317
00318
00319
00320
00321
00322
00323
00324
00325 if( x <= 0 )
00326 {
00327 fprintf( ioQQQ, " DISASTER negative argument in function ee1, x<0\n" );
00328 cdEXIT(EXIT_FAILURE);
00329 }
00330
00331
00332 else if( x < 1. )
00333 {
00334
00335 ans = ((((a[5]*x + a[4])*x + a[3])*x + a[2])*x + a[1])*x + a[0] - log(x);
00336 }
00337
00338
00339 else
00340 {
00341
00342 top = (((x + b[0])*x + b[1])*x + b[2])*x + b[3];
00343 bot = (((x + c[0])*x + c[1])*x + c[2])*x + c[3];
00344 ans = top/bot/x*exp(-x);
00345 }
00346 return ans;
00347 }
00348
00349
00350 double ee1_safe(double x)
00351 {
00352 double ans,
00353 bot,
00354 top;
00355
00356
00357 static double b[4]={8.5733287401,18.0590169730,8.6347608925,.2677737343};
00358 static double c[4]={9.5733223454,25.6329561486,21.0996530827,3.9584969228};
00359
00360 DEBUG_ENTRY( "ee1_safe()" );
00361
00362 ASSERT( x > 1. );
00363
00364
00365
00366 top = (((x + b[0])*x + b[1])*x + b[2])*x + b[3];
00367
00368 bot = (((x + c[0])*x + c[1])*x + c[2])*x + c[3];
00369
00370 ans = top/bot/x;
00371 return ans;
00372 }
00373
00374
00375 double FFmtRead(const char *chCard,
00376 long int *ipnt,
00377
00378 long int last,
00379 bool *lgEOL)
00380 {
00381 DEBUG_ENTRY( "FFmtRead()" );
00382
00383 char chr = '\0';
00384 const char *eol_ptr = &chCard[last];
00385 const char *ptr = min(&chCard[*ipnt-1],eol_ptr);
00386
00387 ASSERT( *ipnt > 0 && last - *ipnt < INPUT_LINE_LENGTH );
00388
00389 while( ptr < eol_ptr && ( chr = *ptr++ ) != '\0' )
00390 {
00391 const char *lptr = ptr;
00392 char lchr = chr;
00393 if( lchr == '-' || lchr == '+' )
00394 lchr = *lptr++;
00395 if( lchr == '.' )
00396 lchr = *lptr;
00397 if( isdigit(lchr) )
00398 break;
00399 }
00400
00401 if( ptr == eol_ptr || chr == '\0' )
00402 {
00403 *ipnt = last+1;
00404 *lgEOL = true;
00405 return 0.;
00406 }
00407
00408 string chNumber;
00409 bool lgCommaFound = false, lgLastComma = false;
00410 do
00411 {
00412 lgCommaFound = lgLastComma;
00413 if( chr != ',' )
00414 {
00415 chNumber += chr;
00416 }
00417 else
00418 {
00419
00420
00421 lgLastComma = true;
00422
00423 }
00424 if( ptr == eol_ptr )
00425 break;
00426 chr = *ptr++;
00427 }
00428 while( isdigit(chr) || chr == '.' || chr == '-' || chr == '+' || chr == ',' || chr == 'e' || chr == 'E' );
00429
00430 if( lgCommaFound )
00431 {
00432 fprintf( ioQQQ, " PROBLEM - a comma was found embedded in a number, this is deprecated.\n" );
00433 fprintf(ioQQQ, "== %-80s ==\n",chCard);
00434 }
00435
00436 double value = atof(chNumber.c_str());
00437
00438 *ipnt = (long)(ptr - chCard);
00439 *lgEOL = false;
00440 return value;
00441 }
00442
00443
00444
00445 long nMatch(const char *chKey,
00446 const char *chCard)
00447 {
00448 const char *ptr;
00449 long Match_v;
00450
00451 DEBUG_ENTRY( "nMatch()" );
00452
00453 ASSERT( strlen(chKey) > 0 );
00454
00455 if( ( ptr = strstr_s( chCard, chKey ) ) == NULL )
00456 {
00457
00458 Match_v = 0L;
00459 }
00460 else
00461 {
00462
00463 Match_v = (long)(ptr-chCard+1);
00464 }
00465 return Match_v;
00466 }
00467
00468
00469
00470
00471
00472
00473
00474
00475 double fudge(long int ipnt)
00476 {
00477 double fudge_v;
00478
00479 DEBUG_ENTRY( "fudge()" );
00480
00481 if( ipnt < 0 )
00482 {
00483
00484 fudge_v = fudgec.nfudge;
00485 fudgec.lgFudgeUsed = true;
00486 }
00487 else if( ipnt >= fudgec.nfudge )
00488 {
00489 fprintf( ioQQQ, " FUDGE factor not entered for array number %3ld\n",
00490 ipnt );
00491 cdEXIT(EXIT_FAILURE);
00492 }
00493 else
00494 {
00495 fudge_v = fudgec.fudgea[ipnt];
00496 fudgec.lgFudgeUsed = true;
00497 }
00498 return fudge_v;
00499 }
00500
00501
00502
00503
00504
00505
00506
00507 int GetQuote(
00508
00509 char *chStringOut,
00510
00511 char *chCard,
00512 char *chCardRaw,
00513
00514
00515 bool lgAbort )
00516 {
00517 char *i0,
00518 *i1,
00519 *iLoc;
00520 size_t len;
00521
00522 DEBUG_ENTRY( "GetQuote()" );
00523
00524
00525 i0 = strchr_s( chCardRaw,'\"' );
00526
00527 if( i0 != NULL )
00528 {
00529
00530 i1 = strchr_s( i0+1 ,'\"' );
00531 }
00532 else
00533 {
00534 i1 = NULL;
00535 }
00536
00537
00538
00539 if( i0 == NULL || i1 == NULL )
00540 {
00541 if( lgAbort )
00542 {
00543
00544 fprintf( ioQQQ,
00545 " A filename or label must be specified within double quotes, but no quotes were encountered on this command.\n" );
00546 fprintf( ioQQQ, " Name must be surrounded by exactly two double quotes, like \"name.txt\". \n" );
00547 fprintf( ioQQQ, " The line image follows:\n" );
00548 fprintf( ioQQQ, " %s\n", chCardRaw);
00549 fprintf( ioQQQ, " Sorry\n" );
00550 cdEXIT(EXIT_FAILURE);
00551 }
00552 else
00553 {
00554
00555 chStringOut[0] = '\0';
00556
00557 return 1;
00558 }
00559 }
00560
00561
00562 len = (size_t)(i1-i0-1);
00563 strncpy(chStringOut,i0+1,len);
00564
00565 chStringOut[len] = '\0';
00566
00567
00568 iLoc = strchr_s( chCard, '\"' );
00569 if( iLoc == NULL )
00570 TotalInsanity();
00571
00572
00573
00574 while( i0 <= i1 )
00575 {
00576 *i0 = ' ';
00577 *iLoc = ' ';
00578 ++i0;
00579 ++iLoc;
00580 }
00581
00582 return 0;
00583 }
00584
00585
00586 #ifndef HAVE_POWI
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598 double powi( double x, long int n )
00599
00600
00601 {
00602 double p;
00603
00604 DEBUG_ENTRY( "powi()" );
00605
00606 if( x == 0 )
00607 return 0.;
00608
00609
00610 if( n < 0 )
00611 {
00612 n = -n;
00613 x = 1/x;
00614 }
00615
00616 p = is_odd(n) ? x : 1;
00617
00618
00619
00620 while( n >>= 1 )
00621 {
00622 x *= x;
00623 if( is_odd(n) )
00624 p *= x;
00625 }
00626
00627
00628 return p;
00629 }
00630
00631 #endif
00632
00633 long ipow( long m, long n )
00634
00635
00636 {
00637 long p;
00638
00639 DEBUG_ENTRY( "ipow()" );
00640
00641 if( m == 0 || (n < 0 && m > 1) )
00642 return 0L;
00643
00644
00645
00646 if( n < 0 )
00647 {
00648 n = -n;
00649 m = 1/m;
00650 }
00651
00652 p = is_odd(n) ? m : 1;
00653
00654
00655
00656 while( n >>= 1 )
00657 {
00658 m *= m;
00659 if( is_odd(n) )
00660 p *= m;
00661 }
00662
00663
00664 return p;
00665 }
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677 #ifdef _MSC_VER
00678
00679
00680
00681
00682
00683
00684 char *PrintEfmt(const char *fmt, double val )
00685 {
00686 static char buf[30];
00687
00688 DEBUG_ENTRY( "PrintEfmt()" );
00689
00690
00691 sprintf(buf, fmt, val);
00692
00693
00694
00695
00696 char *ep , buf2[30];
00697
00698
00699
00700 if( val >= 0.)
00701 {
00702 strcpy(buf2 , " " );
00703 strcat(buf2 , buf);
00704 strcpy( buf , buf2);
00705 }
00706
00707
00708 if((ep = strchr_s(buf, 'e')) == NULL)
00709 {
00710 ep = strchr_s(buf, 'E');
00711 }
00712
00713
00714 if(ep != NULL)
00715 {
00716
00717 ep += 2;
00718
00719
00720 *ep = '\0';
00721
00722
00723 ++ep;
00724
00725
00726 strcat( buf, ep );
00727 }
00728 return buf;
00729 }
00730 #endif
00731
00732
00733 void PrintE82( FILE* ioOUT, double value )
00734 {
00735 double frac , xlog , xfloor , tvalue;
00736 int iExp;
00737
00738 DEBUG_ENTRY( "PrintE82()" );
00739
00740 if( value < 0 )
00741 {
00742 fprintf(ioOUT,"********");
00743 }
00744 else if( value == 0 )
00745 {
00746 fprintf(ioOUT,"0.00E+00");
00747 }
00748 else
00749 {
00750
00751 tvalue = value;
00752 xlog = log10( tvalue );
00753 xfloor = floor(xlog);
00754
00755 frac = tvalue*pow(10.,-xfloor);
00756
00757 iExp = (int)xfloor;
00758 if( frac>9.9945 )
00759 {
00760 frac /= 10.;
00761 iExp += 1;
00762 }
00763
00764 fprintf(ioOUT,"%.2f",frac);
00765
00766 fprintf(ioOUT,"E");
00767
00768 if(iExp>=0 )
00769 {
00770 fprintf(ioOUT,"+");
00771 }
00772 fprintf(ioOUT,"%.2d",iExp);
00773 }
00774 return;
00775 }
00776
00777
00778
00779 void PrintE71( FILE* ioOUT, double value )
00780 {
00781 double frac , xlog , xfloor , tvalue;
00782 int iExp;
00783
00784 DEBUG_ENTRY( "PrintE71()" );
00785
00786 if( value < 0 )
00787 {
00788 fprintf(ioOUT,"*******");
00789 }
00790 else if( value == 0 )
00791 {
00792 fprintf(ioOUT,"0.0E+00");
00793 }
00794 else
00795 {
00796
00797 tvalue = value;
00798 xlog = log10( tvalue );
00799 xfloor = floor(xlog);
00800
00801 frac = tvalue*pow(10.,-xfloor);
00802
00803 iExp = (int)xfloor;
00804 if( frac>9.9945 )
00805 {
00806 frac /= 10.;
00807 iExp += 1;
00808 }
00809
00810 fprintf(ioOUT,"%.1f",frac);
00811
00812 fprintf(ioOUT,"E");
00813
00814 if(iExp>=0 )
00815 {
00816 fprintf(ioOUT,"+");
00817 }
00818 fprintf(ioOUT,"%.2d",iExp);
00819 }
00820 return;
00821 }
00822
00823
00824
00825
00826 void PrintE93( FILE* ioOUT, double value )
00827 {
00828 double frac , xlog , xfloor, tvalue;
00829 int iExp;
00830
00831 DEBUG_ENTRY( "PrintE93()" );
00832
00833 if( value < 0 )
00834 {
00835 fprintf(ioOUT,"*********");
00836 }
00837 else if( value == 0 )
00838 {
00839 fprintf(ioOUT,"0.000E+00");
00840 }
00841 else
00842 {
00843
00844 tvalue = value;
00845 xlog = log10( tvalue );
00846 xfloor = floor(xlog);
00847
00848 frac = tvalue*pow(10.,-xfloor);
00849
00850 iExp = (int)xfloor;
00851 if( frac>9.99949 )
00852 {
00853 frac /= 10.;
00854 iExp += 1;
00855 }
00856
00857 fprintf(ioOUT,"%5.3f",frac);
00858
00859 fprintf(ioOUT,"E");
00860
00861 if(iExp>=0 )
00862 {
00863 fprintf(ioOUT,"+");
00864 }
00865 fprintf(ioOUT,"%.2d",iExp);
00866 }
00867 return;
00868 }
00869
00870
00871 NORETURN void TotalInsanity(void)
00872 {
00873 DEBUG_ENTRY( "TotalInsanity()" );
00874
00875
00876
00877
00878 fprintf( ioQQQ, " Something that cannot happen, has happened.\n" );
00879 fprintf( ioQQQ, " This is TotalInsanity, I live in %s.\n", __FILE__ );
00880 ShowMe();
00881
00882 cdEXIT(EXIT_FAILURE);
00883 }
00884
00885
00886 NORETURN void BadRead(void)
00887 {
00888 DEBUG_ENTRY( "BadRead()" );
00889
00890
00891 fprintf( ioQQQ, " A read of internal input data has failed.\n" );
00892 fprintf( ioQQQ, " This is BadRead, I live in %s.\n", __FILE__ );
00893 ShowMe();
00894
00895 cdEXIT(EXIT_FAILURE);
00896 }
00897
00898
00899 sys_float sexp(sys_float x)
00900 {
00901 sys_float sexp_v;
00902
00903 DEBUG_ENTRY( "sexp()" );
00904
00905
00906 if( x < SEXP_LIMIT )
00907 {
00908 sexp_v = exp(-x);
00909 }
00910 else
00911 {
00912 sexp_v = 0.f;
00913 }
00914 return sexp_v;
00915 }
00916
00917
00918 double sexp(double x)
00919 {
00920 double sexp_v;
00921
00922 DEBUG_ENTRY( "sexp()" );
00923
00924
00925 if( x < SEXP_LIMIT )
00926 {
00927 sexp_v = exp(-x);
00928 }
00929 else
00930 {
00931 sexp_v = 0.;
00932 }
00933 return sexp_v;
00934 }
00935
00936
00937
00938 double dsexp(double x)
00939 {
00940 double dsexp_v;
00941
00942 DEBUG_ENTRY( "dsexp()" );
00943
00944 if( x < DSEXP_LIMIT )
00945 {
00946 dsexp_v = exp(-x);
00947 }
00948 else
00949 {
00950 dsexp_v = 0.;
00951 }
00952 return dsexp_v;
00953 }
00954
00955
00956
00957 void TestCode(void)
00958 {
00959 DEBUG_ENTRY( "TestCode( )" );
00960
00961
00962 lgTestCodeCalled = true;
00963 return;
00964 }
00965
00966
00967 void broken(void)
00968 {
00969 DEBUG_ENTRY( "broken( )" );
00970
00971 broke.lgBroke = true;
00972 return;
00973 }
00974
00975
00976 void fixit(void)
00977 {
00978 DEBUG_ENTRY( "fixit( )" );
00979
00980 broke.lgFixit = true;
00981 return;
00982 }
00983
00984
00985 void CodeReview(void)
00986 {
00987 DEBUG_ENTRY( "CodeReview( )" );
00988
00989 broke.lgCheckit = true;
00990 return;
00991 }
00992
00994 int dprintf(FILE *fp, const char *format, ...)
00995 {
00996 va_list ap;
00997 int i1, i2;
00998
00999 DEBUG_ENTRY( "dprintf()" );
01000 va_start(ap,format);
01001 i1 = fprintf(fp,"DEBUG ");
01002 if (i1 >= 0)
01003 i2 = vfprintf(fp,format,ap);
01004 else
01005 i2 = 0;
01006 if (i2 < 0)
01007 i1 = 0;
01008 va_end(ap);
01009
01010 return i1+i2;
01011 }
01012
01013
01014
01015
01016 int dbg_printf(int debug, const char *fmt, ...)
01017 {
01018 va_list ap;
01019 int i=0;
01020
01021 DEBUG_ENTRY( "dbg_printf()" );
01022
01023
01024 if(debug <= trace.debug_level)
01025 {
01026 va_start(ap, fmt);
01027
01028 i = vfprintf(ioQQQ, fmt, ap);
01029
01030 fflush(ioQQQ);
01031 va_end(ap);
01032 }
01033 return i;
01034 }
01035
01036
01037
01038 double qg32(
01039 double xl,
01040 double xu,
01041
01042 double (*fct)(double) )
01043 {
01044 double a = 0.5*(xu+xl),
01045 b = xu-xl,
01046 y = 0.;
01047
01048 DEBUG_ENTRY( "qg32()" );
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063 double weights[16] = {
01064 .35093050047350483e-2, .81371973654528350e-2, .12696032654631030e-1, .17136931456510717e-1,
01065 .21417949011113340e-1, .25499029631188088e-1, .29342046739267774e-1, .32911111388180923e-1,
01066 .36172897054424253e-1, .39096947893535153e-1, .41655962113473378e-1, .43826046502201906e-1,
01067 .45586939347881942e-1, .46922199540402283e-1, .47819360039637430e-1, .48270044257363900e-1};
01068
01069 double c[16] = {
01070 .498631930924740780, .49280575577263417, .4823811277937532200, .46745303796886984000,
01071 .448160577883026060, .42468380686628499, .3972418979839712000, .36609105937014484000,
01072 .331522133465107600, .29385787862038116, .2534499544661147000, .21067563806531767000,
01073 .165934301141063820, .11964368112606854, .7223598079139825e-1, .24153832843869158e-1};
01074
01075 for( int i=0; i<16; i++)
01076 {
01077 y += b * weights[i] * ((*fct)(a+b*c[i]) + (*fct)(a-b*c[i]));
01078 }
01079
01080
01081 return y;
01082 }
01083
01084
01085 void spsort(
01086
01087 realnum x[],
01088
01089 long int n,
01090
01091 long int iperm[],
01092
01093
01094 int kflag,
01095
01096 int *ier)
01097 {
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169 long int i,
01170 ij,
01171 il[21],
01172 indx,
01173 indx0,
01174 istrt,
01175 istrt_,
01176 iu[21],
01177 j,
01178 k,
01179 kk,
01180 l,
01181 lm,
01182 lmt,
01183 m,
01184 nn;
01185 realnum r,
01186 ttemp;
01187
01188 DEBUG_ENTRY( "spsort()" );
01189
01190
01191
01192
01193
01194
01195
01196 *ier = 0;
01197 nn = n;
01198 if( nn < 1 )
01199 {
01200 *ier = 1;
01201 return;
01202 }
01203 else
01204 {
01205 kk = labs(kflag);
01206 if( kk != 1 && kk != 2 )
01207 {
01208 *ier = 2;
01209 return;
01210 }
01211 else
01212 {
01213
01214
01215
01216 for( i=0; i < nn; i++ )
01217 {
01218 iperm[i] = i+1;
01219 }
01220
01221
01222 if( nn == 1 )
01223 {
01224 --iperm[0];
01225 return;
01226 }
01227
01228
01229 if( kflag <= -1 )
01230 {
01231 for( i=0; i < nn; i++ )
01232 {
01233 x[i] = -x[i];
01234 }
01235 }
01236
01237
01238 m = 1;
01239 i = 1;
01240 j = nn;
01241 r = .375e0;
01242 }
01243 }
01244
01245 while( true )
01246 {
01247 if( i == j )
01248 goto L_80;
01249 if( r <= 0.5898437e0 )
01250 {
01251 r += 3.90625e-2;
01252 }
01253 else
01254 {
01255 r -= 0.21875e0;
01256 }
01257
01258 L_40:
01259 k = i;
01260
01261
01262
01263 ij = i + (long)((j-i)*r);
01264 lm = iperm[ij-1];
01265
01266
01267
01268 if( x[iperm[i-1]-1] > x[lm-1] )
01269 {
01270 iperm[ij-1] = iperm[i-1];
01271 iperm[i-1] = lm;
01272 lm = iperm[ij-1];
01273 }
01274 l = j;
01275
01276
01277
01278 if( x[iperm[j-1]-1] < x[lm-1] )
01279 {
01280 iperm[ij-1] = iperm[j-1];
01281 iperm[j-1] = lm;
01282 lm = iperm[ij-1];
01283
01284
01285
01286
01287 if( x[iperm[i-1]-1] > x[lm-1] )
01288 {
01289 iperm[ij-1] = iperm[i-1];
01290 iperm[i-1] = lm;
01291 lm = iperm[ij-1];
01292 }
01293 }
01294
01295
01296
01297 while( true )
01298 {
01299 l -= 1;
01300 if( x[iperm[l-1]-1] <= x[lm-1] )
01301 {
01302
01303
01304
01305 while( true )
01306 {
01307 k += 1;
01308 if( x[iperm[k-1]-1] >= x[lm-1] )
01309 break;
01310 }
01311
01312
01313 if( k > l )
01314 break;
01315 lmt = iperm[l-1];
01316 iperm[l-1] = iperm[k-1];
01317 iperm[k-1] = lmt;
01318 }
01319 }
01320
01321
01322 if( l - i > j - k )
01323 {
01324 il[m-1] = i;
01325 iu[m-1] = l;
01326 i = k;
01327 m += 1;
01328 }
01329 else
01330 {
01331 il[m-1] = k;
01332 iu[m-1] = j;
01333 j = l;
01334 m += 1;
01335 }
01336
01337 L_90:
01338 if( j - i >= 1 )
01339 goto L_40;
01340 if( i == 1 )
01341 continue;
01342 i -= 1;
01343
01344 while( true )
01345 {
01346 i += 1;
01347 if( i == j )
01348 break;
01349 lm = iperm[i];
01350 if( x[iperm[i-1]-1] > x[lm-1] )
01351 {
01352 k = i;
01353
01354 while( true )
01355 {
01356 iperm[k] = iperm[k-1];
01357 k -= 1;
01358
01359 if( x[lm-1] >= x[iperm[k-1]-1] )
01360 break;
01361 }
01362 iperm[k] = lm;
01363 }
01364 }
01365
01366
01367 L_80:
01368 m -= 1;
01369 if( m == 0 )
01370 break;
01371
01372 i = il[m-1];
01373 j = iu[m-1];
01374
01375 goto L_90;
01376 }
01377
01378
01379 if( kflag <= -1 )
01380 {
01381 for( i=0; i < nn; i++ )
01382 {
01383 x[i] = -x[i];
01384 }
01385 }
01386
01387
01388 if( kk == 2 )
01389 {
01390
01391
01392
01393 for( istrt=1; istrt <= nn; istrt++ )
01394 {
01395 istrt_ = istrt - 1;
01396 if( iperm[istrt_] >= 0 )
01397 {
01398 indx = istrt;
01399 indx0 = indx;
01400 ttemp = x[istrt_];
01401 while( iperm[indx-1] > 0 )
01402 {
01403 x[indx-1] = x[iperm[indx-1]-1];
01404 indx0 = indx;
01405 iperm[indx-1] = -iperm[indx-1];
01406 indx = labs(iperm[indx-1]);
01407 }
01408 x[indx0-1] = ttemp;
01409 }
01410 }
01411
01412
01413 for( i=0; i < nn; i++ )
01414 {
01415 iperm[i] = -iperm[i];
01416 }
01417 }
01418
01419 for( i=0; i < nn; i++ )
01420 {
01421 --iperm[i];
01422 }
01423 return;
01424 }
01425
01426
01427
01428
01429
01430
01431
01432
01433 void *MyMalloc(
01434
01435 size_t size ,
01436 const char *chFile,
01437 int line
01438 )
01439 {
01440 void *ptr;
01441
01442 DEBUG_ENTRY( "MyMalloc()" );
01443
01444 ASSERT( size > 0 );
01445
01446
01447 {
01448 enum{DEBUG_LOC=false};
01449 if( DEBUG_LOC)
01450 {
01451 static long int kount=0, nTot=0;
01452 nTot += (long)size;
01453 fprintf(ioQQQ,"%li\t%li\t%li\n",
01454 kount ,
01455 (long)size ,
01456 nTot );
01457 ++kount;
01458 }
01459 }
01460
01461 if( ( ptr = malloc( size ) ) == NULL )
01462 {
01463 fprintf(ioQQQ,"DISASTER MyMalloc could not allocate %lu bytes. Exit in MyMalloc.",
01464 (unsigned long)size );
01465 fprintf(ioQQQ,"MyMalloc called from file %s at line %i.\n",
01466 chFile , line );
01467
01468 if( struc.nzlim>2000 )
01469 fprintf(ioQQQ,"This may have been caused by the large number of zones."
01470 " %li zones were requested. Is this many zones really necessary?\n",
01471 struc.nzlim );
01472
01473 cdEXIT(EXIT_FAILURE);
01474 }
01475
01476
01477 # if !defined(NDEBUG) && !defined(NOINIT)
01478
01479 size_t nFloat = size/4;
01480 size_t nDouble = size/8;
01481 sys_float *fptr = static_cast<sys_float*>(ptr);
01482 double *dptr = static_cast<double*>(ptr);
01483
01484
01485
01486
01487
01488 if( size == nDouble*8 )
01489 {
01490
01491
01492
01493
01494
01495
01496
01497 set_NaN( dptr, (long)nDouble );
01498 }
01499 else if( size == nFloat*4 )
01500 {
01501
01502 set_NaN( fptr, (long)nFloat );
01503 }
01504 else
01505 {
01506 memset( ptr, 0xff, size );
01507 }
01508
01509 # endif
01510 return ptr;
01511 }
01512
01513
01514
01515
01516
01517
01518 void *MyCalloc(
01519
01520 size_t num ,
01521 size_t size )
01522 {
01523 void *ptr;
01524
01525 DEBUG_ENTRY( "MyCalloc()" );
01526
01527 ASSERT( size > 0 );
01528
01529
01530 {
01531 enum{DEBUG_LOC=false};
01532 if( DEBUG_LOC)
01533 {
01534 static long int kount=0;
01535 fprintf(ioQQQ,"%li\tcall\t%li\tbytes\n", kount,
01536 (long)size );
01537 ++kount;
01538 }
01539 }
01540
01541 if( ( ptr = calloc( num , size ) ) == NULL )
01542 {
01543 fprintf(ioQQQ,"MyCalloc could not allocate %lu bytes. Exit in MyCalloc.",
01544 (unsigned long)size );
01545 cdEXIT(EXIT_FAILURE);
01546 }
01547 return ptr;
01548 }
01549
01550
01551
01552
01553
01554 void *MyRealloc(
01555
01556 void *p ,
01557 size_t size )
01558 {
01559 void *ptr;
01560
01561 DEBUG_ENTRY( "MyRealloc()" );
01562
01563 ASSERT( size > 0 );
01564
01565
01566 {
01567 enum{DEBUG_LOC=false};
01568 if( DEBUG_LOC)
01569 {
01570 static long int kount=0;
01571 fprintf(ioQQQ,"%li\tcall\t%li\tbytes\n", kount,
01572 (long)size );
01573 ++kount;
01574 }
01575 }
01576
01577 if( ( ptr = realloc( p , size ) ) == NULL )
01578 {
01579 fprintf(ioQQQ,"MyRealloc could not allocate %lu bytes. Exit in MyRealloc.",
01580 (unsigned long)size );
01581 cdEXIT(EXIT_FAILURE);
01582 }
01583 return ptr;
01584 }
01585
01586
01587 double csphot(
01588
01589
01590 long int inu,
01591
01592 long int ithr,
01593
01594 long int iofset)
01595 {
01596 double csphot_v;
01597
01598 DEBUG_ENTRY( "csphot()" );
01599
01600 csphot_v = opac.OpacStack[inu-ithr+iofset-1];
01601 return csphot_v;
01602 }
01603
01604
01605
01606
01607
01608
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627
01628 double RandGauss(
01629
01630 double xMean,
01631
01632 double s )
01633 {
01634 double x1, x2, w, yy1;
01635 static double yy2=-BIGDOUBLE;
01636 static int use_last = false;
01637
01638 DEBUG_ENTRY( "RandGauss()" );
01639
01640 if( use_last )
01641 {
01642 yy1 = yy2;
01643 use_last = false;
01644 }
01645 else
01646 {
01647 do {
01648 x1 = 2.*genrand_real3() - 1.;
01649 x2 = 2.*genrand_real3() - 1.;
01650 w = x1 * x1 + x2 * x2;
01651 } while ( w >= 1.0 );
01652
01653 w = sqrt((-2.0*log(w))/w);
01654 yy1 = x1 * w;
01655 yy2 = x2 * w;
01656 use_last = true;
01657 }
01658 return xMean + yy1 * s;
01659 }
01660
01661
01662
01663
01664
01665
01666
01667
01668 double MyGaussRand( double PctUncertainty )
01669 {
01670 double StdDev;
01671 double result;
01672
01673 DEBUG_ENTRY( "MyGaussRand()" );
01674
01675 ASSERT( PctUncertainty < 0.5 );
01676
01677 StdDev = PctUncertainty;
01678
01679 do
01680 {
01681
01682 result = 1.+RandGauss( 0., StdDev );
01683 }
01684
01685 while( (result < 1.-3.*PctUncertainty) || (result > 1.+3.*PctUncertainty) );
01686
01687 ASSERT( result>0. && result<2. );
01688 return result;
01689 }
01690
01691
01692 double plankf(long int ip)
01693 {
01694 double plankf_v;
01695
01696 DEBUG_ENTRY( "plankf()" );
01697
01698
01699
01700
01701 if( rfield.ContBoltz[ip] <= 0. )
01702 {
01703 plankf_v = 1e-35;
01704 }
01705 else
01706 {
01707 plankf_v = 6.991e-21*POW2(FR1RYD*rfield.anu[ip])/
01708 (1./rfield.ContBoltz[ip] - 1.)*FR1RYD*4.;
01709 }
01710 return plankf_v;
01711 }