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