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

Generated for cloudy by doxygen 1.8.3.1