77 if( (chRet = fgets( chLine , nChar, ioIN )) == NULL )
84 chEOL = (
char*)memchr( chLine ,
'\0' , nChar );
89 if( (chEOL==NULL) || (chEOL - chLine)>=nChar-1 )
92 fprintf(
ioQQQ,
"DISASTER PROBLEM read_whole_line found input"
93 " with a line too long to be read.\n" );
112 string::size_type ptr1 = 0;
113 string::size_type ptr2 = str.find( sep );
114 string sstr = str.substr( ptr1, ptr2-ptr1 );
115 if( sstr.length() > 0 )
116 lst.push_back( sstr );
118 if( lgStrict ) lgFail =
true;
119 if( lgKeep ) lst.push_back( sstr );
121 while( ptr2 != string::npos ) {
123 ptr1 = ptr2 + sep.length();
124 if( ptr1 < str.length() ) {
125 ptr2 = str.find( sep, ptr1 );
126 sstr = str.substr( ptr1, ptr2-ptr1 );
127 if( sstr.length() > 0 )
128 lst.push_back( sstr );
130 if( lgStrict ) lgFail =
true;
131 if( lgKeep ) lst.push_back( sstr );
136 if( lgStrict ) lgFail =
true;
137 if( lgKeep ) lst.push_back(
"" );
142 fprintf(
ioQQQ,
" A syntax error occurred while splitting the string: \"%s\"\n", str.c_str() );
143 fprintf(
ioQQQ,
" The separator is \"%s\". Empty substrings are not allowed.\n", sep.c_str() );
153 fprintf(
ioQQQ,
" PROBLEM DISASTER An assert has been thrown, this is bad.\n");
154 fprintf(
ioQQQ,
" It happened in the file %s at line number %i\n", file, line );
155 fprintf(
ioQQQ,
" This is iteration %li, nzone %li, fzone %.2f, lgSearch=%c.\n",
162 # if defined ASSERTDEBUG
176 if( energy_ryd <=0. )
184 AnuUnit_v = energy_ryd;
189 AnuUnit_v =
RYDLAM/energy_ryd*1e-4;
194 AnuUnit_v = energy_ryd*
EVRYD*1.e-3;
199 AnuUnit_v = energy_ryd*
EVRYD;
204 AnuUnit_v =
RYDLAM/energy_ryd;
209 AnuUnit_v =
RYDLAM/energy_ryd*1e-8;
214 AnuUnit_v =
RYD_INF*energy_ryd;
223 fprintf(
ioQQQ,
" insane units in AnuUnit =%4.4s\n",
243 fprintf(
ioQQQ,
" >>> \n >>> \n >>> Cosmic rays are not included and the gas is molecular. "
244 "THIS IS KNOWN TO BE UNSTABLE. Add cosmic rays and try again.\n >>> \n >>>\n\n");
248 fprintf(
ioQQQ,
"\n\n\n" );
249 fprintf(
ioQQQ,
" vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv \n" );
250 fprintf(
ioQQQ,
" > PROBLEM DISASTER PROBLEM DISASTER. <\n" );
251 fprintf(
ioQQQ,
" > Sorry, something bad has happened. <\n" );
252 fprintf(
ioQQQ,
" > Please post this on the Cloudy web site <\n" );
253 fprintf(
ioQQQ,
" > discussion board at www.nublado.org <\n" );
254 fprintf(
ioQQQ,
" > Please send all following information: <\n" );
255 fprintf(
ioQQQ,
" ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ \n" );
256 fprintf(
ioQQQ,
"\n\n" );
259 fprintf(
ioQQQ,
" Cloudy version number is %s\n",
263 fprintf(
ioQQQ,
"%5ld warnings,%3ld cautions,%3ld temperature failures. Messages follow.\n",
279 fprintf(
ioQQQ,
" This input stream included an init file.\n");
280 fprintf(
ioQQQ,
" If this init file is not part of the standard Cloudy distribution\n");
281 fprintf(
ioQQQ,
" then I will need a copy of it too.\n");
300 for( i=0; i < 4; i++ )
303 chCAP[i] = (char)toupper( chLab[i] );
320 while( chCard[i]!=
'\0' )
322 chCard[i] = (char)tolower( chCard[i] );
337 while( chCard[i]!=
'\0' )
339 chCard[i] = (char)toupper( chCard[i] );
354 double hold =
sexp(t) - t*
ee1(t);
357 return max( hold, 0. );
366 static double a[6]={-.57721566,.99999193,-.24991055,.05519968,-.00976004,
368 static double b[4]={8.5733287401,18.0590169730,8.6347608925,.2677737343};
369 static double c[4]={9.5733223454,25.6329561486,21.0996530827,3.9584969228};
382 fprintf(
ioQQQ,
" DISASTER negative argument in function ee1, x<0\n" );
390 ans = ((((a[5]*x + a[4])*x + a[3])*x + a[2])*x + a[1])*x + a[0] - log(x);
397 top = (((x + b[0])*x + b[1])*x + b[2])*x + b[3];
398 bot = (((x + c[0])*x + c[1])*x + c[2])*x + c[3];
399 ans = top/bot/x*exp(-x);
412 static double b[4]={8.5733287401,18.0590169730,8.6347608925,.2677737343};
413 static double c[4]={9.5733223454,25.6329561486,21.0996530827,3.9584969228};
421 top = (((x + b[0])*x + b[1])*x + b[2])*x + b[3];
423 bot = (((x + c[0])*x + c[1])*x + c[2])*x + c[3];
438 const int MAX_SIZE = 1001;
439 char chr =
'\0', chLocal[MAX_SIZE];
440 const char *eol_ptr = &chCard[last];
441 const char *ptr =
min(&chCard[*ipnt-1],eol_ptr);
443 ASSERT( *ipnt > 0 && last - *ipnt + 2 <= MAX_SIZE );
445 while( ptr < eol_ptr && ( chr = *ptr++ ) !=
'\0' )
447 const char *lptr = ptr;
449 if( lchr ==
'-' || lchr ==
'+' )
457 if( ptr == eol_ptr || chr ==
'\0' )
465 char* ptr2 = chLocal;
485 while( isdigit(chr) || chr ==
'.' || chr ==
'-' || chr ==
'+' || chr ==
',' || chr ==
'e' || chr ==
'E' );
491 double value = atof( chLocal );
493 *ipnt = (long)(ptr - chCard);
508 ASSERT( strlen(chKey) > 0 );
510 if( ( ptr = strstr( chCard, chKey ) ) == NULL )
518 Match_v = (long)(ptr-chCard+1);
544 fprintf(
ioQQQ,
" FUDGE factor not entered for array number %3ld\n",
569 for( i=0; i<(int)
LIMELM; ++i )
611 i1 = strchr( i0+1 ,
'\"' );
620 if( i0 == NULL || i1 == NULL )
626 " A filename or label must be specified within double quotes, but no quotes were encountered on this command.\n" );
627 fprintf(
ioQQQ,
" Name must be surrounded by exactly two double quotes, like \"name.txt\". \n" );
628 fprintf(
ioQQQ,
" The line image follows:\n" );
630 fprintf(
ioQQQ,
" Sorry\n" );
643 len = (size_t)(i1-i0-1);
644 strncpy(chLabel,i0+1,len);
649 iLoc = strchr( chCard ,
'\"' );
652 fprintf(
ioQQQ,
" Insanity in GetQuote - line image follows:\n" );
683 double powi(
double x,
long int n )
726 if( m == 0 || (n < 0 && m > 1) )
775 sprintf(buf, fmt, val);
795 if((ep = strchr(buf,
'e')) == NULL)
797 ep = strchr(buf,
'E');
823 double frac , xlog , xfloor , tvalue;
830 fprintf(ioOUT,
"********");
832 else if( value == 0 )
834 fprintf(ioOUT,
"0.00E+00");
840 xlog = log10( tvalue );
841 xfloor = floor(xlog);
843 frac = tvalue*pow(10.,-xfloor);
852 fprintf(ioOUT,
"%.2f",frac);
860 fprintf(ioOUT,
"%.2d",iExp);
869 double frac , xlog , xfloor , tvalue;
876 fprintf(ioOUT,
"*******");
878 else if( value == 0 )
880 fprintf(ioOUT,
"0.0E+00");
886 xlog = log10( tvalue );
887 xfloor = floor(xlog);
889 frac = tvalue*pow(10.,-xfloor);
898 fprintf(ioOUT,
"%.1f",frac);
906 fprintf(ioOUT,
"%.2d",iExp);
916 double frac , xlog , xfloor, tvalue;
923 fprintf(ioOUT,
"*********");
925 else if( value == 0 )
927 fprintf(ioOUT,
"0.000E+00");
933 xlog = log10( tvalue );
934 xfloor = floor(xlog);
936 frac = tvalue*pow(10.,-xfloor);
945 fprintf(ioOUT,
"%5.3f",frac);
953 fprintf(ioOUT,
"%.2d",iExp);
967 fprintf(
ioQQQ,
" Something that cannot happen, has happened.\n" );
968 fprintf(
ioQQQ,
" This is TotalInsanity, I live in service.cpp.\n" );
982 fprintf(
ioQQQ,
" A read of internal input data has failed.\n" );
983 fprintf(
ioQQQ,
" This is BadRead, I live in service.c.\n" );
996 fprintf(
ioQQQ,
" An attempt at opening a files has failed.\n" );
997 fprintf(
ioQQQ,
" This is BadOpen, I live in service.c.\n" );
1010 fprintf(
ioQQQ,
" There is a problem on the following command line:\n" );
1011 fprintf(
ioQQQ,
" %s\n", chCard );
1012 fprintf(
ioQQQ,
" There should have been a number on this line. Sorry.\n" );
1120 va_start(ap,format);
1121 i1 = fprintf(fp,
"DEBUG ");
1123 i2 = vfprintf(fp,format,ap);
1148 i = vfprintf(
ioQQQ, fmt, ap);
1162 double (*fct)(
double) )
1186 c = .498631930924740780*b;
1187 y = .35093050047350483e-2*((*fct)(a+c) + (*fct)(a-c));
1188 c = b*.49280575577263417;
1189 y += .8137197365452835e-2*((*fct)(a+c) + (*fct)(a-c));
1190 c = b*.48238112779375322;
1191 y += .1269603265463103e-1*((*fct)(a+c) + (*fct)(a-c));
1192 c = b*.46745303796886984;
1193 y += .17136931456510717e-1*((*fct)(a+c) + (*fct)(a-c));
1194 c = b*.44816057788302606;
1195 y += .21417949011113340e-1*((*fct)(a+c) + (*fct)(a-c));
1196 c = b*.42468380686628499;
1197 y += .25499029631188088e-1*((*fct)(a+c) + (*fct)(a-c));
1198 c = b*.3972418979839712;
1199 y += .29342046739267774e-1*((*fct)(a+c) + (*fct)(a-c));
1200 c = b*.36609105937014484;
1201 y += .32911111388180923e-1*((*fct)(a+c) + (*fct)(a-c));
1202 c = b*.3315221334651076;
1203 y += .36172897054424253e-1*((*fct)(a+c) + (*fct)(a-c));
1204 c = b*.29385787862038116;
1205 y += .39096947893535153e-1*((*fct)(a+c) + (*fct)(a-c));
1206 c = b*.2534499544661147;
1207 y += .41655962113473378e-1*((*fct)(a+c) + (*fct)(a-c));
1208 c = b*.21067563806531767;
1209 y += .43826046502201906e-1*((*fct)(a+c) + (*fct)(a-c));
1210 c = b*.16593430114106382;
1211 y += .45586939347881942e-1*((*fct)(a+c) + (*fct)(a-c));
1212 c = b*.11964368112606854;
1213 y += .46922199540402283e-1*((*fct)(a+c) + (*fct)(a-c));
1214 c = b*.7223598079139825e-1;
1215 y += .47819360039637430e-1*((*fct)(a+c) + (*fct)(a-c));
1216 c = b*.24153832843869158e-1;
1217 y = b*(y + .482700442573639e-1*((*fct)(a+c) + (*fct)(a-c)));
1344 if( kk != 1 && kk != 2 )
1354 for( i=0; i < nn; i++ )
1369 for( i=0; i < nn; i++ )
1387 if( r <= 0.5898437e0 )
1401 ij = i + (long)((j-i)*r);
1406 if( x[iperm[i-1]-1] > x[lm-1] )
1408 iperm[ij-1] = iperm[i-1];
1416 if( x[iperm[j-1]-1] < x[lm-1] )
1418 iperm[ij-1] = iperm[j-1];
1425 if( x[iperm[i-1]-1] > x[lm-1] )
1427 iperm[ij-1] = iperm[i-1];
1438 if( x[iperm[l-1]-1] <= x[lm-1] )
1446 if( x[iperm[k-1]-1] >= x[lm-1] )
1454 iperm[l-1] = iperm[k-1];
1488 if( x[iperm[i-1]-1] > x[lm-1] )
1494 iperm[k] = iperm[k-1];
1497 if( x[lm-1] >= x[iperm[k-1]-1] )
1519 for( i=0; i < nn; i++ )
1531 for( istrt=1; istrt <= nn; istrt++ )
1534 if( iperm[istrt_] >= 0 )
1539 while( iperm[indx-1] > 0 )
1541 x[indx-1] = x[iperm[indx-1]-1];
1543 iperm[indx-1] = -iperm[indx-1];
1544 indx = labs(iperm[indx-1]);
1551 for( i=0; i < nn; i++ )
1553 iperm[i] = -iperm[i];
1557 for( i=0; i < nn; i++ )
1587 enum{DEBUG_LOC=
false};
1591 static long int kount=0, nTot=0;
1593 fprintf(
ioQQQ,
"%li\t%li\t%li\n",
1601 if( ( ptr = malloc( size ) ) == NULL )
1603 fprintf(
ioQQQ,
"DISASTER MyMalloc could not allocate %lu bytes. Exit in MyMalloc.",
1604 (
unsigned long)size );
1605 fprintf(
ioQQQ,
"MyMalloc called from file %s at line %i.\n",
1611 # if !defined(NDEBUG) && !defined(NOINIT)
1613 size_t nFloat = size/4;
1614 size_t nDouble = size/8;
1616 double *dptr =
static_cast<double*
>(ptr);
1622 if( size == nDouble*8 )
1631 set_NaN( dptr, (
long)nDouble );
1633 else if( size == nFloat*4 )
1636 set_NaN( fptr, (
long)nFloat );
1640 memset( ptr, 0xff, size );
1666 enum{DEBUG_LOC=
false};
1670 static long int kount=0;
1671 fprintf(
ioQQQ,
"%li\tcall\t%li\tbytes\n", kount,
1677 if( ( ptr = calloc( num , size ) ) == NULL )
1679 fprintf(
ioQQQ,
"MyCalloc could not allocate %lu bytes. Exit in MyCalloc.",
1680 (
unsigned long)size );
1704 enum{DEBUG_LOC=
false};
1708 static long int kount=0;
1709 fprintf(
ioQQQ,
"%li\tcall\t%li\tbytes\n", kount,
1715 if( ( ptr = realloc( p , size ) ) == NULL )
1717 fprintf(
ioQQQ,
"MyRealloc could not allocate %lu bytes. Exit in MyRealloc.",
1718 (
unsigned long)size );
1772 double x1,
x2, w, yy1;
1774 static int use_last =
false;
1788 w = x1 * x1 + x2 *
x2;
1789 }
while ( w >= 1.0 );
1791 w = sqrt((-2.0*log(w))/w);
1796 return xMean + yy1 * s;
1813 ASSERT( PctUncertainty < 0.5 );
1815 StdDev = 0.5 * PctUncertainty;
1824 while( (result < 1.-PctUncertainty) || (result > 1+PctUncertainty) );
1826 ASSERT( result>0. && result<2. );