77 memset( chLine, 0, (
size_t)nChar );
81 if( (chRet = fgets( chLine, nChar, ioIN )) == NULL )
84 long lineLength = strlen( chRet );
91 if( lineLength>=nChar-1 )
94 fprintf(
ioQQQ,
"DISASTER PROBLEM read_whole_line found input"
95 " with a line too long to be read, limit is %i char. "
96 "Start of line follows:\n%s\n",
116 string::size_type ptr1 = 0;
117 string::size_type ptr2 = str.find( sep );
118 string sstr = str.substr( ptr1, ptr2-ptr1 );
119 if( sstr.length() > 0 )
120 lst.push_back( sstr );
122 if( lgStrict ) lgFail =
true;
123 if( lgKeep ) lst.push_back( sstr );
125 while( ptr2 != string::npos ) {
127 ptr1 = ptr2 + sep.length();
128 if( ptr1 < str.length() ) {
129 ptr2 = str.find( sep, ptr1 );
130 sstr = str.substr( ptr1, ptr2-ptr1 );
131 if( sstr.length() > 0 )
132 lst.push_back( sstr );
134 if( lgStrict ) lgFail =
true;
135 if( lgKeep ) lst.push_back( sstr );
140 if( lgStrict ) lgFail =
true;
141 if( lgKeep ) lst.push_back(
"" );
146 fprintf(
ioQQQ,
" A syntax error occurred while splitting the string: \"%s\"\n", str.c_str() );
147 fprintf(
ioQQQ,
" The separator is \"%s\". Empty substrings are not allowed.\n", sep.c_str() );
153void MyAssert(
const char *file,
int line,
const char *comment)
157 fprintf(
ioQQQ,
"\n\n\n PROBLEM DISASTER\n An assert has been thrown, this is bad.\n");
158 fprintf(
ioQQQ,
" %s\n",comment);
159 fprintf(
ioQQQ,
" It happened in the file %s at line number %i\n", file, line );
160 fprintf(
ioQQQ,
" This is iteration %li, nzone %li, fzone %.2f, lgSearch=%c.\n",
192 if( (
hextra.cryden == 0.) &&
h2 != NULL &&
h2->xFracLim > 0.1 )
194 fprintf(
ioQQQ,
" >>> \n >>> \n >>> Cosmic rays are not included and the gas is molecular. "
195 "THIS IS KNOWN TO BE UNSTABLE. Add cosmic rays and try again.\n >>> \n >>>\n\n");
199 fprintf(
ioQQQ,
"\n\n\n" );
200 fprintf(
ioQQQ,
" vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv \n" );
201 fprintf(
ioQQQ,
" > PROBLEM DISASTER PROBLEM DISASTER. <\n" );
202 fprintf(
ioQQQ,
" > Sorry, something bad has happened. <\n" );
203 fprintf(
ioQQQ,
" > Please post this on the Cloudy web site <\n" );
204 fprintf(
ioQQQ,
" > discussion board at www.nublado.org <\n" );
205 fprintf(
ioQQQ,
" > Please send all following information: <\n" );
206 fprintf(
ioQQQ,
" ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ \n" );
207 fprintf(
ioQQQ,
"\n\n" );
210 fprintf(
ioQQQ,
" Cloudy version number is %s\n",
214 fprintf(
ioQQQ,
"%5ld warnings,%3ld cautions,%3ld temperature failures. Messages follow.\n",
230 fprintf(
ioQQQ,
" This input stream included an init file.\n");
231 fprintf(
ioQQQ,
" If this init file is not part of the standard Cloudy distribution\n");
232 fprintf(
ioQQQ,
" then I will need a copy of it too.\n");
251 for( i=0; i < 4; i++ )
254 chCAP[i] =
toupper( chLab[i] );
271 while( chCard[i]!=
'\0' )
273 chCard[i] =
tolower( chCard[i] );
288 while( chCard[i]!=
'\0' )
290 chCard[i] =
toupper( chCard[i] );
305 double hold =
sexp(t) - t*
ee1(t);
308 return max( hold, 0. );
317 static double a[6]={-.57721566,.99999193,-.24991055,.05519968,-.00976004,
319 static double b[4]={8.5733287401,18.0590169730,8.6347608925,.2677737343};
320 static double c[4]={9.5733223454,25.6329561486,21.0996530827,3.9584969228};
333 fprintf(
ioQQQ,
" DISASTER negative argument in function ee1, x<0\n" );
341 ans = ((((a[5]*x + a[4])*x + a[3])*x + a[2])*x + a[1])*x + a[0] - log(x);
348 top = (((x + b[0])*x + b[1])*x + b[2])*x + b[3];
349 bot = (((x + c[0])*x + c[1])*x + c[2])*x + c[3];
350 ans = top/bot/x*exp(-x);
363 static double b[4]={8.5733287401,18.0590169730,8.6347608925,.2677737343};
364 static double c[4]={9.5733223454,25.6329561486,21.0996530827,3.9584969228};
372 top = (((x + b[0])*x + b[1])*x + b[2])*x + b[3];
374 bot = (((x + c[0])*x + c[1])*x + c[2])*x + c[3];
390 const char *eol_ptr = &chCard[last];
391 const char *ptr =
min(&chCard[*ipnt-1],eol_ptr);
393 ASSERT( *ipnt > 0 && *ipnt < last );
395 while( ptr < eol_ptr && ( chr = *ptr++ ) !=
'\0' )
397 const char *lptr = ptr;
399 if( lchr ==
'-' || lchr ==
'+' )
407 if( ptr >= eol_ptr || chr ==
'\0' )
415 bool lgCommaFound =
false, lgLastComma =
false;
418 lgCommaFound = lgLastComma;
434 while( isdigit(chr) || chr ==
'.' || chr ==
'-' || chr ==
'+' || chr ==
',' || chr ==
'e' || chr ==
'E' );
438 fprintf(
ioQQQ,
" PROBLEM - a comma was found embedded in a number, this is deprecated.\n" );
439 fprintf(
ioQQQ,
"== %-80s ==\n",chCard);
442 double value = atof(chNumber.c_str());
444 *ipnt = (long)(ptr - chCard);
459 ASSERT( strlen(chKey) > 0 );
461 if( ( ptr =
strstr_s( chCard, chKey ) ) == NULL )
469 Match_v = (long)(ptr-chCard+1);
491 fudgec.lgFudgeUsed =
true;
493 else if( ipnt >=
fudgec.nfudge )
495 fprintf(
ioQQQ,
" FUDGE factor not entered for array number %3ld\n",
501 fudge_v =
fudgec.fudgea[ipnt];
502 fudgec.lgFudgeUsed =
true;
545 if( i0 == NULL || i1 == NULL )
551 " A filename or label must be specified within double quotes, but no quotes were encountered on this command.\n" );
552 fprintf(
ioQQQ,
" Name must be surrounded by exactly two double quotes, like \"name.txt\". \n" );
553 fprintf(
ioQQQ,
" The line image follows:\n" );
554 fprintf(
ioQQQ,
" %s\n", chCardRaw);
555 fprintf(
ioQQQ,
" Sorry\n" );
561 chStringOut[0] =
'\0';
569 strncpy(chStringOut,i0+1,len);
571 chStringOut[len] =
'\0';
604double powi(
double x,
long int n )
647 if( m == 0 || (n < 0 && m > 1) )
690char *
PrintEfmt(
const char *fmt,
double val )
697 sprintf(buf, fmt, val);
714 if((ep =
strchr_s(buf,
'e')) == NULL)
741 double frac , xlog , xfloor , tvalue;
748 fprintf(ioOUT,
"********");
750 else if( value <= DBL_MIN )
752 fprintf(ioOUT,
"0.00E+00");
758 xlog = log10( tvalue );
759 xfloor = floor(xlog);
762 frac = tvalue*pow(10.,-xfloor);
764 frac = (10.*tvalue)*pow(10.,-(xfloor+1.));
773 fprintf(ioOUT,
"%.2f",frac);
781 fprintf(ioOUT,
"%.2d",iExp);
790 double frac , xlog , xfloor , tvalue;
797 fprintf(ioOUT,
"*******");
799 else if( value <= DBL_MIN )
801 fprintf(ioOUT,
"0.0E+00");
807 xlog = log10( tvalue );
808 xfloor = floor(xlog);
811 frac = tvalue*pow(10.,-xfloor);
813 frac = (10.*tvalue)*pow(10.,-(xfloor+1.));
822 fprintf(ioOUT,
"%.1f",frac);
830 fprintf(ioOUT,
"%.2d",iExp);
840 double frac , xlog , xfloor, tvalue;
847 fprintf(ioOUT,
"*********");
849 else if( value <= DBL_MIN )
851 fprintf(ioOUT,
"0.000E+00");
857 xlog = log10( tvalue );
858 xfloor = floor(xlog);
861 frac = tvalue*pow(10.,-xfloor);
863 frac = (10.*tvalue)*pow(10.,-(xfloor+1.));
872 fprintf(ioOUT,
"%5.3f",frac);
880 fprintf(ioOUT,
"%.2d",iExp);
893 fprintf(
ioQQQ,
" Something that cannot happen, has happened.\n" );
894 fprintf(
ioQQQ,
" This is TotalInsanity, I live in %s.\n", __FILE__ );
906 fprintf(
ioQQQ,
" A read of internal input data has failed.\n" );
907 fprintf(
ioQQQ,
" This is BadRead, I live in %s.\n", __FILE__ );
986 broke.lgBroke =
true;
995 broke.lgFixit =
true;
1004 broke.lgCheckit =
true;
1015 va_start(ap,format);
1016 i1 = fprintf(fp,
"DEBUG ");
1018 i2 = vfprintf(fp,format,ap);
1039 if(debug <=
trace.debug_level)
1043 i = vfprintf(
ioQQQ, fmt, ap);
1057 double (*fct)(
double) )
1059 double a = 0.5*(xu+xl),
1078 double weights[16] = {
1079 .35093050047350483e-2, .81371973654528350e-2, .12696032654631030e-1, .17136931456510717e-1,
1080 .21417949011113340e-1, .25499029631188088e-1, .29342046739267774e-1, .32911111388180923e-1,
1081 .36172897054424253e-1, .39096947893535153e-1, .41655962113473378e-1, .43826046502201906e-1,
1082 .45586939347881942e-1, .46922199540402283e-1, .47819360039637430e-1, .48270044257363900e-1};
1085 .498631930924740780, .49280575577263417, .4823811277937532200, .46745303796886984000,
1086 .448160577883026060, .42468380686628499, .3972418979839712000, .36609105937014484000,
1087 .331522133465107600, .29385787862038116, .2534499544661147000, .21067563806531767000,
1088 .165934301141063820, .11964368112606854, .7223598079139825e-1, .24153832843869158e-1};
1090 for(
int i=0; i<16; i++)
1092 y += b * weights[i] * ((*fct)(a+b*c[i]) + (*fct)(a-b*c[i]));
1221 if( kk != 1 && kk != 2 )
1231 for( i=0; i < nn; i++ )
1246 for( i=0; i < nn; i++ )
1264 if( r <= 0.5898437e0 )
1278 ij = i + (long)((j-i)*r);
1283 if( x[iperm[i-1]-1] > x[lm-1] )
1285 iperm[ij-1] = iperm[i-1];
1293 if( x[iperm[j-1]-1] < x[lm-1] )
1295 iperm[ij-1] = iperm[j-1];
1302 if( x[iperm[i-1]-1] > x[lm-1] )
1304 iperm[ij-1] = iperm[i-1];
1315 if( x[iperm[l-1]-1] <= x[lm-1] )
1323 if( x[iperm[k-1]-1] >= x[lm-1] )
1331 iperm[l-1] = iperm[k-1];
1365 if( x[iperm[i-1]-1] > x[lm-1] )
1371 iperm[k] = iperm[k-1];
1374 if( x[lm-1] >= x[iperm[k-1]-1] )
1396 for( i=0; i < nn; i++ )
1408 for( istrt=1; istrt <= nn; istrt++ )
1411 if( iperm[istrt_] >= 0 )
1416 while( iperm[indx-1] > 0 )
1418 x[indx-1] = x[iperm[indx-1]-1];
1420 iperm[indx-1] = -iperm[indx-1];
1421 indx = labs(iperm[indx-1]);
1428 for( i=0; i < nn; i++ )
1430 iperm[i] = -iperm[i];
1434 for( i=0; i < nn; i++ )
1463 enum{DEBUG_LOC=
false};
1466 static long int kount=0, nTot=0;
1468 fprintf(
ioQQQ,
"%li\t%li\t%li\n",
1476 if( ( ptr = malloc( size ) ) == NULL )
1478 fprintf(
ioQQQ,
"DISASTER MyMalloc could not allocate %lu bytes. Exit in MyMalloc.",
1479 (
unsigned long)size );
1480 fprintf(
ioQQQ,
"MyMalloc called from file %s at line %i.\n",
1483 if(
struc.nzlim>2000 )
1484 fprintf(
ioQQQ,
"This may have been caused by the large number of zones."
1485 " %li zones were requested. Is this many zones really necessary?\n",
1492# if !defined(NDEBUG) && !defined(NOINIT)
1494 size_t nFloat = size/4;
1495 size_t nDouble = size/8;
1497 double *dptr =
static_cast<double*
>(ptr);
1503 if( size == nDouble*8 )
1512 set_NaN( dptr, (
long)nDouble );
1514 else if( size == nFloat*4 )
1517 set_NaN( fptr, (
long)nFloat );
1521 memset( ptr, 0xff, size );
1546 enum{DEBUG_LOC=
false};
1549 static long int kount=0;
1550 fprintf(
ioQQQ,
"%li\tcall\t%li\tbytes\n", kount,
1556 if( ( ptr = calloc( num , size ) ) == NULL )
1558 fprintf(
ioQQQ,
"MyCalloc could not allocate %lu bytes. Exit in MyCalloc.",
1559 (
unsigned long)size );
1582 enum{DEBUG_LOC=
false};
1585 static long int kount=0;
1586 fprintf(
ioQQQ,
"%li\tcall\t%li\tbytes\n", kount,
1592 if( ( ptr = realloc( p , size ) ) == NULL )
1594 fprintf(
ioQQQ,
"MyRealloc could not allocate %lu bytes. Exit in MyRealloc.",
1595 (
unsigned long)size );
1615 csphot_v =
opac.OpacStack[inu-ithr+iofset-1];
1649 double x1,
x2, w, yy1;
1651 static int use_last =
false;
1666 }
while ( w >= 1.0 );
1668 w = sqrt((-2.0*log(w))/w);
1673 return xMean + yy1 * s;
1690 ASSERT( PctUncertainty < 0.5 );
1692 StdDev = PctUncertainty;
1700 while( (result < 1.-3.*PctUncertainty) || (result > 1.+3.*PctUncertainty) );
1702 ASSERT( result>0. && result<2. );
1716 if(
rfield.ContBoltz[ip] <= 0. )
1735 if( line[0] ==
'#' )
1738 size_t p = line.find(
"XXXX" );
1739 if( p != string::npos )
1741 fprintf(
ioQQQ,
"%s\n", line.c_str() );
1752 if( line[0] ==
'#' )
1755 size_t p = line.find(
"XXXX" );
1756 if( p != string::npos )
1757 line.replace( p, 4,
atmdat.chVersion );
1758 fprintf(
ioQQQ,
"%s\n", line.c_str() );
1780 istream::sentry se(is,
true);
1781 streambuf* sb = is.rdbuf();
1785 int c = sb->sbumpc();
1789 if( sb->sgetc() == EOF )
1790 is.setstate(ios::eofbit);
1793 if( sb->sgetc() ==
'\n' )
1795 if( sb->sgetc() == EOF )
1796 is.setstate(ios::eofbit);
1800 is.setstate(ios::eofbit);
const char * strstr_s(const char *haystack, const char *needle)
const char * strchr_s(const char *s, int c)
#define DEBUG_ENTRY(funcname)
void cdCautions(FILE *ioOUT)
void cdPrintCommands(FILE *ioOUT)
void cdWarnings(FILE *ioPNT)
double get(const char *unit) const
static t_version & Inst()
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
void set_NaN(sys_float &x)
const ios_base::openmode mode_r
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
molezone * findspecieslocal(const char buf[])
UNUSED const double FR1RYD
double csphot(long int inu, long int ithr, long int iofset)
double plankf(long int ip)
sys_float sexp(sys_float x)
double MyGaussRand(double PctUncertainty)
long ipow(long m, long n)
long nMatch(const char *chKey, const char *chCard)
double ee1_safe(double x)
void * MyCalloc(size_t num, size_t size)
double powi(double x, long int n)
double qg32(double xl, double xu, double(*fct)(double))
void PrintE82(FILE *ioOUT, double value)
NORETURN void BadRead(void)
double fudge(long int ipnt)
int dbg_printf(int debug, const char *fmt,...)
int dprintf(FILE *fp, const char *format,...)
double RandGauss(double xMean, double s)
void cap4(char *chCAP, const char *chLab)
void MyAssert(const char *file, int line, const char *comment)
int GetQuote(char *chStringOut, char *chCard, char *chCardRaw, bool lgAbort)
void Split(const string &str, const string &sep, vector< string > &lst, split_mode mode)
void PrintE71(FILE *ioOUT, double value)
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
void * MyRealloc(void *p, size_t size)
void * MyMalloc(size_t size, const char *chFile, int line)
void CloudyPrintReference()
void PrintE93(FILE *ioOUT, double value)
void DatabasePrintReference()
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
void spsort(realnum x[], long int n, long int iperm[], int kflag, int *ier)
NORETURN void TotalInsanity(void)
double AnuUnit(realnum energy_ryd)
void uncaps(char *chCard)
istream & SafeGetline(istream &is, string &t)