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 <ctype.h>
00045 #include <stdarg.h>
00046
00047
00048
00049
00050
00051
00052
00053
00054 #include "cddefines.h"
00055 #include "physconst.h"
00056 #include "cddrive.h"
00057 #include "opacity.h"
00058 #include "rfield.h"
00059 #include "hextra.h"
00060 #include "hmi.h"
00061 #include "fudgec.h"
00062 #include "broke.h"
00063 #include "trace.h"
00064 #include "path.h"
00065 #include "input.h"
00066 #include "elementnames.h"
00067 #include "punch.h"
00068 #include "version.h"
00069 #include "warnings.h"
00070 #include "conv.h"
00071
00072
00073 NORETURN void MyAssert(const char *file, int line)
00074 {
00075 fprintf(ioQQQ," PROBLEM DISASTER An assert has been thrown, this is bad.\n");
00076 fprintf(ioQQQ," It happened in the file %s at line number %i\n", file, line );
00077 fprintf(ioQQQ," This is iteration %li, nzone %li, fzone %.2f, lgSearch=%c.\n",
00078 iteration ,
00079 nzone ,
00080 fnzone ,
00081 TorF(conv.lgSearch) );
00082
00083 ShowMe();
00084 if( lgAssertFPE )
00085 {
00086
00087
00088
00089 float a = 1.f/(float)ZeroNum;
00090 fprintf(ioQQQ," when I div by 0 I get %f\n", a);
00091 }
00092
00093 fflush( ioQQQ );
00094 puts( "[Stop in MyAssert]" );
00095 cdEXIT(EXIT_FAILURE);
00096 }
00097
00098
00099 void path_not_set( void )
00100 {
00101
00102 DEBUG_ENTRY( "AnuUnit()" );
00103
00104 fprintf(ioQQQ, "\n\nAlthough there may be other reasons you have received this error,\n");
00105 fprintf(ioQQQ, "the most likely reason is that the path has not been properly set.\n");
00106 fprintf(ioQQQ, "Please have a look at the file path.cpp in the source directory.\n");
00107 fprintf(ioQQQ, "Check how the variable chDataPath is set - it should give the location of the data files I need.\n");
00108 fprintf(ioQQQ, "These are the files in the data download from the nublado web site.\n");
00109 fprintf(ioQQQ, "I bet that the path has not been set properly - I could not find what I expected to find.\n");
00110 fprintf(ioQQQ, "The internal flags I use suggest that the path has");
00111 if( !lgDataPathSet )
00112 {
00113 fprintf(ioQQQ, " not");
00114 }
00115 fprintf(ioQQQ, " been set. Its value is %s\nSorry.\n\n\n", chDataPath);
00116
00117 DEBUG_EXIT( "AnuUnit()" );
00118
00119 return;
00120 }
00121
00122
00123 double AnuUnit(float energy_ryd)
00124
00125 {
00126 double AnuUnit_v;
00127
00128 DEBUG_ENTRY( "AnuUnit()" );
00129
00130
00131 if( energy_ryd <=0. )
00132 {
00133
00134 AnuUnit_v = 0.;
00135 }
00136 else if( strcmp(punch.chConPunEnr[punch.ipConPun],"ryd ") == 0 )
00137 {
00138
00139 AnuUnit_v = energy_ryd;
00140 }
00141 else if( strcmp(punch.chConPunEnr[punch.ipConPun],"micr") == 0 )
00142 {
00143
00144 AnuUnit_v = RYDLAM/energy_ryd*1e-4;
00145 }
00146 else if( strcmp(punch.chConPunEnr[punch.ipConPun]," kev") == 0 )
00147 {
00148
00149 AnuUnit_v = energy_ryd*EVRYD*1.e-3;
00150 }
00151 else if( strcmp(punch.chConPunEnr[punch.ipConPun]," ev ") == 0 )
00152 {
00153
00154 AnuUnit_v = energy_ryd*EVRYD;
00155 }
00156 else if( strcmp(punch.chConPunEnr[punch.ipConPun],"angs") == 0 )
00157 {
00158
00159 AnuUnit_v = RYDLAM/energy_ryd;
00160 }
00161 else if( strcmp(punch.chConPunEnr[punch.ipConPun],"cent") == 0 )
00162 {
00163
00164 AnuUnit_v = RYDLAM/energy_ryd*1e-8;
00165 }
00166 else if( strcmp(punch.chConPunEnr[punch.ipConPun],"wave") == 0 )
00167 {
00168
00169 AnuUnit_v = RYD_INF*energy_ryd;
00170 }
00171 else if( strcmp(punch.chConPunEnr[punch.ipConPun]," mhz") == 0 )
00172 {
00173
00174 AnuUnit_v = RYD_INF*energy_ryd*SPEEDLIGHT/1e6;
00175 }
00176 else
00177 {
00178 fprintf( ioQQQ, " insane units in AnuUnit =%4.4s\n",
00179 punch.chConPunEnr[punch.ipConPun] );
00180 puts( "[Stop in AnuUnit]" );
00181 cdEXIT(EXIT_FAILURE);
00182 }
00183
00184 DEBUG_EXIT( "AnuUnit()" );
00185
00186 return AnuUnit_v;
00187 }
00188
00189
00190 void ShowMe(void)
00191 {
00192
00193 DEBUG_ENTRY( "ShowMe()" );
00194
00195
00196 if( ioQQQ != NULL )
00197 {
00198
00199 if( (hextra.cryden == 0.) && hmi.BiggestH2 > 0.1 )
00200 {
00201 fprintf( ioQQQ, " >>> \n >>> \n >>> Cosmic rays are not included and the gas is molecular. "
00202 "THIS IS KNOWN TO BE UNSTABLE. Add cosmic rays and try again.\n >>> \n >>>\n\n");
00203 }
00204 else
00205 {
00206 fprintf( ioQQQ, "\n\n" );
00207 fprintf( ioQQQ, "\n\n" );
00208 fprintf( ioQQQ, " vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv \n" );
00209 fprintf( ioQQQ, " > PROBLEM DISASTER PROBLEM DISASTER. <\n" );
00210 fprintf( ioQQQ, " > Sorry, something bad has happened. <\n" );
00211 fprintf( ioQQQ, " > Please post this on the Cloudy web site <\n" );
00212 fprintf( ioQQQ, " > discussion board at www.nublado.org <\n" );
00213 fprintf( ioQQQ, " > Please send all following information: <\n" );
00214 fprintf( ioQQQ, " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ \n" );
00215 fprintf( ioQQQ, "\n\n" );
00216
00217
00218 fprintf( ioQQQ, " Cloudy version number is %s\n",
00219 version.chVersion );
00220 fprintf( ioQQQ, " %s\n\n", version.chInfo );
00221
00222 fprintf( ioQQQ, "%5ld warnings,%3ld cautions,%3ld temperature failures. Messages follow.\n",
00223 warnings.nwarn, warnings.ncaun, conv.nTeFail );
00224
00225
00226 cdWarnings(ioQQQ);
00227
00228
00229 cdCautions(ioQQQ);
00230
00231
00232 cdPrintCommands(ioQQQ);
00233
00234
00235
00236 if( input.nSaveIni )
00237 {
00238 fprintf(ioQQQ," This input stream included an init file.\n");
00239 fprintf(ioQQQ," If this init file is not part of the standard Cloudy distribution\n");
00240 fprintf(ioQQQ," then I will need a copy of it too.\n");
00241 }
00242 }
00243 }
00244
00245 DEBUG_EXIT( "ShowMe()" );
00246
00247 return;
00248 }
00249
00250
00251 void cap4(
00252 char *chCAP ,
00253
00254 char *chLab)
00255 {
00256 long int
00257 i;
00258
00259 DEBUG_ENTRY( "cap4()" );
00260
00261
00262 for( i=0; i < 4; i++ )
00263 {
00264
00265 chCAP[i] = (char)toupper( chLab[i] );
00266 }
00267
00268
00269 chCAP[4] = '\0';
00270
00271
00272 DEBUG_EXIT( "cap4()" );
00273
00274 return;
00275 }
00276
00277
00278 void caps(char *chCard )
00279 {
00280 long int i;
00281
00282 DEBUG_ENTRY( "caps()" );
00283
00284
00285 i = 0;
00286 while( chCard[i]!= '\0' )
00287 {
00288 chCard[i] = (char)toupper( chCard[i] );
00289 ++i;
00290 }
00291
00292 DEBUG_EXIT( "caps()" );
00293
00294 return;
00295 }
00296
00297
00298
00299
00300 double e2(
00301
00302 double t )
00303 {
00304
00305
00306 double hold = sexp(t) - t*ee1(t);
00307
00308 return max( hold, 0. );
00309 }
00310
00311
00312 double ee1(double x)
00313 {
00314 double ans,
00315 bot,
00316 top;
00317 static double a[6]={-.57721566,.99999193,-.24991055,.05519968,-.00976004,
00318 .00107857};
00319 static double b[4]={8.5733287401,18.0590169730,8.6347608925,.2677737343};
00320 static double c[4]={9.5733223454,25.6329561486,21.0996530827,3.9584969228};
00321
00322 DEBUG_ENTRY( "ee1()" );
00323
00324
00325
00326
00327
00328
00329
00330
00331 if( x <= 0 )
00332 {
00333 fprintf( ioQQQ, " DISASTER negative argument in function ee1, x<0\n" );
00334 puts( "[Stop in ee1]" );
00335 cdEXIT(EXIT_FAILURE);
00336 }
00337
00338
00339 else if( x < 1. )
00340 {
00341
00342 ans = ((((a[5]*x + a[4])*x + a[3])*x + a[2])*x + a[1])*x + a[0] - log(x);
00343 }
00344
00345
00346 else
00347 {
00348
00349 top = (((x + b[0])*x + b[1])*x + b[2])*x + b[3];
00350 bot = (((x + c[0])*x + c[1])*x + c[2])*x + c[3];
00351 ans = top/bot/x*exp(-x);
00352 }
00353
00354 DEBUG_EXIT( "ee1()" );
00355
00356 return ans;
00357 }
00358
00359
00360 double ee1_safe(double x)
00361 {
00362 double ans,
00363 bot,
00364 top;
00365
00366
00367 static double b[4]={8.5733287401,18.0590169730,8.6347608925,.2677737343};
00368 static double c[4]={9.5733223454,25.6329561486,21.0996530827,3.9584969228};
00369
00370 DEBUG_ENTRY( "ee1_safe()" );
00371
00372 ASSERT( x > 1. );
00373
00374
00375
00376 top = (((x + b[0])*x + b[1])*x + b[2])*x + b[3];
00377
00378 bot = (((x + c[0])*x + c[1])*x + c[2])*x + c[3];
00379
00380 ans = top/bot/x;
00381
00382 DEBUG_EXIT( "ee1_safe()" );
00383
00384 return ans;
00385 }
00386
00387
00388
00389 #if 1
00390
00391
00392
00393 static void nxtchr(char *chr,
00394 long int *ichr,
00395 long int ipnt,
00396 const char *chCard,
00397 long int last,
00398 bool *lgEOL)
00399 {
00400
00401 DEBUG_ENTRY( "nxtchr()" );
00402
00403
00404
00405
00406 *chr = chCard[ipnt - 1];
00407 *ichr = (*chr);
00408
00409
00410
00411
00412
00413
00414
00415 if( ipnt > last || *chr=='\0' )
00416 {
00417 *lgEOL = true;
00418 }
00419 else
00420 {
00421 *lgEOL = false;
00422 }
00423
00424 DEBUG_EXIT( "nxtchr()" );
00425
00426 return;
00427 }
00428
00429
00430
00431
00432
00433 double FFmtRead(const char *chCard,
00434 long int *ipnt,
00435
00436 long int last,
00437 bool *lgEOL)
00438 {
00439 bool lgNFound,
00440 lgNumber;
00441 char chr;
00442 long int ichr,
00443 l1,
00444 l2;
00445 double FFmtRead_v,
00446 expn,
00447 PM1,
00448 value;
00449
00450 DEBUG_ENTRY( "FFmtRead()" );
00451
00452
00453
00454
00455
00456 ASSERT(*ipnt > 0 );
00457
00458 l1 = '0';
00459 l2 = '9';
00460
00461 L_999:
00462 lgNumber = false;
00463 lgNFound = false;
00464 PM1 = 1.;
00465 value = 0.;
00466 FFmtRead_v = 0.;
00467
00468 nxtchr(&chr,&ichr,*ipnt,chCard,last,lgEOL);
00469 if( *lgEOL )
00470 {
00471 DEBUG_EXIT( "FFmtRead()" );
00472 return FFmtRead_v;
00473 }
00474
00475
00476
00477
00478 L_4:
00479 if( chr == '.' )
00480 goto L_5;
00481
00482
00483 if( (ichr >= l1) && (ichr <= l2) )
00484 goto L_1;
00485
00486 *ipnt += 1;
00487 nxtchr(&chr,&ichr,*ipnt,chCard,last,lgEOL);
00488 if( *lgEOL )
00489 {
00490 DEBUG_EXIT( "FFmtRead()" );
00491 return FFmtRead_v;
00492 }
00493 goto L_4;
00494
00495
00496
00497
00498 L_1:
00499
00500 if( *ipnt>1 )
00501 {
00502 if( chCard[*ipnt - 2] == '-' )
00503 PM1 = -1.;
00504 lgNFound = true;
00505 }
00506
00507
00508
00509
00510 L_6:
00511 value = 10.*value + (double)(labs(ichr-l1));
00512 lgNumber = true;
00513 L_7:
00514 *ipnt += 1;
00515 nxtchr(&chr,&ichr,*ipnt,chCard,last,lgEOL);
00516 if( *lgEOL )
00517 {
00518 FFmtRead_v = value*PM1;
00519
00520 *lgEOL = false;
00521
00522 DEBUG_EXIT( "FFmtRead()" );
00523 return FFmtRead_v;
00524 }
00525 if( chr == ',' )
00526 goto L_7;
00527
00528
00529 if( (ichr >= l1) && (ichr <= l2) )
00530 goto L_6;
00531
00532 if( *ipnt > 0 )
00533 {
00534 if( chCard[*ipnt - 1] != '.' )
00535 {
00536 FFmtRead_v = value*PM1;
00537
00538 DEBUG_EXIT( "FFmtRead()" );
00539 return FFmtRead_v;
00540 }
00541 }
00542
00543
00544 L_5:
00545 expn = 1.;
00546 if( *ipnt > 1 )
00547 {
00548 if( chCard[*ipnt - 2] == '-' )
00549 PM1 = -1.;
00550 }
00551
00552
00553 L_3:
00554 *ipnt += 1;
00555 nxtchr(&chr,&ichr,*ipnt,chCard,last,lgEOL);
00556 if( chr == ',' )
00557 goto L_3;
00558 if( (*lgEOL || ichr < l1) || ichr > l2 )
00559 {
00560
00561
00562
00563 if( !lgNFound && value==0. )
00564 goto L_999;
00565 FFmtRead_v = value*PM1;
00566 if( lgNumber )
00567 *lgEOL = false;
00568
00569 DEBUG_EXIT( "FFmtRead()" );
00570 return FFmtRead_v;
00571 }
00572 expn *= 0.1;
00573 lgNFound = true;
00574 lgNumber = true;
00575 value += expn*(double)(ichr-l1);
00576 goto L_3;
00577
00578 }
00579
00580 #else
00581 double FFmtRead(const char *chCard,
00582 long int *ipnt,
00583
00584 long int last,
00585 bool *lgEOL)
00586 {
00587 DEBUG_ENTRY( "FFmtRead()" );
00588
00589 const int MAX_SIZE = 1001;
00590 char chr = '\0', chLocal[MAX_SIZE];
00591 const char *eol_ptr = &chCard[last];
00592 const char *ptr = min(&chCard[*ipnt-1],eol_ptr);
00593
00594 ASSERT( *ipnt > 0 && last - *ipnt + 2 <= MAX_SIZE );
00595
00596 while( ptr < eol_ptr && ( chr = *ptr++ ) != '\0' )
00597 {
00598 if( isdigit(chr) || chr == '-' )
00599 break;
00600 }
00601
00602 if( ptr == eol_ptr || chr == '\0' )
00603 {
00604 *ipnt = last+1;
00605 *lgEOL = true;
00606 DEBUG_EXIT( "FFmtRead()" );
00607 return 0.;
00608 }
00609
00610 char* ptr2 = chLocal;
00611
00612 do
00613 {
00614 if( chr != ',' )
00615 *ptr2++ = chr;
00616 if( ptr == eol_ptr )
00617 break;
00618 chr = *ptr++;
00619 }
00620 while( chr != '\0' && !isblank(chr) );
00621 *ptr2 = '\0';
00622
00623 double value;
00625 if( sscanf( chLocal, "%le", &value ) != 1 )
00626 value = 0.;
00627
00628 *ipnt = ptr - chCard;
00629 *lgEOL = false;
00630
00631 DEBUG_EXIT( "FFmtRead()" );
00632 return value;
00633 }
00634 #endif
00635
00636
00637
00638 long nMatch(const char *chKey,
00639 const char *chCard)
00640 {
00641 const char *ptr;
00642 long Match_v;
00643
00644 DEBUG_ENTRY( "nMatch()" );
00645
00646 ASSERT( strlen(chKey) > 0 );
00647
00648 if( ( ptr = strstr( chCard, chKey ) ) == NULL )
00649 {
00650
00651 Match_v = 0L;
00652 }
00653 else
00654 {
00655
00656 Match_v = (long)(ptr-chCard+1);
00657 }
00658
00659 DEBUG_EXIT( "nMatch()" );
00660 return Match_v;
00661 }
00662
00663
00664
00665
00666
00667
00668
00669
00670 double fudge(long int ipnt)
00671 {
00672 double fudge_v;
00673
00674 DEBUG_ENTRY( "fudge()" );
00675
00676 if( ipnt < 0 )
00677 {
00678
00679 fudge_v = fudgec.nfudge;
00680 }
00681 else if( ipnt >= fudgec.nfudge )
00682 {
00683 fprintf( ioQQQ, " FUDGE factor not entered for array number %3ld\n",
00684 ipnt );
00685 puts( "[Stop in fudge]" );
00686 cdEXIT(EXIT_FAILURE);
00687 }
00688 else
00689 {
00690 fudge_v = fudgec.fudgea[ipnt];
00691 }
00692
00693 DEBUG_EXIT( "fudge()" );
00694
00695 return fudge_v;
00696 }
00697
00698
00699
00700 long int GetElem(char *chCARD_CAPS )
00701 {
00702 int i;
00703
00704 DEBUG_ENTRY( "GetElem()" );
00705
00706
00707
00708
00709
00710
00711 for( i=0; i<(int)LIMELM; ++i )
00712 {
00713 if( nMatch( elementnames.chElementNameShort[i], chCARD_CAPS ) )
00714 {
00715
00716 DEBUG_EXIT( "GetElem()" );
00717
00718
00719 return i;
00720 }
00721 }
00722
00723 DEBUG_EXIT( "GetElem()" );
00724
00725
00726 return (-1 );
00727 }
00728
00729
00730
00731
00732 int GetQuote(
00733
00734 char *chLabel,
00735
00736 char *chCard ,
00737
00738
00739 bool lgABORT )
00740 {
00741 char *i0,
00742 *i1,
00743 *iLoc;
00744 size_t len;
00745
00746 DEBUG_ENTRY( "GetQuote()" );
00747
00748
00749
00750
00751
00752
00753 i0 = strchr( input.chOrgCard,'\"' );
00754
00755 if( i0 != NULL )
00756 {
00757
00758
00759 i1 = strchr( i0+1 ,'\"' );
00760 }
00761 else
00762 {
00763 i1 = NULL;
00764 }
00765
00766
00767
00768 if( i0 == NULL || i1 == NULL )
00769 {
00770 if( lgABORT )
00771 {
00772
00773 fprintf( ioQQQ,
00774 " A filename or label must be specified within double quotes, but no quotes were encountered on this command.\n" );
00775 fprintf( ioQQQ, " Name must be surrounded by exactly two double quotes, like \"name.txt\". \n" );
00776 fprintf( ioQQQ, " The line image follows:\n" );
00777 fprintf( ioQQQ, " %s\n", input.chOrgCard);
00778 fprintf( ioQQQ, " Sorry\n" );
00779 puts( "[Stop in getquote]" );
00780 cdEXIT(EXIT_FAILURE);
00781 }
00782 else
00783 {
00784
00785 chLabel[0] = '\0';
00786
00787 return 1;
00788 }
00789 }
00790
00791
00792 len = (size_t)(i1-i0-1);
00793 strncpy(chLabel,i0+1,len);
00794
00795 chLabel[len] = '\0';
00796
00797
00798 iLoc = strchr( chCard ,'\"' );
00799 if( iLoc == NULL )
00800 {
00801 fprintf( ioQQQ, " Insanity in GetQuote - line image follows:\n" );
00802 fprintf( ioQQQ, " %s\n", input.chOrgCard);
00803 puts( "[Stop in getquote]" );
00804 cdEXIT(EXIT_FAILURE);
00805 }
00806
00807
00808
00809 while( i0 <= i1 )
00810 {
00811 *i0 = ' ';
00812 *iLoc = ' ';
00813 ++i0;
00814 ++iLoc;
00815 }
00816
00817 DEBUG_EXIT( "GetQuote()" );
00818
00819
00820 return 0;
00821 }
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834 #if !HAVE_POWI
00835
00836 double powi( double x, long int n )
00837
00838
00839 {
00840 double p;
00841
00842 DEBUG_ENTRY( "powi()" );
00843
00844 if( x == 0 )
00845 return 0.;
00846
00847
00848 if( n < 0 )
00849 {
00850 n = -n;
00851 x = 1/x;
00852 }
00853
00854 p = is_odd(n) ? x : 1;
00855
00856
00857
00858 while( n >>= 1 )
00859 {
00860 x *= x;
00861 if( is_odd(n) )
00862 p *= x;
00863 }
00864
00865
00866
00867 DEBUG_EXIT( "powi()" );
00868
00869 return p;
00870 }
00871
00872 #endif
00873
00874
00875 long ipow( long m, long n )
00876
00877
00878 {
00879 long p;
00880
00881 DEBUG_ENTRY( "ipow()" );
00882
00883 if( m == 0 || (n < 0 && m > 1) )
00884 return 0L;
00885
00886
00887
00888 if( n < 0 )
00889 {
00890 n = -n;
00891 m = 1/m;
00892 }
00893
00894 p = is_odd(n) ? m : 1;
00895
00896
00897
00898 while( n >>= 1 )
00899 {
00900 m *= m;
00901 if( is_odd(n) )
00902 p *= m;
00903 }
00904
00905
00906
00907 DEBUG_EXIT( "ipow()" );
00908
00909 return p;
00910 }
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928 char *PrintEfmt(const char *fmt, double val )
00929 {
00930 static char buf[30];
00931
00932 DEBUG_ENTRY( "PrintEfmt()" );
00933
00934
00935 sprintf(buf, fmt, val);
00936
00937
00938 # ifdef _MSC_VER
00939 {
00940
00941
00942
00943 char *ep , buf2[30];
00944
00945
00946
00947 if( val >= 0.)
00948 {
00949 strcpy(buf2 , " " );
00950 strcat(buf2 , buf);
00951 strcpy( buf , buf2);
00952 }
00953
00954
00955 if((ep = strchr(buf, 'e')) == NULL)
00956 {
00957 ep = strchr(buf, 'E');
00958 }
00959
00960
00961 if(ep != NULL)
00962 {
00963
00964 ep += 2;
00965
00966
00967 *ep = '\0';
00968
00969
00970 ++ep;
00971
00972
00973 strcat( buf, ep );
00974 }
00975 }
00976 # endif
00977
00978 DEBUG_EXIT( "PrintEfmt()" );
00979
00980 return buf;
00981 }
00982
00983 void PrintE82( FILE* ioOUT, double value )
00984 {
00985 double frac , xlog , xfloor , tvalue;
00986 int iExp;
00987
00988 DEBUG_ENTRY( "PrintE82()" );
00989
00990 if( value < 0 )
00991 {
00992 fprintf(ioOUT,"********");
00993 }
00994 else if( value == 0 )
00995 {
00996 fprintf(ioOUT,"0.00E+00");
00997 }
00998 else
00999 {
01000
01001 tvalue = value;
01002 xlog = log10( tvalue );
01003 xfloor = floor(xlog);
01004
01005 frac = tvalue*pow(10.,-xfloor);
01006
01007 iExp = (int)xfloor;
01008 if( frac>9.9945 )
01009 {
01010 frac /= 10.;
01011 iExp += 1;
01012 }
01013
01014 fprintf(ioOUT,"%.2f",frac);
01015
01016 fprintf(ioOUT,"E");
01017
01018 if(iExp>=0 )
01019 {
01020 fprintf(ioOUT,"+");
01021 }
01022 fprintf(ioOUT,"%.2d",iExp);
01023 }
01024
01025 DEBUG_EXIT( "PrintE82()" );
01026
01027 return;
01028 }
01029
01030
01031
01032 void PrintE71( FILE* ioOUT, double value )
01033 {
01034 double frac , xlog , xfloor , tvalue;
01035 int iExp;
01036
01037 DEBUG_ENTRY( "PrintE71()" );
01038
01039 if( value < 0 )
01040 {
01041 fprintf(ioOUT,"*******");
01042 }
01043 else if( value == 0 )
01044 {
01045 fprintf(ioOUT,"0.0E+00");
01046 }
01047 else
01048 {
01049
01050 tvalue = value;
01051 xlog = log10( tvalue );
01052 xfloor = floor(xlog);
01053
01054 frac = tvalue*pow(10.,-xfloor);
01055
01056 iExp = (int)xfloor;
01057 if( frac>9.9945 )
01058 {
01059 frac /= 10.;
01060 iExp += 1;
01061 }
01062
01063 fprintf(ioOUT,"%.1f",frac);
01064
01065 fprintf(ioOUT,"E");
01066
01067 if(iExp>=0 )
01068 {
01069 fprintf(ioOUT,"+");
01070 }
01071 fprintf(ioOUT,"%.2d",iExp);
01072 }
01073
01074 DEBUG_EXIT( "PrintE71()" );
01075
01076 return;
01077 }
01078
01079
01080
01081
01082 void PrintE93( FILE* ioOUT, double value )
01083 {
01084 double frac , xlog , xfloor, tvalue;
01085 int iExp;
01086
01087 DEBUG_ENTRY( "PrintE93()" );
01088
01089 if( value < 0 )
01090 {
01091 fprintf(ioOUT,"*********");
01092 }
01093 else if( value == 0 )
01094 {
01095 fprintf(ioOUT,"0.000E+00");
01096 }
01097 else
01098 {
01099
01100 tvalue = value;
01101 xlog = log10( tvalue );
01102 xfloor = floor(xlog);
01103
01104 frac = tvalue*pow(10.,-xfloor);
01105
01106 iExp = (int)xfloor;
01107 if( frac>9.99949 )
01108 {
01109 frac /= 10.;
01110 iExp += 1;
01111 }
01112
01113 fprintf(ioOUT,"%5.3f",frac);
01114
01115 fprintf(ioOUT,"E");
01116
01117 if(iExp>=0 )
01118 {
01119 fprintf(ioOUT,"+");
01120 }
01121 fprintf(ioOUT,"%.2d",iExp);
01122 }
01123
01124 DEBUG_EXIT( "PrintE93()" );
01125
01126 return;
01127 }
01128
01129
01130 NORETURN void TotalInsanity(void)
01131 {
01132
01133 DEBUG_ENTRY( "TotalInsanity()" );
01134
01135
01136
01137
01138 fprintf( ioQQQ, " Something that cannot happen, has happened.\n" );
01139 fprintf( ioQQQ, " This is TotalInsanity, I live in service.c.\n" );
01140 ShowMe();
01141
01142 puts( "[Stop in TotalInsanity]" );
01143
01144 DEBUG_EXIT( "TotalInsanity()" );
01145
01146 cdEXIT(EXIT_FAILURE);
01147 }
01148
01149
01150
01151 NORETURN void BadRead(void)
01152 {
01153
01154 DEBUG_ENTRY( "BadRead()" );
01155
01156
01157 fprintf( ioQQQ, " A read of internal input data has failed.\n" );
01158 fprintf( ioQQQ, " This is BadRead, I live in service.c.\n" );
01159 ShowMe();
01160
01161 puts( "[Stop in BadRead]" );
01162
01163 DEBUG_EXIT( "BadRead()" );
01164
01165 cdEXIT(EXIT_FAILURE);
01166 }
01167
01168
01169 NORETURN void BadOpen(void)
01170 {
01171
01172 DEBUG_ENTRY( "BadOpen()" );
01173
01174
01175 fprintf( ioQQQ, " An attempt at opening a files has failed.\n" );
01176 fprintf( ioQQQ, " This is BadOpen, I live in service.c.\n" );
01177 ShowMe();
01178
01179 puts( "[Stop in BadOpen]" );
01180
01181 DEBUG_EXIT( "BadOpen()" );
01182
01183 cdEXIT(EXIT_FAILURE);
01184 }
01185
01186
01187 NORETURN void NoNumb(char *chCard)
01188 {
01189
01190 DEBUG_ENTRY( "NoNumb()" );
01191
01192
01193 fprintf( ioQQQ, " There is a problem on the following command line:\n" );
01194 fprintf( ioQQQ, " %s\n", chCard );
01195 fprintf( ioQQQ, " There should have been a number on this line. Sorry.\n" );
01196 puts( "[Stop in nonumb]" );
01197
01198 DEBUG_EXIT( "NoNumb()" );
01199
01200 cdEXIT(EXIT_FAILURE);
01201 }
01202
01203
01204 double sexp(double x)
01205 {
01206 double sexp_v;
01207
01208 DEBUG_ENTRY( "sexp()" );
01209
01210
01211 if( x < SEXP_LIMIT )
01212 {
01213 sexp_v = exp(-x);
01214 }
01215 else
01216 {
01217 sexp_v = 0.;
01218 }
01219
01220 DEBUG_EXIT( "sexp()" );
01221
01222 return sexp_v;
01223 }
01224
01225
01226
01227 double dsexp(double x)
01228 {
01229 double dsexp_v;
01230
01231 DEBUG_ENTRY( "dsexp()" );
01232
01233 if( x < 680. )
01234 {
01235 dsexp_v = exp(-x);
01236 }
01237 else
01238 {
01239 dsexp_v = 0.;
01240 }
01241
01242 DEBUG_EXIT( "dsexp()" );
01243
01244 return dsexp_v;
01245 }
01246
01247
01248
01249 void TestCode(void)
01250 {
01251
01252 DEBUG_ENTRY( "TestCode( )" );
01253
01254
01255 lgTestCodeCalled = true;
01256
01257 DEBUG_EXIT( "TestCode( )" );
01258
01259 return;
01260 }
01261
01262
01263 void broken(void)
01264 {
01265
01266 broke.lgBroke = true;
01267
01268 return;
01269 }
01270
01271
01272 void fixit(void)
01273 {
01274
01275 broke.lgFixit = true;
01276
01277 return;
01278 }
01279
01280
01281 void CodeReview(void)
01282 {
01283
01284 broke.lgCheckit = true;
01285
01286 return;
01287 }
01288
01289
01290
01291
01292 void dbg_printf(int debug, const char *fmt, ...)
01293 {
01294 va_list ap;
01295
01296 DEBUG_ENTRY( "dbg_printf()" );
01297
01298
01299 if(debug <= trace.debug_level)
01300 {
01301 va_start(ap, fmt);
01302
01303
01304 vfprintf(ioQQQ, fmt, ap);
01305
01306
01307 fflush(ioQQQ);
01308 va_end(ap);
01309 }
01310 DEBUG_EXIT( "dbg_printf()" );
01311
01312 return;
01313 }
01314
01315
01316
01317 double qg32(
01318 double xl,
01319 double xu,
01320
01321 double (*fct)(double) )
01322 {
01323 double a,
01324 b,
01325 c,
01326 y;
01327
01328 DEBUG_ENTRY( "qg32()" );
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344 a = .5*(xu + xl);
01345 b = xu - xl;
01346 c = .498631930924740780*b;
01347 y = .35093050047350483e-2*((*fct)(a+c) + (*fct)(a-c));
01348 c = b*.49280575577263417;
01349 y += .8137197365452835e-2*((*fct)(a+c) + (*fct)(a-c));
01350 c = b*.48238112779375322;
01351 y += .1269603265463103e-1*((*fct)(a+c) + (*fct)(a-c));
01352 c = b*.46745303796886984;
01353 y += .17136931456510717e-1*((*fct)(a+c) + (*fct)(a-c));
01354 c = b*.44816057788302606;
01355 y += .21417949011113340e-1*((*fct)(a+c) + (*fct)(a-c));
01356 c = b*.42468380686628499;
01357 y += .25499029631188088e-1*((*fct)(a+c) + (*fct)(a-c));
01358 c = b*.3972418979839712;
01359 y += .29342046739267774e-1*((*fct)(a+c) + (*fct)(a-c));
01360 c = b*.36609105937014484;
01361 y += .32911111388180923e-1*((*fct)(a+c) + (*fct)(a-c));
01362 c = b*.3315221334651076;
01363 y += .36172897054424253e-1*((*fct)(a+c) + (*fct)(a-c));
01364 c = b*.29385787862038116;
01365 y += .39096947893535153e-1*((*fct)(a+c) + (*fct)(a-c));
01366 c = b*.2534499544661147;
01367 y += .41655962113473378e-1*((*fct)(a+c) + (*fct)(a-c));
01368 c = b*.21067563806531767;
01369 y += .43826046502201906e-1*((*fct)(a+c) + (*fct)(a-c));
01370 c = b*.16593430114106382;
01371 y += .45586939347881942e-1*((*fct)(a+c) + (*fct)(a-c));
01372 c = b*.11964368112606854;
01373 y += .46922199540402283e-1*((*fct)(a+c) + (*fct)(a-c));
01374 c = b*.7223598079139825e-1;
01375 y += .47819360039637430e-1*((*fct)(a+c) + (*fct)(a-c));
01376 c = b*.24153832843869158e-1;
01377 y = b*(y + .482700442573639e-1*((*fct)(a+c) + (*fct)(a-c)));
01378
01379 DEBUG_EXIT( "qg32()" );
01380
01381
01382 return y;
01383 }
01384
01385
01386 void spsort(
01387
01388 float x[],
01389
01390 long int n,
01391
01392 long int iperm[],
01393
01394
01395 int kflag,
01396
01397 int *ier)
01398 {
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469
01470 long int i,
01471 ij,
01472 il[21],
01473 indx,
01474 indx0,
01475 istrt,
01476 istrt_,
01477 iu[21],
01478 j,
01479 k,
01480 kk,
01481 l,
01482 lm,
01483 lmt,
01484 m,
01485 nn;
01486 float r,
01487 ttemp;
01488
01489 DEBUG_ENTRY( "spsort()" );
01490
01491
01492
01493
01494
01495
01496
01497 *ier = 0;
01498 nn = n;
01499 if( nn < 1 )
01500 {
01501 *ier = 1;
01502
01503 DEBUG_EXIT( "spsort()" );
01504 return;
01505 }
01506 else
01507 {
01508 kk = labs(kflag);
01509 if( kk != 1 && kk != 2 )
01510 {
01511 *ier = 2;
01512
01513 DEBUG_EXIT( "spsort()" );
01514 return;
01515 }
01516 else
01517 {
01518
01519
01520
01521 for( i=0; i < nn; i++ )
01522 {
01523 iperm[i] = i+1;
01524 }
01525
01526
01527
01528 if( nn == 1 )
01529 {
01530 --iperm[0];
01531 DEBUG_EXIT( "spsort()" );
01532 return;
01533 }
01534
01535
01536
01537 if( kflag <= -1 )
01538 {
01539 for( i=0; i < nn; i++ )
01540 {
01541 x[i] = -x[i];
01542 }
01543 }
01544
01545
01546
01547 m = 1;
01548 i = 1;
01549 j = nn;
01550 r = .375e0;
01551 }
01552 }
01553
01554 while( true )
01555 {
01556 if( i == j )
01557 goto L_80;
01558 if( r <= 0.5898437e0 )
01559 {
01560 r += 3.90625e-2;
01561 }
01562 else
01563 {
01564 r -= 0.21875e0;
01565 }
01566
01567 L_40:
01568 k = i;
01569
01570
01571
01572 ij = i + (long)((j-i)*r);
01573 lm = iperm[ij-1];
01574
01575
01576
01577 if( x[iperm[i-1]-1] > x[lm-1] )
01578 {
01579 iperm[ij-1] = iperm[i-1];
01580 iperm[i-1] = lm;
01581 lm = iperm[ij-1];
01582 }
01583 l = j;
01584
01585
01586
01587 if( x[iperm[j-1]-1] < x[lm-1] )
01588 {
01589 iperm[ij-1] = iperm[j-1];
01590 iperm[j-1] = lm;
01591 lm = iperm[ij-1];
01592
01593
01594
01595
01596 if( x[iperm[i-1]-1] > x[lm-1] )
01597 {
01598 iperm[ij-1] = iperm[i-1];
01599 iperm[i-1] = lm;
01600 lm = iperm[ij-1];
01601 }
01602 }
01603
01604
01605
01606
01607 while( true )
01608 {
01609 l -= 1;
01610 if( x[iperm[l-1]-1] <= x[lm-1] )
01611 {
01612
01613
01614
01615
01616 while( true )
01617 {
01618 k += 1;
01619 if( x[iperm[k-1]-1] >= x[lm-1] )
01620 break;
01621 }
01622
01623
01624
01625 if( k > l )
01626 break;
01627 lmt = iperm[l-1];
01628 iperm[l-1] = iperm[k-1];
01629 iperm[k-1] = lmt;
01630 }
01631 }
01632
01633
01634
01635 if( l - i > j - k )
01636 {
01637 il[m-1] = i;
01638 iu[m-1] = l;
01639 i = k;
01640 m += 1;
01641 }
01642 else
01643 {
01644 il[m-1] = k;
01645 iu[m-1] = j;
01646 j = l;
01647 m += 1;
01648 }
01649
01650 L_90:
01651 if( j - i >= 1 )
01652 goto L_40;
01653 if( i == 1 )
01654 continue;
01655 i -= 1;
01656
01657 while( true )
01658 {
01659 i += 1;
01660 if( i == j )
01661 break;
01662 lm = iperm[i];
01663 if( x[iperm[i-1]-1] > x[lm-1] )
01664 {
01665 k = i;
01666
01667 while( true )
01668 {
01669 iperm[k] = iperm[k-1];
01670 k -= 1;
01671
01672 if( x[lm-1] >= x[iperm[k-1]-1] )
01673 break;
01674 }
01675 iperm[k] = lm;
01676 }
01677 }
01678
01679
01680
01681 L_80:
01682 m -= 1;
01683 if( m == 0 )
01684 break;
01685
01686 i = il[m-1];
01687 j = iu[m-1];
01688
01689 goto L_90;
01690 }
01691
01692
01693
01694 if( kflag <= -1 )
01695 {
01696 for( i=0; i < nn; i++ )
01697 {
01698 x[i] = -x[i];
01699 }
01700 }
01701
01702
01703
01704 if( kk == 2 )
01705 {
01706
01707
01708
01709
01710 for( istrt=1; istrt <= nn; istrt++ )
01711 {
01712 istrt_ = istrt - 1;
01713 if( iperm[istrt_] >= 0 )
01714 {
01715 indx = istrt;
01716 indx0 = indx;
01717 ttemp = x[istrt_];
01718 while( iperm[indx-1] > 0 )
01719 {
01720 x[indx-1] = x[iperm[indx-1]-1];
01721 indx0 = indx;
01722 iperm[indx-1] = -iperm[indx-1];
01723 indx = labs(iperm[indx-1]);
01724 }
01725 x[indx0-1] = ttemp;
01726 }
01727 }
01728
01729
01730
01731 for( i=0; i < nn; i++ )
01732 {
01733 iperm[i] = -iperm[i];
01734 }
01735 }
01736
01737
01738 for( i=0; i < nn; i++ )
01739 {
01740 --iperm[i];
01741 }
01742
01743 DEBUG_EXIT( "spsort()" );
01744
01745 return;
01746 }
01747
01748
01749
01750
01751
01752
01753
01754
01755 void *MyMalloc(
01756
01757 size_t size ,
01758 const char *chFile,
01759 int line
01760 )
01761 {
01762 void *ptr;
01763
01764 DEBUG_ENTRY( "MyMalloc()" );
01765
01766 ASSERT( size > 0 );
01767
01768
01769 {
01770
01771 enum{DEBUG_LOC=false};
01772
01773 if( DEBUG_LOC)
01774 {
01775 static long int kount=0, nTot=0;
01776 nTot += (long)size;
01777 fprintf(ioQQQ,"%li\t%li\t%li\n",
01778 kount ,
01779 (long)size ,
01780 nTot );
01781 ++kount;
01782 }
01783 }
01784
01785 if( ( ptr = malloc( size ) ) == NULL )
01786 {
01787 fprintf(ioQQQ,"DISASTER MyMalloc could not allocate %lu bytes. Exit in MyMalloc.",
01788 (unsigned long)size );
01789 fprintf(ioQQQ,"MyMalloc called from file %s at line %i.\n",
01790 chFile , line );
01791 puts( "[Stop in MyMalloc]" );
01792 cdEXIT(EXIT_FAILURE);
01793 }
01794
01795
01796 # if !defined(NDEBUG) && !defined(NOINIT)
01797
01798 size_t nFloat = size/4;
01799 size_t nDouble = size/8;
01800 float *fptr = static_cast<float*>(ptr);
01801 double *dptr = static_cast<double*>(ptr);
01802
01803
01804
01805
01806
01807 if( size == nDouble*8 )
01808 {
01809
01810
01811
01812
01813
01814
01815
01816 set_nan( dptr, (long)nDouble );
01817 }
01818 else if( size == nFloat*4 )
01819 {
01820
01821 set_nanf( fptr, (long)nFloat );
01822 }
01823 else
01824 {
01825 memset( ptr, 0xff, size );
01826 }
01827
01828 # endif
01829
01830 DEBUG_EXIT( "MyMalloc()" );
01831
01832 return ptr;
01833 }
01834
01835
01836
01837
01838
01839
01840 void *MyCalloc(
01841
01842 size_t num ,
01843 size_t size )
01844 {
01845 void *ptr;
01846
01847 DEBUG_ENTRY( "MyCalloc()" );
01848
01849 ASSERT( size > 0 );
01850
01851
01852 {
01853
01854 enum{DEBUG_LOC=false};
01855
01856 if( DEBUG_LOC)
01857 {
01858 static long int kount=0;
01859 fprintf(ioQQQ,"%li\tcall\t%li\tbytes\n", kount,
01860 (long)size );
01861 ++kount;
01862 }
01863 }
01864
01865 if( ( ptr = calloc( num , size ) ) == NULL )
01866 {
01867 fprintf(ioQQQ,"MyCalloc could not allocate %lu bytes. Exit in MyCalloc.",
01868 (unsigned long)size );
01869 puts( "[Stop in MyCalloc]" );
01870 cdEXIT(EXIT_FAILURE);
01871 }
01872
01873 DEBUG_EXIT( "MyCalloc()" );
01874
01875 return ptr;
01876 }
01877
01878
01879
01880
01881
01882 void *MyRealloc(
01883
01884 void *p ,
01885 size_t size )
01886 {
01887 void *ptr;
01888
01889 DEBUG_ENTRY( "MyRealloc()" );
01890
01891 ASSERT( size > 0 );
01892
01893
01894 {
01895
01896 enum{DEBUG_LOC=false};
01897
01898 if( DEBUG_LOC)
01899 {
01900 static long int kount=0;
01901 fprintf(ioQQQ,"%li\tcall\t%li\tbytes\n", kount,
01902 (long)size );
01903 ++kount;
01904 }
01905 }
01906
01907 if( ( ptr = realloc( p , size ) ) == NULL )
01908 {
01909 fprintf(ioQQQ,"MyRealloc could not allocate %lu bytes. Exit in MyRealloc.",
01910 (unsigned long)size );
01911 puts( "[Stop in MyRealloc]" );
01912 cdEXIT(EXIT_FAILURE);
01913 }
01914
01915 DEBUG_EXIT( "MyRealloc()" );
01916
01917 return ptr;
01918 }
01919
01920
01921 double csphot(
01922
01923
01924 long int inu,
01925
01926 long int ithr,
01927
01928 long int iofset)
01929 {
01930 double csphot_v;
01931
01932 DEBUG_ENTRY( "csphot()" );
01933
01934 csphot_v = opac.OpacStack[inu-ithr+iofset-1];
01935
01936 DEBUG_EXIT( "csphot()" );
01937
01938 return csphot_v;
01939 }
01940
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960
01961
01962
01963
01964
01965
01966
01967
01968
01969 double RandGauss(
01970
01971 double xMean,
01972
01973 double s )
01974 {
01975 double x1, x2, w, yy1;
01976 static double yy2=-BIGDOUBLE;
01977 static int use_last = false;
01978
01979 DEBUG_ENTRY( "RandGauss()" );
01980
01981 if( use_last )
01982 {
01983 yy1 = yy2;
01984 use_last = false;
01985 }
01986 else
01987 {
01988 do {
01989 x1 = 2.*((double)rand()/ (double)RAND_MAX) - 1.;
01990 x2 = 2.*((double)rand() / (double)RAND_MAX) - 1.;
01991 w = x1 * x1 + x2 * x2;
01992 } while ( w >= 1.0 );
01993
01994 w = sqrt((-2.0*log(w))/w);
01995 yy1 = x1 * w;
01996 yy2 = x2 * w;
01997 use_last = true;
01998 }
01999
02000 DEBUG_EXIT( "RandGauss()" );
02001
02002 return xMean + yy1 * s;
02003 }
02004
02005
02006
02007
02008
02009
02010
02011
02012 double MyGaussRand( double PctUncertainty )
02013 {
02014 double StdDev;
02015 double result;
02016
02017 DEBUG_ENTRY( "MyGaussRand()" );
02018
02019 ASSERT( PctUncertainty < 0.5 );
02020
02021 StdDev = 0.5 * PctUncertainty;
02022
02023 do
02024 {
02025
02026 result = 1.+RandGauss( 0., StdDev );
02027 }
02028
02029
02030 while( (result < 1.-PctUncertainty) || (result > 1+PctUncertainty) );
02031
02032 ASSERT( result>0. && result<2. );
02033
02034 DEBUG_EXIT( "MyGaussRand()" );
02035
02036 return result;
02037 }
02038
02039
02040 double plankf(long int ip)
02041 {
02042 double plankf_v;
02043
02044 DEBUG_ENTRY( "plankf()" );
02045
02046
02047
02048
02049 if( rfield.ContBoltz[ip] <= 0. )
02050 {
02051 plankf_v = 1e-35;
02052 }
02053 else
02054 {
02055 plankf_v = 6.991e-21*POW2(FR1RYD*rfield.anu[ip])/
02056 (1./rfield.ContBoltz[ip] - 1.)*FR1RYD*4.;
02057 }
02058
02059 DEBUG_EXIT( "plankf()" );
02060
02061 return plankf_v;
02062 }
02063