cloudy trunk
Loading...
Searching...
No Matches
service.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*
4* a set of routines that are widely used across the code for various
5* housekeeping chores. These do not do any physics and are unlikely to
6* change over time. The prototypes are in cddefines.h and so are
7* automatically picked up by all routines
8*/
9/*FFmtRead scan input line for free format number */
10/*e2 second exponential integral */
11/*caps convert input command line (through eol) to ALL CAPS */
12/*ShowMe produce request to send information to GJF after a crash */
13/*AnuUnit produce continuum energy in arbitrary units */
14/*cap4 convert first 4 char of input line chLab into chCAP all in caps, null termination */
15/*insane set flag saying that insanity has occurred */
16/*nMatch determine whether match to a keyword occurs on command line,
17 * return value is 0 if no match, and position of match within string if hit */
18/*fudge enter fudge factors, or some arbitrary number, with fudge command*/
19/*GetQuote get any name between double quotes off command line
20 * return string as chLabel, is null terminated */
21/*qip compute pow(x,n) for positive integer n through repeated squares */
22/*dsexp safe exponential function for doubles */
23/*sexp safe exponential function */
24/*TestCode set flag saying that test code is in place */
25/*CodeReview - placed next to code that needs to be checked */
26/*fixit - say that code needs to be fixed */
27/*broken set flag saying that the code is broken, */
28/*dbg_printf is a debug print routine that was provided by Peter Teuben,
29 * as a component from his NEMO package. It offers run-time specification
30 * of the level of debugging */
31/*qg32 32 point Gaussian quadrature, original Fortran given to Gary F by Jim Lattimer */
32/*TotalInsanity general error handler for something that cannot happen */
33/*BadRead general error handler for trying to read data, but failing */
34/*MyMalloc wrapper for malloc(). Returns a good pointer or dies. */
35/*MyCalloc wrapper for calloc(). Returns a good pointer or dies. */
36/*spsort netlib routine to sort array returning sorted indices */
37/*chLineLbl use information in line transfer arrays to generate a line label *
38 * this label is null terminated */
39/*chIonLbl use information in line array to generate a null terminated ion label in "Fe 2" */
40/*csphot returns photoionization cross section from opacity stage using std pointers */
41/*MyAssert a version of assert that fails gracefully */
42/*RandGauss normal random variate generator */
43/*MyGaussRand a wrapper for RandGauss, see below */
44
45#include "cdstd.h"
46#include <cstdarg> /* ANSI variable arg macros */
47#include "cddefines.h"
48#include "physconst.h"
49#include "cddrive.h"
50#include "called.h"
51#include "opacity.h"
52#include "rfield.h"
53#include "hextra.h"
54#include "struc.h"
55#include "hmi.h"
56#include "fudgec.h"
57#include "broke.h"
58#include "trace.h"
59#include "input.h"
60#include "save.h"
61#include "version.h"
62#include "warnings.h"
63#include "conv.h"
64#include "thirdparty.h"
65#include "mole.h"
66#include "atmdat.h"
67
68/*read_whole_line - safe version of fgets - read a line,
69 * return null if cannot read line or if input line is too long */
70char *read_whole_line( char *chLine , int nChar , FILE *ioIN )
71{
72 char *chRet;
73
74 DEBUG_ENTRY( "read_whole_line()" );
75
76 /* wipe the buffer to prevent the code from accidentally picking up on previous input */
77 memset( chLine, 0, (size_t)nChar );
78
79 /* this always writes a '\0' character, even if line does not fit in buffer
80 * the terminating newline is copied only if the line does fit in the buffer */
81 if( (chRet = fgets( chLine, nChar, ioIN )) == NULL )
82 return NULL;
83
84 long lineLength = strlen( chRet );
85 //fprintf(ioQQQ , "DEBUG reading:%s\n" , chLine);
86 //fprintf(ioQQQ , "DEBUG length is %li nChar is %i \n", lineLength , nChar);
87
88 /* return null if input string is longer than nChar-1 (including terminating newline),
89 * the longest we can read. Print and return null but chLine still has as much of
90 * the line as could be placed in the buffer */
91 if( lineLength>=nChar-1 )
92 {
93 if( called.lgTalk )
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",
97 nChar , chLine );
98
99 lgAbort = true;
100 return NULL;
101 }
102 return chRet;
103}
104
106void Split(const string& str, // input string
107 const string& sep, // separator, may be multiple characters
108 vector<string>& lst, // the separated items will be appended here
109 split_mode mode) // SPM_RELAX, SPM_KEEP_EMPTY, or SPM_STRICT; see cddefines.h
110{
111 DEBUG_ENTRY( "Split()" );
112
113 bool lgStrict = ( mode == SPM_STRICT );
114 bool lgKeep = ( mode == SPM_KEEP_EMPTY );
115 bool lgFail = false;
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 );
121 else {
122 if( lgStrict ) lgFail = true;
123 if( lgKeep ) lst.push_back( sstr );
124 }
125 while( ptr2 != string::npos ) {
126 // the separator is skipped
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 );
133 else {
134 if( lgStrict ) lgFail = true;
135 if( lgKeep ) lst.push_back( sstr );
136 }
137 }
138 else {
139 ptr2 = string::npos;
140 if( lgStrict ) lgFail = true;
141 if( lgKeep ) lst.push_back( "" );
142 }
143 }
144 if( lgFail )
145 {
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() );
149 }
150}
151
152/* a version of assert that fails gracefully */
153void MyAssert(const char *file, int line, const char *comment)
154{
155 DEBUG_ENTRY( "MyAssert()" );
156
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",
161 iteration ,
162 nzone ,
163 fnzone ,
164 TorF(conv.lgSearch) );
165
166 ShowMe();
167# ifdef OLD_ASSERT
169# endif
170}
171
172/*AnuUnit produce continuum energy in arbitrary units, as determined by ChkUnits() */
173double AnuUnit(realnum energy_ryd)
174{
175 DEBUG_ENTRY( "AnuUnit()" );
176
177 return Energy((double)energy_ryd).get(save.chConPunEnr[save.ipConPun]);
178}
179
180/*ShowMe produce request to send information to GJF after a crash */
181void ShowMe(void)
182{
183
184 DEBUG_ENTRY( "ShowMe()" );
185
186 /* print info if output unit is defined */
187 if( ioQQQ != NULL )
188 {
189 /* >>chng 06 mar 02 - check if molecular but cosmic rays are ignored */
191 // molecular species may not be set up yet, so check for NULL pointer...
192 if( (hextra.cryden == 0.) && h2 != NULL && h2->xFracLim > 0.1 )
193 {
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");
196 }
197 else
198 {
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" );
208
209
210 fprintf( ioQQQ, " Cloudy version number is %s\n",
211 t_version::Inst().chVersion );
212 fprintf( ioQQQ, " %s\n\n", t_version::Inst().chInfo );
213
214 fprintf( ioQQQ, "%5ld warnings,%3ld cautions,%3ld temperature failures. Messages follow.\n",
215 warnings.nwarn, warnings.ncaun, conv.nTeFail );
216
217 /* print the warnings first */
219
220 /* now print the cautions */
222
223 /* now output the commands */
225
226 /* if init command was present, this is the number of lines in it -
227 * if no init then still set to zero as done in cdInit */
228 if( input.nSaveIni )
229 {
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");
233 }
234 }
235 }
236 return;
237}
238
239/*cap4 convert first 4 char of input line chLab into chCAP all in caps, null termination */
240void cap4(
241 char *chCAP , /* output string, cap'd first 4 char of chLab, */
242 /* with null terminating */
243 const char *chLab) /* input string ending with eol*/
244{
245 long int /*chr,*/
246 i;
247
248 DEBUG_ENTRY( "cap4()" );
249
250 /* convert 4 character string in chLab to ALL CAPS in chCAP */
251 for( i=0; i < 4; i++ )
252 {
253 /* toupper is function in ctype that converts to upper case */
254 chCAP[i] = toupper( chLab[i] );
255 }
256
257 /* now end string with eol */
258 chCAP[4] = '\0';
259 return;
260}
261
262/*uncaps convert input command line (through eol) to all lowercase */
263void uncaps(char *chCard )
264{
265 long int i;
266
267 DEBUG_ENTRY( "caps()" );
268
269 /* convert full character string in chCard to ALL CAPS */
270 i = 0;
271 while( chCard[i]!= '\0' )
272 {
273 chCard[i] = tolower( chCard[i] );
274 ++i;
275 }
276 return;
277}
278
279/*caps convert input command line (through eol) to ALL CAPS */
280void caps(char *chCard )
281{
282 long int i;
283
284 DEBUG_ENTRY( "caps()" );
285
286 /* convert full character string in chCard to ALL CAPS */
287 i = 0;
288 while( chCard[i]!= '\0' )
289 {
290 chCard[i] = toupper( chCard[i] );
291 ++i;
292 }
293 return;
294}
295
296/*e2 second exponential integral */
297/*>>chng 07 jan 17, PvH discover that exp-t is not really
298 * exp-t - this changed results in several tests */
299double e2(
300 /* the argument to E2 */
301 double t )
302{
303 /* use recurrence relation */
304 /* ignore exp_mt, it is *very* unreliable */
305 double hold = sexp(t) - t*ee1(t);
306 DEBUG_ENTRY( "e2()" );
307 /* guard against negative results, this can happen for very large t */
308 return max( hold, 0. );
309}
310
311/*ee1 first exponential integral */
312double ee1(double x)
313{
314 double ans,
315 bot,
316 top;
317 static double a[6]={-.57721566,.99999193,-.24991055,.05519968,-.00976004,
318 .00107857};
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};
321
322 DEBUG_ENTRY( "ee1()" );
323
324 /* computes the exponential integral E1(x),
325 * from Abramowitz and Stegun
326 * stops with error condition for negative argument,
327 * returns zero in large x limit
328 * */
329
330 /* error - does not accept negative arguments */
331 if( x <= 0 )
332 {
333 fprintf( ioQQQ, " DISASTER negative argument in function ee1, x<0\n" );
335 }
336
337 /* branch for x less than 1 */
338 else if( x < 1. )
339 {
340 /* abs. accuracy better than 2e-7 */
341 ans = ((((a[5]*x + a[4])*x + a[3])*x + a[2])*x + a[1])*x + a[0] - log(x);
342 }
343
344 /* branch for x greater than, or equal to, one */
345 else
346 {
347 /* abs. accuracy better than 2e-8 */
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);
351 }
352 return ans;
353}
354
355/* same as ee1, except without factor of exp(x) in returned value */
356double ee1_safe(double x)
357{
358 double ans,
359 bot,
360 top;
361 /*static double a[6]={-.57721566,.99999193,-.24991055,.05519968,-.00976004,
362 .00107857};*/
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};
365
366 DEBUG_ENTRY( "ee1_safe()" );
367
368 ASSERT( x > 1. );
369
370 /* abs. accuracy better than 2e-8 */
371 /* top = powi(x,4) + b[0]*powi(x,3) + b[1]*x*x + b[2]*x + b[3]; */
372 top = (((x + b[0])*x + b[1])*x + b[2])*x + b[3];
373 /* bot = powi(x,4) + c[0]*powi(x,3) + c[1]*x*x + c[2]*x + c[3]; */
374 bot = (((x + c[0])*x + c[1])*x + c[2])*x + c[3];
375
376 ans = top/bot/x;
377 return ans;
378}
379
380/*FFmtRead scan input line for free format number */
381double FFmtRead(const char *chCard,
382 long int *ipnt,
383 /* the contents of this array element is the last that will be read */
384 long int last,
385 bool *lgEOL)
386{
387 DEBUG_ENTRY( "FFmtRead()" );
388
389 char chr = '\0';
390 const char *eol_ptr = &chCard[last]; // eol_ptr points one beyond last valid char
391 const char *ptr = min(&chCard[*ipnt-1],eol_ptr); // ipnt is on fortran scale
392
393 ASSERT( *ipnt > 0 && *ipnt < last );
394
395 while( ptr < eol_ptr && ( chr = *ptr++ ) != '\0' )
396 {
397 const char *lptr = ptr;
398 char lchr = chr;
399 if( lchr == '-' || lchr == '+' )
400 lchr = *lptr++;
401 if( lchr == '.' )
402 lchr = *lptr;
403 if( isdigit(lchr) )
404 break;
405 }
406
407 if( ptr >= eol_ptr || chr == '\0' )
408 {
409 *ipnt = last+1;
410 *lgEOL = true;
411 return 0.;
412 }
413
414 string chNumber;
415 bool lgCommaFound = false, lgLastComma = false;
416 do
417 {
418 lgCommaFound = lgLastComma;
419 if( chr != ',' )
420 {
421 chNumber += chr;
422 }
423 else
424 {
425 /* don't complain about comma if it appears after number,
426 as determined by exiting loop before this sets lgCommaFound */
427 lgLastComma = true;
428
429 }
430 if( ptr == eol_ptr )
431 break;
432 chr = *ptr++;
433 }
434 while( isdigit(chr) || chr == '.' || chr == '-' || chr == '+' || chr == ',' || chr == 'e' || chr == 'E' );
435
436 if( lgCommaFound )
437 {
438 fprintf( ioQQQ, " PROBLEM - a comma was found embedded in a number, this is deprecated.\n" );
439 fprintf(ioQQQ, "== %-80s ==\n",chCard);
440 }
441
442 double value = atof(chNumber.c_str());
443
444 *ipnt = (long)(ptr - chCard); // ptr already points 1 beyond where next read should start
445 *lgEOL = false;
446 return value;
447}
448
449/*nMatch determine whether match to a keyword occurs on command line,
450 * return value is 0 if no match, and position of match within string if hit */
451long nMatch(const char *chKey,
452 const char *chCard)
453{
454 const char *ptr;
455 long Match_v;
456
457 DEBUG_ENTRY( "nMatch()" );
458
459 ASSERT( strlen(chKey) > 0 );
460
461 if( ( ptr = strstr_s( chCard, chKey ) ) == NULL )
462 {
463 /* did not find match, return 0 */
464 Match_v = 0L;
465 }
466 else
467 {
468 /* return position within chCard (fortran scale) */
469 Match_v = (long)(ptr-chCard+1);
470 }
471 return Match_v;
472}
473
474/* fudge enter fudge factors, or some arbitrary number, with fudge command
475 * other sections of the code access these numbers by calling fudge
476 * fudge(0) returns the first number that was entered
477 * prototype for this function is in cddefines.h so that it can be used without
478 * declarations
479 * fudge(-1) queries the routine for the number of fudge parameters that were entered,
480 * zero returned if none */
481double fudge(long int ipnt)
482{
483 double fudge_v;
484
485 DEBUG_ENTRY( "fudge()" );
486
487 if( ipnt < 0 )
488 {
489 /* this is special case, return number of arguments */
490 fudge_v = fudgec.nfudge;
491 fudgec.lgFudgeUsed = true;
492 }
493 else if( ipnt >= fudgec.nfudge )
494 {
495 fprintf( ioQQQ, " FUDGE factor not entered for array number %3ld\n",
496 ipnt );
498 }
499 else
500 {
501 fudge_v = fudgec.fudgea[ipnt];
502 fudgec.lgFudgeUsed = true;
503 }
504 return fudge_v;
505}
506
507
508/* GetQuote get any name between double quotes off command line
509 * sets chLabel - null terminated string as chLabel
510 * sets string between quotes to spaces
511 * returns zero for success, 1 for did not find double quotes
512 * lgAbort - true then above, false only set flag*/
514 /* we will generate a label and stuff it here */
515 char *chStringOut,
516 /* line image, we will blank out anything between quotes */
517 char *chCard,
518 char *chCardRaw,
519 /* if true then abort if no double quotes found,
520 * if false then return null string in this case */
521 bool lgAbort )
522{
523 char *i0, /* pointer to first " */
524 *i1, /* pointer to second ", name is in between */
525 *iLoc; /* pointer to first " in local version of card in calling routine */
526 size_t len;
527
528 DEBUG_ENTRY( "GetQuote()" );
529
530 /* find first quote start of string, string begins and ends with quotes */
531 i0 = strchr_s( chCardRaw,'\"' );
532
533 if( i0 != NULL )
534 {
535 /* get pointer to next quote */
536 i1 = strchr_s( i0+1 ,'\"' );
537 }
538 else
539 {
540 i1 = NULL;
541 }
542
543 /* check that pointers were not NULL */
544 /* >>chng 00 apr 27, check for i0 and i1 equal not needed anymore, by PvH */
545 if( i0 == NULL || i1 == NULL )
546 {
547 if( lgAbort )
548 {
549 /* lgAbort true, must abort if no string found */
550 fprintf( ioQQQ,
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" );
557 }
558 else
559 {
560 /* this branch, ok if not present, return null string in that case */
561 chStringOut[0] = '\0';
562 /* return value of 1 indicates did not find double quotes */
563 return 1;
564 }
565 }
566
567 /* now copy the text in between quotes */
568 len = (size_t)(i1-i0-1);
569 strncpy(chStringOut,i0+1,len);
570 /* strncpy doesn't terminate the label */
571 chStringOut[len] = '\0';
572
573 /* get pointer to first quote in local copy of line image in calling routine */
574 iLoc = strchr_s( chCard, '\"' );
575 if( iLoc == NULL )
577
578 // blank out label once finished, to not be picked up later
579 // erase quotes as well, so that we can find second label, by PvH
580 while( i0 <= i1 )
581 {
582 *i0 = ' ';
583 *iLoc = ' ';
584 ++i0;
585 ++iLoc;
586 }
587 /* return condition of 0 indicates success */
588 return 0;
589}
590
591/* want to define this only if no native os support exists */
592#ifndef HAVE_POWI
593
594/* powi.c - calc x^n, where n is an integer! */
595
596/* Very slightly modified version of power() from Computer Language, Sept. 86,
597 pg 87, by Jon Snader (who extracted the binary algorithm from Donald Knuth,
598 "The Art of Computer Programming", vol 2, 1969).
599 powi() will only be called when an exponentiation with an integer
600 exponent is performed, thus tests & code for fractional exponents were
601 removed.
602 */
603
604double powi( double x, long int n ) /* returns: x^n */
605/* x; base */
606/* n; exponent */
607{
608 double p; /* holds partial product */
609
610 DEBUG_ENTRY( "powi()" );
611
612 if( x == 0 )
613 return 0.;
614
615 /* test for negative exponent */
616 if( n < 0 )
617 {
618 n = -n;
619 x = 1/x;
620 }
621
622 p = is_odd(n) ? x : 1; /* test & set zero power */
623
624 /*lint -e704 shift right of signed quantity */
625 /*lint -e720 Boolean test of assignment */
626 while( n >>= 1 )
627 { /* now do the other powers */
628 x *= x; /* sq previous power of x */
629 if( is_odd(n) ) /* if low order bit set */
630 p *= x; /* then, multiply partial product by latest power of x */
631 }
632 /*lint +e704 shift right of signed quantity */
633 /*lint +e720 Boolean test of assignment */
634 return p;
635}
636
637#endif /* HAVE_POWI */
638
639long ipow( long m, long n ) /* returns: m^n */
640/* m; base */
641/* n; exponent */
642{
643 long p; /* holds partial product */
644
645 DEBUG_ENTRY( "ipow()" );
646
647 if( m == 0 || (n < 0 && m > 1) )
648 return 0L;
649 /* NOTE: negative exponent always results in 0 for integers!
650 * (except for the case when m==1 or -1) */
651
652 if( n < 0 )
653 { /* test for negative exponent */
654 n = -n;
655 m = 1/m;
656 }
657
658 p = is_odd(n) ? m : 1; /* test & set zero power */
659
660 /*lint -e704 shift right of signed quantity */
661 /*lint -e720 Boolean test of assignment */
662 while( n >>= 1 )
663 { /* now do the other powers */
664 m *= m; /* sq previous power of m */
665 if( is_odd(n) ) /* if low order bit set */
666 p *= m; /* then, multiply partial product by latest power of m */
667 }
668 /*lint +e704 shift right of signed quantity */
669 /*lint +e720 Boolean test of assignment */
670 return p;
671}
672
673/*PrintE82 - series of routines to mimic 1p, e8.2 fortran output */
674/***********************************************************
675 * contains the following sets of routines to get around *
676 * the MS C++ compilers unusual exponential output. *
677 * PrintEfmt <= much faster, no overhead with unix *
678 * PrintE93 *
679 * PrintE82 *
680 * PrintE71 *
681 **********************************************************/
682
683#ifdef _MSC_VER
684/**********************************************************/
685/*
686 * Instead of printf'ing with %e or %.5e or whatever, call
687 * efmt("%e", val) and print the result with %s. This lets
688 * us work around bugs in VS C 6.0.
689 */
690char *PrintEfmt(const char *fmt, double val /* , char *buf */)
691{
692 static char buf[30]; /* or pass it in */
693
694 DEBUG_ENTRY( "PrintEfmt()" );
695
696 /* create the string */
697 sprintf(buf, fmt, val);
698
699 /* code to fix incorrect ms v e format. works only for case where there is
700 * a leading space in the format - for formats like 7.1, 8.2, 9.3, 10.4, etc
701 * result will have 1 too many characters */
702 char *ep , buf2[30];
703
704 /* msvc behaves badly in different ways for positive vs negative sign vals,
705 * if val is positive must create a leading space */
706 if( val >= 0.)
707 {
708 strcpy(buf2 , " " );
709 strcat(buf2 , buf);
710 strcpy( buf , buf2);
711 }
712
713 /* allow for both e and E formats */
714 if((ep = strchr_s(buf, 'e')) == NULL)
715 {
716 ep = strchr_s(buf, 'E');
717 }
718
719 /* ep can still be NULL if val is Inf or NaN */
720 if(ep != NULL)
721 {
722 /* move pointer two char past the e, to pick up the e and sign */
723 ep += 2;
724
725 /* terminate buf where the e is, *ep points to this location */
726 *ep = '\0';
727
728 /* skip next char, */
729 ++ep;
730
731 /* copy resulting string to return string */
732 strcat( buf, ep );
733 }
734 return buf;
735}
736#endif
737
738/**********************************************************/
739void PrintE82( FILE* ioOUT, double value )
740{
741 double frac , xlog , xfloor , tvalue;
742 int iExp;
743
744 DEBUG_ENTRY( "PrintE82()" );
745
746 if( value < 0 )
747 {
748 fprintf(ioOUT,"********");
749 }
750 else if( value <= DBL_MIN )
751 {
752 fprintf(ioOUT,"0.00E+00");
753 }
754 else
755 {
756 /* round number off for 8.2 format (not needed since can't be negative) */
757 tvalue = value;
758 xlog = log10( tvalue );
759 xfloor = floor(xlog);
760 /* this is now the fractional part */
761 if (xfloor < 0.)
762 frac = tvalue*pow(10.,-xfloor);
763 else
764 frac = (10.*tvalue)*pow(10.,-(xfloor+1.));
765 /*this is the possibly signed exponential part */
766 iExp = (int)xfloor;
767 if( frac>9.9945 )
768 {
769 frac /= 10.;
770 iExp += 1;
771 }
772 /* print the fractional part*/
773 fprintf(ioOUT,"%.2f",frac);
774 /* E for exponent */
775 fprintf(ioOUT,"E");
776 /* if positive throw a + sign*/
777 if(iExp>=0 )
778 {
779 fprintf(ioOUT,"+");
780 }
781 fprintf(ioOUT,"%.2d",iExp);
782 }
783 return;
784}
785/*
786 *==============================================================================
787 */
788void PrintE71( FILE* ioOUT, double value )
789{
790 double frac , xlog , xfloor , tvalue;
791 int iExp;
792
793 DEBUG_ENTRY( "PrintE71()" );
794
795 if( value < 0 )
796 {
797 fprintf(ioOUT,"*******");
798 }
799 else if( value <= DBL_MIN )
800 {
801 fprintf(ioOUT,"0.0E+00");
802 }
803 else
804 {
805 /* round number off for 8.2 format (not needed since can't be negative) */
806 tvalue = value;
807 xlog = log10( tvalue );
808 xfloor = floor(xlog);
809 /* this is now the fractional part */
810 if (xfloor < 0.)
811 frac = tvalue*pow(10.,-xfloor);
812 else
813 frac = (10.*tvalue)*pow(10.,-(xfloor+1.));
814 /*this is the possibly signed exponential part */
815 iExp = (int)xfloor;
816 if( frac>9.9945 )
817 {
818 frac /= 10.;
819 iExp += 1;
820 }
821 /* print the fractional part*/
822 fprintf(ioOUT,"%.1f",frac);
823 /* E for exponent */
824 fprintf(ioOUT,"E");
825 /* if positive throw a + sign*/
826 if(iExp>=0 )
827 {
828 fprintf(ioOUT,"+");
829 }
830 fprintf(ioOUT,"%.2d",iExp);
831 }
832 return;
833}
834
835/*
836 *==============================================================================
837 */
838void PrintE93( FILE* ioOUT, double value )
839{
840 double frac , xlog , xfloor, tvalue;
841 int iExp;
842
843 DEBUG_ENTRY( "PrintE93()" );
844
845 if( value < 0 )
846 {
847 fprintf(ioOUT,"*********");
848 }
849 else if( value <= DBL_MIN )
850 {
851 fprintf(ioOUT,"0.000E+00");
852 }
853 else
854 {
855 /* round number off for 9.3 format, neg numb not possible */
856 tvalue = value;
857 xlog = log10( tvalue );
858 xfloor = floor(xlog);
859 /* this is now the fractional part */
860 if (xfloor < 0.)
861 frac = tvalue*pow(10.,-xfloor);
862 else
863 frac = (10.*tvalue)*pow(10.,-(xfloor+1.));
864 /*this is the possibly signed exponential part */
865 iExp = (int)xfloor;
866 if( frac>9.99949 )
867 {
868 frac /= 10.;
869 iExp += 1;
870 }
871 /* print the fractional part*/
872 fprintf(ioOUT,"%5.3f",frac);
873 /* E for exponent */
874 fprintf(ioOUT,"E");
875 /* if positive throw a + sign*/
876 if(iExp>=0 )
877 {
878 fprintf(ioOUT,"+");
879 }
880 fprintf(ioOUT,"%.2d",iExp);
881 }
882 return;
883}
884
885/*TotalInsanity general error handler for something that cannot happen */
887{
888 DEBUG_ENTRY( "TotalInsanity()" );
889
890 /* something that cannot happen, happened,
891 * if this message is triggered, simply place a breakpoint here
892 * and debug the error */
893 fprintf( ioQQQ, " Something that cannot happen, has happened.\n" );
894 fprintf( ioQQQ, " This is TotalInsanity, I live in %s.\n", __FILE__ );
895 ShowMe();
896
898}
899
900/*BadRead general error handler for trying to read data, but failing */
902{
903 DEBUG_ENTRY( "BadRead()" );
904
905 /* read failed */
906 fprintf( ioQQQ, " A read of internal input data has failed.\n" );
907 fprintf( ioQQQ, " This is BadRead, I live in %s.\n", __FILE__ );
908 ShowMe();
909
911}
912
913/*sexp safe exponential function */
915{
916 sys_float sexp_v;
917
918 DEBUG_ENTRY( "sexp()" );
919
920 /* SEXP_LIMIT is 84 in cddefines.h */
921 if( x < SEXP_LIMIT )
922 {
923 sexp_v = exp(-x);
924 }
925 else
926 {
927 sexp_v = 0.f;
928 }
929 return sexp_v;
930}
931
932/*sexp safe exponential function */
933double sexp(double x)
934{
935 double sexp_v;
936
937 DEBUG_ENTRY( "sexp()" );
938
939 /* SEXP_LIMIT is 84 in cddefines.h */
940 if( x < SEXP_LIMIT )
941 {
942 sexp_v = exp(-x);
943 }
944 else
945 {
946 sexp_v = 0.;
947 }
948 return sexp_v;
949}
950
951
952/*dsexp safe exponential function for doubles */
953double dsexp(double x)
954{
955 double dsexp_v;
956
957 DEBUG_ENTRY( "dsexp()" );
958
959 if( x < DSEXP_LIMIT )
960 {
961 dsexp_v = exp(-x);
962 }
963 else
964 {
965 dsexp_v = 0.;
966 }
967 return dsexp_v;
968}
969
970/*TestCode set flag saying that test code is in place
971 * prototype in cddefines.h */
972void TestCode(void)
973{
974 DEBUG_ENTRY( "TestCode( )" );
975
976 /* called if test code is in place */
977 lgTestCodeCalled = true;
978 return;
979}
980
981/*broken set flag saying that the code is broken, */
982void broken(void)
983{
984 DEBUG_ENTRY( "broken( )" );
985
986 broke.lgBroke = true;
987 return;
988}
989
990/*fixit say that code needs to be fixed */
991void fixit(void)
992{
993 DEBUG_ENTRY( "fixit( )" );
994
995 broke.lgFixit = true;
996 return;
997}
998
999/*CodeReview placed next to code that needs to be checked */
1000void CodeReview(void)
1001{
1002 DEBUG_ENTRY( "CodeReview( )" );
1003
1004 broke.lgCheckit = true;
1005 return;
1006}
1007
1009int dprintf(FILE *fp, const char *format, ...)
1010{
1011 va_list ap;
1012 int i1, i2;
1013
1014 DEBUG_ENTRY( "dprintf()" );
1015 va_start(ap,format);
1016 i1 = fprintf(fp,"DEBUG ");
1017 if (i1 >= 0)
1018 i2 = vfprintf(fp,format,ap);
1019 else
1020 i2 = 0;
1021 if (i2 < 0)
1022 i1 = 0;
1023 va_end(ap);
1024
1025 return i1+i2;
1026}
1027
1028/* dbg_printf is a debug print routine that was provided by Peter Teuben,
1029 * as a component from his NEMO package. It offers run-time specification
1030 * of the level of debugging */
1031int dbg_printf(int debug, const char *fmt, ...)
1032{
1033 va_list ap;
1034 int i=0;
1035
1036 DEBUG_ENTRY( "dbg_printf()" );
1037
1038 /* print this debug message? (debug_level not currently used)*/
1039 if(debug <= trace.debug_level)
1040 {
1041 va_start(ap, fmt);
1042
1043 i = vfprintf(ioQQQ, fmt, ap);
1044 /* drain ioQQQ */
1045 fflush(ioQQQ);
1046 va_end(ap);
1047 }
1048 return i;
1049}
1050
1051
1052/*qg32 32 point Gaussian quadrature, originally given to Gary F by Jim Lattimer */
1053double qg32(
1054 double xl, /*lower limit to integration range*/
1055 double xu, /*upper limit to integration range*/
1056 /*following is the pointer to the function that will be evaluated*/
1057 double (*fct)(double) )
1058{
1059 double a = 0.5*(xu+xl),
1060 b = xu-xl,
1061 y = 0.;
1062
1063 DEBUG_ENTRY( "qg32()" );
1064
1065 /********************************************************************************
1066 * *
1067 * 32-point Gaussian quadrature *
1068 * xl : the lower limit of integration *
1069 * xu : the upper limit *
1070 * fct : the (external) function *
1071 * returns the value of the integral *
1072 * *
1073 * simple call to integrate sine from 0 to pi *
1074 * double agn = qg32( 0., 3.141592654 , sin ); *
1075 * *
1076 *******************************************************************************/
1077
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};
1083
1084 double c[16] = {
1085 .498631930924740780, .49280575577263417, .4823811277937532200, .46745303796886984000,
1086 .448160577883026060, .42468380686628499, .3972418979839712000, .36609105937014484000,
1087 .331522133465107600, .29385787862038116, .2534499544661147000, .21067563806531767000,
1088 .165934301141063820, .11964368112606854, .7223598079139825e-1, .24153832843869158e-1};
1089
1090 for( int i=0; i<16; i++)
1091 {
1092 y += b * weights[i] * ((*fct)(a+b*c[i]) + (*fct)(a-b*c[i]));
1093 }
1094
1095 /* the answer */
1096 return y;
1097}
1098
1099/*spsort netlib routine to sort array returning sorted indices */
1101 /* input array to be sorted */
1102 realnum x[],
1103 /* number of values in x */
1104 long int n,
1105 /* permutation output array */
1106 long int iperm[],
1107 /* flag saying what to do - 1 sorts into increasing order, not changing
1108 * the original vector, -1 sorts into decreasing order. 2, -2 change vector */
1109 int kflag,
1110 /* error condition, should be 0 */
1111 int *ier)
1112{
1113 /*
1114 ****BEGIN PROLOGUE SPSORT
1115 ****PURPOSE Return the permutation vector generated by sorting a given
1116 * array and, optionally, rearrange the elements of the array.
1117 * The array may be sorted in increasing or decreasing order.
1118 * A slightly modified quicksort algorithm is used.
1119 ****LIBRARY SLATEC
1120 ****CATEGORY N6A1B, N6A2B
1121 ****TYPE SINGLE PRECISION (SPSORT-S, DPSORT-D, IPSORT-I, HPSORT-H)
1122 ****KEY WORDS NUMBER SORTING, PASSIVE SORTING, SINGLETON QUICKSORT, SORT
1123 ****AUTHOR Jones, R. E., (SNLA)
1124 * Rhoads, G. S., (NBS)
1125 * Wisniewski, J. A., (SNLA)
1126 ****DESCRIPTION
1127 *
1128 * SPSORT returns the permutation vector IPERM generated by sorting
1129 * the array X and, optionally, rearranges the values in X. X may
1130 * be sorted in increasing or decreasing order. A slightly modified
1131 * quicksort algorithm is used.
1132 *
1133 * IPERM is such that X(IPERM(I)) is the Ith value in the rearrangement
1134 * of X. IPERM may be applied to another array by calling IPPERM,
1135 * SPPERM, DPPERM or HPPERM.
1136 *
1137 * The main difference between SPSORT and its active sorting equivalent
1138 * SSORT is that the data are referenced indirectly rather than
1139 * directly. Therefore, SPSORT should require approximately twice as
1140 * long to execute as SSORT. However, SPSORT is more general.
1141 *
1142 * Description of Parameters
1143 * X - input/output -- real array of values to be sorted.
1144 * If ABS(KFLAG) = 2, then the values in X will be
1145 * rearranged on output; otherwise, they are unchanged.
1146 * N - input -- number of values in array X to be sorted.
1147 * IPERM - output -- permutation array such that IPERM(I) is the
1148 * index of the value in the original order of the
1149 * X array that is in the Ith location in the sorted
1150 * order.
1151 * KFLAG - input -- control parameter:
1152 * = 2 means return the permutation vector resulting from
1153 * sorting X in increasing order and sort X also.
1154 * = 1 means return the permutation vector resulting from
1155 * sorting X in increasing order and do not sort X.
1156 * = -1 means return the permutation vector resulting from
1157 * sorting X in decreasing order and do not sort X.
1158 * = -2 means return the permutation vector resulting from
1159 * sorting X in decreasing order and sort X also.
1160 * IER - output -- error indicator:
1161 * = 0 if no error,
1162 * = 1 if N is zero or negative,
1163 * = 2 if KFLAG is not 2, 1, -1, or -2.
1164 ****REFERENCES R. C. Singleton, Algorithm 347, An efficient algorithm
1165 * for sorting with minimal storage, Communications of
1166 * the ACM, 12, 3 (1969), pp. 185-187.
1167 ****ROUTINES CALLED XERMSG
1168 ****REVISION HISTORY (YYMMDD)
1169 * 761101 DATE WRITTEN
1170 * 761118 Modified by John A. Wisniewski to use the Singleton
1171 * quicksort algorithm.
1172 * 870423 Modified by Gregory S. Rhoads for passive sorting with the
1173 * option for the rearrangement of the original data.
1174 * 890620 Algorithm for rearranging the data vector corrected by R.
1175 * Boisvert.
1176 * 890622 Prologue upgraded to Version 4.0 style by D. Lozier.
1177 * 891128 Error when KFLAG.LT.0 and N=1 corrected by R. Boisvert.
1178 * 920507 Modified by M. McClain to revise prologue text.
1179 * 920818 Declarations section rebuilt and code restructured to use
1180 * IF-THEN-ELSE-ENDIF. (SMR, WRB)
1181 ****END PROLOGUE SPSORT
1182 * .. Scalar Arguments ..
1183 */
1184 long int i,
1185 ij,
1186 il[21],
1187 indx,
1188 indx0,
1189 istrt,
1190 istrt_,
1191 iu[21],
1192 j,
1193 k,
1194 kk,
1195 l,
1196 lm,
1197 lmt,
1198 m,
1199 nn;
1200 realnum r,
1201 ttemp;
1202
1203 DEBUG_ENTRY( "spsort()" );
1204
1205 /* .. Array Arguments .. */
1206 /* .. Local Scalars .. */
1207 /* .. Local Arrays .. */
1208 /* .. External Subroutines .. */
1209 /* .. Intrinsic Functions .. */
1210 /****FIRST EXECUTABLE STATEMENT SPSORT */
1211 *ier = 0;
1212 nn = n;
1213 if( nn < 1 )
1214 {
1215 *ier = 1;
1216 return;
1217 }
1218 else
1219 {
1220 kk = labs(kflag);
1221 if( kk != 1 && kk != 2 )
1222 {
1223 *ier = 2;
1224 return;
1225 }
1226 else
1227 {
1228
1229 /* Initialize permutation vector to index on f scale
1230 * */
1231 for( i=0; i < nn; i++ )
1232 {
1233 iperm[i] = i+1;
1234 }
1235
1236 /* Return if only one value is to be sorted */
1237 if( nn == 1 )
1238 {
1239 --iperm[0];
1240 return;
1241 }
1242
1243 /* Alter array X to get decreasing order if needed */
1244 if( kflag <= -1 )
1245 {
1246 for( i=0; i < nn; i++ )
1247 {
1248 x[i] = -x[i];
1249 }
1250 }
1251
1252 /* Sort X only */
1253 m = 1;
1254 i = 1;
1255 j = nn;
1256 r = .375e0;
1257 }
1258 }
1259
1260 while( true )
1261 {
1262 if( i == j )
1263 goto L_80;
1264 if( r <= 0.5898437e0 )
1265 {
1266 r += 3.90625e-2;
1267 }
1268 else
1269 {
1270 r -= 0.21875e0;
1271 }
1272
1273L_40:
1274 k = i;
1275
1276 /* Select a central element of the array and save it in location L
1277 * */
1278 ij = i + (long)((j-i)*r);
1279 lm = iperm[ij-1];
1280
1281 /* If first element of array is greater than LM, interchange with LM
1282 * */
1283 if( x[iperm[i-1]-1] > x[lm-1] )
1284 {
1285 iperm[ij-1] = iperm[i-1];
1286 iperm[i-1] = lm;
1287 lm = iperm[ij-1];
1288 }
1289 l = j;
1290
1291 /* If last element of array is less than LM, interchange with LM
1292 * */
1293 if( x[iperm[j-1]-1] < x[lm-1] )
1294 {
1295 iperm[ij-1] = iperm[j-1];
1296 iperm[j-1] = lm;
1297 lm = iperm[ij-1];
1298
1299 /* If first element of array is greater than LM, interchange
1300 * with LM
1301 * */
1302 if( x[iperm[i-1]-1] > x[lm-1] )
1303 {
1304 iperm[ij-1] = iperm[i-1];
1305 iperm[i-1] = lm;
1306 lm = iperm[ij-1];
1307 }
1308 }
1309
1310 /* Find an element in the second half of the array which is smaller
1311 * than LM */
1312 while( true )
1313 {
1314 l -= 1;
1315 if( x[iperm[l-1]-1] <= x[lm-1] )
1316 {
1317
1318 /* Find an element in the first half of the array which is greater
1319 * than LM */
1320 while( true )
1321 {
1322 k += 1;
1323 if( x[iperm[k-1]-1] >= x[lm-1] )
1324 break;
1325 }
1326
1327 /* Interchange these elements */
1328 if( k > l )
1329 break;
1330 lmt = iperm[l-1];
1331 iperm[l-1] = iperm[k-1];
1332 iperm[k-1] = lmt;
1333 }
1334 }
1335
1336 /* Save upper and lower subscripts of the array yet to be sorted */
1337 if( l - i > j - k )
1338 {
1339 il[m-1] = i;
1340 iu[m-1] = l;
1341 i = k;
1342 m += 1;
1343 }
1344 else
1345 {
1346 il[m-1] = k;
1347 iu[m-1] = j;
1348 j = l;
1349 m += 1;
1350 }
1351
1352L_90:
1353 if( j - i >= 1 )
1354 goto L_40;
1355 if( i == 1 )
1356 continue;
1357 i -= 1;
1358
1359 while( true )
1360 {
1361 i += 1;
1362 if( i == j )
1363 break;
1364 lm = iperm[i];
1365 if( x[iperm[i-1]-1] > x[lm-1] )
1366 {
1367 k = i;
1368
1369 while( true )
1370 {
1371 iperm[k] = iperm[k-1];
1372 k -= 1;
1373
1374 if( x[lm-1] >= x[iperm[k-1]-1] )
1375 break;
1376 }
1377 iperm[k] = lm;
1378 }
1379 }
1380
1381 /* Begin again on another portion of the unsorted array */
1382L_80:
1383 m -= 1;
1384 if( m == 0 )
1385 break;
1386 /*lint -e644 not explicitly initialized */
1387 i = il[m-1];
1388 j = iu[m-1];
1389 /*lint +e644 not explicitly initialized */
1390 goto L_90;
1391 }
1392
1393 /* Clean up */
1394 if( kflag <= -1 )
1395 {
1396 for( i=0; i < nn; i++ )
1397 {
1398 x[i] = -x[i];
1399 }
1400 }
1401
1402 /* Rearrange the values of X if desired */
1403 if( kk == 2 )
1404 {
1405
1406 /* Use the IPERM vector as a flag.
1407 * If IPERM(I) < 0, then the I-th value is in correct location */
1408 for( istrt=1; istrt <= nn; istrt++ )
1409 {
1410 istrt_ = istrt - 1;
1411 if( iperm[istrt_] >= 0 )
1412 {
1413 indx = istrt;
1414 indx0 = indx;
1415 ttemp = x[istrt_];
1416 while( iperm[indx-1] > 0 )
1417 {
1418 x[indx-1] = x[iperm[indx-1]-1];
1419 indx0 = indx;
1420 iperm[indx-1] = -iperm[indx-1];
1421 indx = labs(iperm[indx-1]);
1422 }
1423 x[indx0-1] = ttemp;
1424 }
1425 }
1426
1427 /* Revert the signs of the IPERM values */
1428 for( i=0; i < nn; i++ )
1429 {
1430 iperm[i] = -iperm[i];
1431 }
1432 }
1433
1434 for( i=0; i < nn; i++ )
1435 {
1436 --iperm[i];
1437 }
1438 return;
1439}
1440
1441/*MyMalloc wrapper for malloc(). Returns a good pointer or dies.
1442 * memory is filled with NaN
1443 * >>chng 05 dec 14, do not set to NaN since tricks debugger
1444 * routines within code do not call this or malloc, but rather MALLOC
1445 * which is resolved into MyMalloc or malloc depending on whether
1446 * NDEBUG is set by the compiler to indicate "not debugging",
1447 * in typical negative C style */
1449 /*use same type as library function MALLOC*/
1450 size_t size ,
1451 const char *chFile,
1452 int line
1453 )
1454{
1455 void *ptr;
1456
1457 DEBUG_ENTRY( "MyMalloc()" );
1458
1459 ASSERT( size > 0 );
1460
1461 /* debug branch for printing malloc args */
1462 {
1463 enum{DEBUG_LOC=false};
1464 if( DEBUG_LOC)
1465 {
1466 static long int kount=0, nTot=0;
1467 nTot += (long)size;
1468 fprintf(ioQQQ,"%li\t%li\t%li\n",
1469 kount ,
1470 (long)size ,
1471 nTot );
1472 ++kount;
1473 }
1474 }
1475
1476 if( ( ptr = malloc( size ) ) == NULL )
1477 {
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",
1481 chFile , line );
1482
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",
1486 struc.nzlim );
1487
1489 }
1490
1491 /* flag -DNOINIT will turn off this initialization which can fool valgrind/purify */
1492# if !defined(NDEBUG) && !defined(NOINIT)
1493
1494 size_t nFloat = size/4;
1495 size_t nDouble = size/8;
1496 sys_float *fptr = static_cast<sys_float*>(ptr);
1497 double *dptr = static_cast<double*>(ptr);
1498
1499 /* >>chng 04 feb 03, fill memory with invalid numbers, PvH */
1500 /* on IA32/AMD64 processors this will generate NaN's for both float and double;
1501 * on most other (modern) architectures it is likely to do the same... */
1502 /* >>chng 05 dec 14, change code to generate signaling NaN's for most cases (but not all!) */
1503 if( size == nDouble*8 )
1504 {
1505 /* this could be an array of doubles as well as floats -> we will hedge our bets
1506 * we will fill the array with a pattern that is interpreted as all signaling
1507 * NaN's for doubles, and alternating signaling and quiet NaN's for floats:
1508 * byte offset: 0 4 8 12 16
1509 * double | SNaN | SNan |
1510 * float | SNaN | QNaN | SNan | QNaN | (little-endian, e.g. Intel, AMD, alpha)
1511 * float | QNaN | SNaN | QNan | SNaN | (big-endian, e.g. Sparc, PowerPC, MIPS) */
1512 set_NaN( dptr, (long)nDouble );
1513 }
1514 else if( size == nFloat*4 )
1515 {
1516 /* this could be an arrays of floats, but not doubles -> init to all float SNaN */
1517 set_NaN( fptr, (long)nFloat );
1518 }
1519 else
1520 {
1521 memset( ptr, 0xff, size );
1522 }
1523
1524# endif /* !defined(NDEBUG) && !defined(NOINIT) */
1525 return ptr;
1526}
1527
1528
1529/* wrapper for calloc(). Returns a good pointer or dies.
1530 * routines within code do not call this or malloc, but rather CALLOC
1531 * which is resolved into MyCalloc or calloc depending on whether
1532 * NDEBUG is set in cddefines. \h */
1534 /*use same type as library function CALLOC*/
1535 size_t num ,
1536 size_t size )
1537{
1538 void *ptr;
1539
1540 DEBUG_ENTRY( "MyCalloc()" );
1541
1542 ASSERT( size > 0 );
1543
1544 /* debug branch for printing malloc args */
1545 {
1546 enum{DEBUG_LOC=false};
1547 if( DEBUG_LOC)
1548 {
1549 static long int kount=0;
1550 fprintf(ioQQQ,"%li\tcall\t%li\tbytes\n", kount,
1551 (long)size );
1552 ++kount;
1553 }
1554 }
1555
1556 if( ( ptr = calloc( num , size ) ) == NULL )
1557 {
1558 fprintf(ioQQQ,"MyCalloc could not allocate %lu bytes. Exit in MyCalloc.",
1559 (unsigned long)size );
1561 }
1562 return ptr;
1563}
1564
1565/* wrapper for realloc(). Returns a good pointer or dies.
1566 * routines within code do not call this or malloc, but rather REALLOC
1567 * which is resolved into MyRealloc or realloc depending on whether
1568 * NDEBUG is set in cddefines.h */
1570 /*use same type as library function realloc */
1571 void *p ,
1572 size_t size )
1573{
1574 void *ptr;
1575
1576 DEBUG_ENTRY( "MyRealloc()" );
1577
1578 ASSERT( size > 0 );
1579
1580 /* debug branch for printing malloc args */
1581 {
1582 enum{DEBUG_LOC=false};
1583 if( DEBUG_LOC)
1584 {
1585 static long int kount=0;
1586 fprintf(ioQQQ,"%li\tcall\t%li\tbytes\n", kount,
1587 (long)size );
1588 ++kount;
1589 }
1590 }
1591
1592 if( ( ptr = realloc( p , size ) ) == NULL )
1593 {
1594 fprintf(ioQQQ,"MyRealloc could not allocate %lu bytes. Exit in MyRealloc.",
1595 (unsigned long)size );
1597 }
1598 return ptr;
1599}
1600
1601/* function to facilitate addressing opacity array */
1602double csphot(
1603 /* INU is array index pointing to frequency where opacity is to be evaluated
1604 * on f not c scale */
1605 long int inu,
1606 /* ITHR is pointer to threshold*/
1607 long int ithr,
1608 /* IOFSET is offset as defined in opac0*/
1609 long int iofset)
1610{
1611 double csphot_v;
1612
1613 DEBUG_ENTRY( "csphot()" );
1614
1615 csphot_v = opac.OpacStack[inu-ithr+iofset-1];
1616 return csphot_v;
1617}
1618
1619/*RandGauss normal Gaussian random number generator
1620 * the user must call srand to set the seed before using this routine.
1621
1622 * the random numbers will have a mean of xMean and a standard deviation of s
1623
1624 * The convention is for srand to be called when the command setting
1625 * the noise is parsed
1626
1627 * for very small dispersion there are no issues, but when the dispersion becomes
1628 * large the routine will find negative values - so most often used in this case
1629 * to find dispersion in log space
1630 * this routine will return a normal Gaussian - must be careful in how this is
1631 * used when adding noise to physical quantity */
1632/*
1633NB - following from Ryan Porter:
1634I discovered that I unintentionally created an antisymmetric skew in my
1635Monte Carlo. RandGauss is symmetric in log space, which means it is not
1636symmetric in linear space. But to get the right standard deviation you
1637have to take 10^x, where x is the return from RandGauss. The problem is
163810^x will happen less frequently than 10^-x, so without realizing it, the
1639average "tweak" to every piece of atomic data in my Monte Carlo run was
1640not 1.0 but something greater than 1.0, causing every single line to have
1641an average Monte Carlo emissivity greater than the regular value. Any place
1642that takes 10^RandGauss() needs to be adjusted if what is intended is +/- x. */
1644 /* mean value */
1645 double xMean,
1646 /*standard deviation s */
1647 double s )
1648{
1649 double x1, x2, w, yy1;
1650 static double yy2=-BIGDOUBLE;
1651 static int use_last = false;
1652
1653 DEBUG_ENTRY( "RandGauss()" );
1654
1655 if( use_last )
1656 {
1657 yy1 = yy2;
1658 use_last = false;
1659 }
1660 else
1661 {
1662 do {
1663 x1 = 2.*genrand_real3() - 1.;
1664 x2 = 2.*genrand_real3() - 1.;
1665 w = x1 * x1 + x2 * x2;
1666 } while ( w >= 1.0 );
1667
1668 w = sqrt((-2.0*log(w))/w);
1669 yy1 = x1 * w;
1670 yy2 = x2 * w;
1671 use_last = true;
1672 }
1673 return xMean + yy1 * s;
1674}
1675
1676/* MyGaussRand takes as input a percent uncertainty less than 50%
1677 * (expressed as 0.5). The routine then assumes this input variable represents one
1678 * standard deviation about a mean of unity, and returns a random number within
1679 * that range. A hard cutoff is imposed at two standard deviations, which
1680 * eliminates roughly 5% of the normal distribution. In other words, the routine
1681 * returns a number in a normal distribution with standard deviation equal to
1682 * the input. The number will be between 1-3*stdev and 1+3*stdev. */
1683double MyGaussRand( double PctUncertainty )
1684{
1685 double StdDev;
1686 double result;
1687
1688 DEBUG_ENTRY( "MyGaussRand()" );
1689
1690 ASSERT( PctUncertainty < 0.5 );
1691 /* We want this "percent uncertainty" to represent one standard deviation */
1692 StdDev = PctUncertainty;
1693
1694 do
1695 {
1696 /*result = pow( 10., RandGauss( 0., logStdDev ) );*/
1697 result = 1.+RandGauss( 0., StdDev );
1698 }
1699 /* only allow values that are within 3 standard deviations */
1700 while( (result < 1.-3.*PctUncertainty) || (result > 1.+3.*PctUncertainty) );
1701
1702 ASSERT( result>0. && result<2. );
1703 return result;
1704}
1705
1706/*plankf evaluate Planck function for any cell at current electron temperature */
1707double plankf(long int ip)
1708{
1709 double plankf_v;
1710
1711 DEBUG_ENTRY( "plankf()" );
1712
1713 /* evaluate Planck function
1714 * argument is pointer to cell energy in ANU
1715 * return photon flux for cell IP */
1716 if( rfield.ContBoltz[ip] <= 0. )
1717 {
1718 plankf_v = 1e-35;
1719 }
1720 else
1721 {
1722 plankf_v = 6.991e-21*POW2(FR1RYD*rfield.anu[ip])/
1723 (1./rfield.ContBoltz[ip] - 1.)*FR1RYD*4.;
1724 }
1725 return plankf_v;
1726}
1727
1729{
1730 fstream io;
1731 string line;
1732 open_data( io, "citation_cloudy.txt", mode_r, AS_DATA_ONLY );
1733 while( SafeGetline( io, line ) )
1734 {
1735 if( line[0] == '#' )
1736 continue;
1737 // replace XXXX with actual version number
1738 size_t p = line.find( "XXXX" );
1739 if( p != string::npos )
1740 line.replace( p, 4, t_version::Inst().chVersion );
1741 fprintf( ioQQQ, "%s\n", line.c_str() );
1742 }
1743}
1744
1746{
1747 fstream io;
1748 string line;
1749 open_data( io, "citation_data.txt", mode_r, AS_DATA_ONLY );
1750 while( SafeGetline( io, line ) )
1751 {
1752 if( line[0] == '#' )
1753 continue;
1754 // replace XXXX with actual version number
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() );
1759 }
1760}
1761
1762// this routine was taken from
1763// http://stackoverflow.com/questions/6089231/getting-std-ifstream-to-handle-lf-cr-and-crlf
1764// it is copyrighted by a creative commons license
1765// http://creativecommons.org/licenses/by-sa/3.0/
1766//
1767// safe version of getline() that correctly handles all types of EOL lf, crlf and cr...
1768// it has been modified such that it does not produce a spurious empty line at the end of a file
1769// this way it is compatible with the standard getline() (at least with g++/linux).
1770istream& SafeGetline(istream& is, string& t)
1771{
1772 t.clear();
1773
1774 // The characters in the stream are read one-by-one using a streambuf.
1775 // That is faster than reading them one-by-one using the istream.
1776 // Code that uses streambuf this way must be guarded by a sentry object.
1777 // The sentry object performs various tasks,
1778 // such as thread synchronization and updating the stream state.
1779
1780 istream::sentry se(is, true);
1781 streambuf* sb = is.rdbuf();
1782
1783 while( true )
1784 {
1785 int c = sb->sbumpc();
1786 switch (c)
1787 {
1788 case '\n':
1789 if( sb->sgetc() == EOF )
1790 is.setstate(ios::eofbit);
1791 return is;
1792 case '\r':
1793 if( sb->sgetc() == '\n' )
1794 sb->sbumpc();
1795 if( sb->sgetc() == EOF )
1796 is.setstate(ios::eofbit);
1797 return is;
1798 case EOF:
1799 // Also handle the case when the last line has no line ending
1800 is.setstate(ios::eofbit);
1801 return is;
1802 default:
1803 t += (char)c;
1804 }
1805 }
1806}
t_atmdat atmdat
Definition atmdat.cpp:6
static double x2[63]
static double x1[83]
t_broke broke
Definition broke.cpp:5
t_called called
Definition called.cpp:5
long int nzone
Definition cddefines.cpp:14
bool lgTestCodeCalled
Definition cddefines.cpp:11
FILE * ioQQQ
Definition cddefines.cpp:7
bool lgAbort
Definition cddefines.cpp:10
long int iteration
Definition cddefines.cpp:16
double fnzone
Definition cddefines.cpp:15
float sys_float
Definition cddefines.h:106
#define ASSERT(exp)
Definition cddefines.h:578
bool is_odd(int j)
Definition cddefines.h:714
#define PrintEfmt(F, V)
Definition cddefines.h:1472
const double DSEXP_LIMIT
Definition cddefines.h:1478
const double SEXP_LIMIT
Definition cddefines.h:1476
char tolower(char c)
Definition cddefines.h:691
#define POW2
Definition cddefines.h:929
#define EXIT_FAILURE
Definition cddefines.h:140
char toupper(char c)
Definition cddefines.h:700
char TorF(bool l)
Definition cddefines.h:710
#define cdEXIT(FAIL)
Definition cddefines.h:434
const char * strstr_s(const char *haystack, const char *needle)
Definition cddefines.h:1429
float realnum
Definition cddefines.h:103
const char * strchr_s(const char *s, int c)
Definition cddefines.h:1439
long max(int a, long b)
Definition cddefines.h:775
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long min(int a, long b)
Definition cddefines.h:723
split_mode
Definition cddefines.h:1321
@ SPM_STRICT
Definition cddefines.h:1321
@ SPM_KEEP_EMPTY
Definition cddefines.h:1321
void cdCautions(FILE *ioOUT)
Definition cddrive.cpp:226
void cdPrintCommands(FILE *ioOUT)
Definition cddrive.cpp:513
void cdWarnings(FILE *ioPNT)
Definition cddrive.cpp:198
Definition energy.h:8
double get(const char *unit) const
Definition energy.cpp:141
static t_version & Inst()
Definition cddefines.h:175
t_conv conv
Definition conv.cpp:5
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
void set_NaN(sys_float &x)
Definition cpu.cpp:682
const ios_base::openmode mode_r
Definition cpu.h:212
@ AS_DATA_ONLY
Definition cpu.h:207
const double BIGDOUBLE
Definition cpu.h:194
#define NORETURN
Definition cpu.h:383
t_fudgec fudgec
Definition fudgec.cpp:5
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
t_hextra hextra
Definition hextra.cpp:5
t_input input
Definition input.cpp:12
molezone * findspecieslocal(const char buf[])
t_opac opac
Definition opacity.cpp:5
UNUSED const double FR1RYD
Definition physconst.h:195
t_rfield rfield
Definition rfield.cpp:8
t_save save
Definition save.cpp:5
double csphot(long int inu, long int ithr, long int iofset)
Definition service.cpp:1602
double plankf(long int ip)
Definition service.cpp:1707
sys_float sexp(sys_float x)
Definition service.cpp:914
double MyGaussRand(double PctUncertainty)
Definition service.cpp:1683
long ipow(long m, long n)
Definition service.cpp:639
long nMatch(const char *chKey, const char *chCard)
Definition service.cpp:451
double ee1_safe(double x)
Definition service.cpp:356
void * MyCalloc(size_t num, size_t size)
Definition service.cpp:1533
double powi(double x, long int n)
Definition service.cpp:604
double qg32(double xl, double xu, double(*fct)(double))
Definition service.cpp:1053
void PrintE82(FILE *ioOUT, double value)
Definition service.cpp:739
void broken(void)
Definition service.cpp:982
NORETURN void BadRead(void)
Definition service.cpp:901
void TestCode(void)
Definition service.cpp:972
void CodeReview(void)
Definition service.cpp:1000
double fudge(long int ipnt)
Definition service.cpp:481
double ee1(double x)
Definition service.cpp:312
int dbg_printf(int debug, const char *fmt,...)
Definition service.cpp:1031
int dprintf(FILE *fp, const char *format,...)
Definition service.cpp:1009
double RandGauss(double xMean, double s)
Definition service.cpp:1643
void cap4(char *chCAP, const char *chLab)
Definition service.cpp:240
void MyAssert(const char *file, int line, const char *comment)
Definition service.cpp:153
int GetQuote(char *chStringOut, char *chCard, char *chCardRaw, bool lgAbort)
Definition service.cpp:513
void Split(const string &str, const string &sep, vector< string > &lst, split_mode mode)
Definition service.cpp:106
void PrintE71(FILE *ioOUT, double value)
Definition service.cpp:788
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition service.cpp:381
void * MyRealloc(void *p, size_t size)
Definition service.cpp:1569
void * MyMalloc(size_t size, const char *chFile, int line)
Definition service.cpp:1448
void CloudyPrintReference()
Definition service.cpp:1728
void PrintE93(FILE *ioOUT, double value)
Definition service.cpp:838
void DatabasePrintReference()
Definition service.cpp:1745
void caps(char *chCard)
Definition service.cpp:280
double e2(double t)
Definition service.cpp:299
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
void spsort(realnum x[], long int n, long int iperm[], int kflag, int *ier)
Definition service.cpp:1100
NORETURN void TotalInsanity(void)
Definition service.cpp:886
double AnuUnit(realnum energy_ryd)
Definition service.cpp:173
double dsexp(double x)
Definition service.cpp:953
void uncaps(char *chCard)
Definition service.cpp:263
void ShowMe(void)
Definition service.cpp:181
void fixit(void)
Definition service.cpp:991
istream & SafeGetline(istream &is, string &t)
Definition service.cpp:1770
t_struc struc
Definition struc.cpp:6
double genrand_real3()
t_trace trace
Definition trace.cpp:5
t_warnings warnings
Definition warnings.cpp:11