cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /* 00004 * a set of routines that are widely used across the code for various 00005 * housekeeping chores. These do not do any physics and are unlikely to 00006 * change over time. The prototypes are in cddefines.h and so are 00007 * automatically picked up by all routines 00008 */ 00009 /*FFmtRead scan input line for free format number */ 00010 /*e2 second exponential integral */ 00011 /*caps convert input command line (through eol) to ALL CAPS */ 00012 /*ShowMe produce request to send information to GJF after a crash */ 00013 /*AnuUnit produce continuum energy in arbitrary units */ 00014 /*cap4 convert first 4 char of input line chLab into chCAP all in caps, null termination */ 00015 /*insane set flag saying that insanity has occurred */ 00016 /*nMatch determine whether match to a keyword occurs on command line, 00017 * return value is 0 if no match, and position of match within string if hit */ 00018 /*fudge enter fudge factors, or some arbitrary number, with fudge command*/ 00019 /*GetElem scans line image, finds element. returns atomic number j, on C scale */ 00020 /*GetQuote get any name between double quotes off command line 00021 * return string as chLabel, is null terminated */ 00022 /*qip compute pow(x,n) for positive integer n through repeated squares */ 00023 /*NoNumb general error handler for no numbers on input line */ 00024 /*dsexp safe exponential function for doubles */ 00025 /*sexp safe exponential function */ 00026 /*TestCode set flag saying that test code is in place */ 00027 /*CodeReview - placed next to code that needs to be checked */ 00028 /*fixit - say that code needs to be fixed */ 00029 /*broken set flag saying that the code is broken, */ 00030 /*dbg_printf is a debug print routine that was provided by Peter Teuben, 00031 * as a component from his NEMO package. It offers run-time specification 00032 * of the level of debugging */ 00033 /*qg32 32 point Gaussian quadrature, original Fortran given to Gary F by Jim Lattimer */ 00034 /*TotalInsanity general error handler for something that cannot happen */ 00035 /*BadRead general error handler for trying to read data, but failing */ 00036 /*BadOpen general error handler for trying to open file, but failing */ 00037 /*MyMalloc wrapper for malloc(). Returns a good pointer or dies. */ 00038 /*MyCalloc wrapper for calloc(). Returns a good pointer or dies. */ 00039 /*spsort netlib routine to sort array returning sorted indices */ 00040 /*chLineLbl use information in line transfer arrays to generate a line label * 00041 * this label is null terminated */ 00042 /*chIonLbl use information in line array to generate a null terminated ion label in "Fe 2" */ 00043 /*csphot returns photoionization cross section from opacity stage using std pointers */ 00044 /*MyAssert a version of assert that fails gracefully */ 00045 /*RandGauss normal random variate generator */ 00046 /*MyGaussRand a wrapper for RandGauss, see below */ 00047 #include <ctype.h> 00048 #include <stdarg.h> /* ANSI variable arg macros */ 00049 #include "cddefines.h" 00050 #include "physconst.h" 00051 #include "cddrive.h" 00052 #include "called.h" 00053 #include "opacity.h" 00054 #include "rfield.h" 00055 #include "hextra.h" 00056 #include "hmi.h" 00057 #include "fudgec.h" 00058 #include "broke.h" 00059 #include "trace.h" 00060 #include "input.h" 00061 #include "elementnames.h" 00062 #include "punch.h" 00063 #include "version.h" 00064 #include "warnings.h" 00065 #include "conv.h" 00066 #include "thirdparty.h" 00067 00068 /*read_whole_line - safe version of fgets - read a line, 00069 * return null if cannot read line or if input line is too long */ 00070 char *read_whole_line( char *chLine , int nChar , FILE *ioIN ) 00071 { 00072 char *chRet , 00073 *chEOL; 00074 00075 DEBUG_ENTRY( "read_whole_line()" ); 00076 00077 if( (chRet = fgets( chLine , nChar, ioIN )) == NULL ) 00078 { 00079 return NULL; 00080 } 00081 00082 /* check that there are less than nChar characters in the line */ 00083 /*chEOL = strchr(chLine , '\0' );*/ 00084 chEOL = (char*)memchr( chLine , '\0' , nChar ); 00085 00086 /* return null if input string longer than nChar, the longest we can read. 00087 * Print and return null but chLine still has as much of the line as 00088 * could be placed in cdLine */ 00089 if( (chEOL==NULL) || (chEOL - chLine)>=nChar-1 ) 00090 { 00091 if( called.lgTalk ) 00092 fprintf( ioQQQ, "DISASTER PROBLEM read_whole_line found input" 00093 " with a line too long to be read.\n" ); 00094 00095 lgAbort = true; 00096 return NULL; 00097 } 00098 return chRet; 00099 } 00100 00102 void Split(const string& str, // input string 00103 const string& sep, // separator, may be multiple characters 00104 vector<string>& lst, // the separated items will be appended here 00105 split_mode mode) // SPM_RELAX, SPM_KEEP_EMPTY, or SPM_STRICT; see cddefines.h 00106 { 00107 DEBUG_ENTRY( "Split()" ); 00108 00109 bool lgStrict = ( mode == SPM_STRICT ); 00110 bool lgKeep = ( mode == SPM_KEEP_EMPTY ); 00111 bool lgFail = false; 00112 string::size_type ptr1 = 0; 00113 string::size_type ptr2 = str.find( sep ); 00114 string sstr = str.substr( ptr1, ptr2-ptr1 ); 00115 if( sstr.length() > 0 ) 00116 lst.push_back( sstr ); 00117 else { 00118 if( lgStrict ) lgFail = true; 00119 if( lgKeep ) lst.push_back( sstr ); 00120 } 00121 while( ptr2 != string::npos ) { 00122 // the separator is skipped 00123 ptr1 = ptr2 + sep.length(); 00124 if( ptr1 < str.length() ) { 00125 ptr2 = str.find( sep, ptr1 ); 00126 sstr = str.substr( ptr1, ptr2-ptr1 ); 00127 if( sstr.length() > 0 ) 00128 lst.push_back( sstr ); 00129 else { 00130 if( lgStrict ) lgFail = true; 00131 if( lgKeep ) lst.push_back( sstr ); 00132 } 00133 } 00134 else { 00135 ptr2 = string::npos; 00136 if( lgStrict ) lgFail = true; 00137 if( lgKeep ) lst.push_back( "" ); 00138 } 00139 } 00140 if( lgFail ) 00141 { 00142 fprintf( ioQQQ, " A syntax error occurred while splitting the string: \"%s\"\n", str.c_str() ); 00143 fprintf( ioQQQ, " The separator is \"%s\". Empty substrings are not allowed.\n", sep.c_str() ); 00144 cdEXIT(EXIT_FAILURE); 00145 } 00146 } 00147 00148 /* a version of assert that fails gracefully */ 00149 void MyAssert(const char *file, int line) 00150 { 00151 DEBUG_ENTRY( "MyAssert()" ); 00152 00153 fprintf(ioQQQ," PROBLEM DISASTER An assert has been thrown, this is bad.\n"); 00154 fprintf(ioQQQ," It happened in the file %s at line number %i\n", file, line ); 00155 fprintf(ioQQQ," This is iteration %li, nzone %li, fzone %.2f, lgSearch=%c.\n", 00156 iteration , 00157 nzone , 00158 fnzone , 00159 TorF(conv.lgSearch) ); 00160 00161 ShowMe(); 00162 # if defined ASSERTDEBUG 00163 cdEXIT(EXIT_FAILURE); 00164 # endif 00165 } 00166 00167 /*AnuUnit produce continuum energy in arbitrary units */ 00168 double AnuUnit(realnum energy_ryd) 00169 /*static double AnuUnit(long int ip)*/ 00170 { 00171 double AnuUnit_v; 00172 00173 DEBUG_ENTRY( "AnuUnit()" ); 00174 00175 /* energy comes in in Ryd */ 00176 if( energy_ryd <=0. ) 00177 { 00178 /* this is insanity */ 00179 AnuUnit_v = 0.; 00180 } 00181 else if( strcmp(punch.chConPunEnr[punch.ipConPun],"ryd ") == 0 ) 00182 { 00183 /* energy in Ryd */ 00184 AnuUnit_v = energy_ryd; 00185 } 00186 else if( strcmp(punch.chConPunEnr[punch.ipConPun],"micr") == 0 ) 00187 { 00188 /* wavelength in microns */ 00189 AnuUnit_v = RYDLAM/energy_ryd*1e-4; 00190 } 00191 else if( strcmp(punch.chConPunEnr[punch.ipConPun]," kev") == 0 ) 00192 { 00193 /* energy in keV */ 00194 AnuUnit_v = energy_ryd*EVRYD*1.e-3; 00195 } 00196 else if( strcmp(punch.chConPunEnr[punch.ipConPun]," ev ") == 0 ) 00197 { 00198 /* energy in eV */ 00199 AnuUnit_v = energy_ryd*EVRYD; 00200 } 00201 else if( strcmp(punch.chConPunEnr[punch.ipConPun],"angs") == 0 ) 00202 { 00203 /* wavelength in Angstroms */ 00204 AnuUnit_v = RYDLAM/energy_ryd; 00205 } 00206 else if( strcmp(punch.chConPunEnr[punch.ipConPun],"cent") == 0 ) 00207 { 00208 /* wavelength in centimeters */ 00209 AnuUnit_v = RYDLAM/energy_ryd*1e-8; 00210 } 00211 else if( strcmp(punch.chConPunEnr[punch.ipConPun],"wave") == 0 ) 00212 { 00213 /* wavenumbers */ 00214 AnuUnit_v = RYD_INF*energy_ryd; 00215 } 00216 else if( strcmp(punch.chConPunEnr[punch.ipConPun]," mhz") == 0 ) 00217 { 00218 /* MHz */ 00219 AnuUnit_v = RYD_INF*energy_ryd*SPEEDLIGHT/1e6; 00220 } 00221 else 00222 { 00223 fprintf( ioQQQ, " insane units in AnuUnit =%4.4s\n", 00224 punch.chConPunEnr[punch.ipConPun] ); 00225 cdEXIT(EXIT_FAILURE); 00226 } 00227 00228 return AnuUnit_v; 00229 } 00230 00231 /*ShowMe produce request to send information to GJF after a crash */ 00232 void ShowMe(void) 00233 { 00234 00235 DEBUG_ENTRY( "ShowMe()" ); 00236 00237 /* print info if output unit is defined */ 00238 if( ioQQQ != NULL ) 00239 { 00240 /* >>chng 06 mar 02 - check if molecular but cosmic rays are ignored */ 00241 if( (hextra.cryden == 0.) && hmi.BiggestH2 > 0.1 ) 00242 { 00243 fprintf( ioQQQ, " >>> \n >>> \n >>> Cosmic rays are not included and the gas is molecular. " 00244 "THIS IS KNOWN TO BE UNSTABLE. Add cosmic rays and try again.\n >>> \n >>>\n\n"); 00245 } 00246 else 00247 { 00248 fprintf( ioQQQ, "\n\n\n" ); 00249 fprintf( ioQQQ, " vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv \n" ); 00250 fprintf( ioQQQ, " > PROBLEM DISASTER PROBLEM DISASTER. <\n" ); 00251 fprintf( ioQQQ, " > Sorry, something bad has happened. <\n" ); 00252 fprintf( ioQQQ, " > Please post this on the Cloudy web site <\n" ); 00253 fprintf( ioQQQ, " > discussion board at www.nublado.org <\n" ); 00254 fprintf( ioQQQ, " > Please send all following information: <\n" ); 00255 fprintf( ioQQQ, " ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ \n" ); 00256 fprintf( ioQQQ, "\n\n" ); 00257 00258 00259 fprintf( ioQQQ, " Cloudy version number is %s\n", 00260 t_version::Inst().chVersion ); 00261 fprintf( ioQQQ, " %s\n\n", t_version::Inst().chInfo ); 00262 00263 fprintf( ioQQQ, "%5ld warnings,%3ld cautions,%3ld temperature failures. Messages follow.\n", 00264 warnings.nwarn, warnings.ncaun, conv.nTeFail ); 00265 00266 /* print the warnings first */ 00267 cdWarnings(ioQQQ); 00268 00269 /* now print the cautions */ 00270 cdCautions(ioQQQ); 00271 00272 /* now output the commands */ 00273 cdPrintCommands(ioQQQ); 00274 00275 /* if init command was present, this is the number of lines in it - 00276 * if no init then still set to zero as done in cdInit */ 00277 if( input.nSaveIni ) 00278 { 00279 fprintf(ioQQQ," This input stream included an init file.\n"); 00280 fprintf(ioQQQ," If this init file is not part of the standard Cloudy distribution\n"); 00281 fprintf(ioQQQ," then I will need a copy of it too.\n"); 00282 } 00283 } 00284 } 00285 return; 00286 } 00287 00288 /*cap4 convert first 4 char of input line chLab into chCAP all in caps, null termination */ 00289 void cap4( 00290 char *chCAP , /* output string, cap'd first 4 char of chLab, */ 00291 /* with null terminating */ 00292 char *chLab) /* input string ending with eol*/ 00293 { 00294 long int /*chr,*/ 00295 i; 00296 00297 DEBUG_ENTRY( "cap4()" ); 00298 00299 /* convert 4 character string in chLab to ALL CAPS in chCAP */ 00300 for( i=0; i < 4; i++ ) 00301 { 00302 /* toupper is function in ctype that converts to upper case */ 00303 chCAP[i] = (char)toupper( chLab[i] ); 00304 } 00305 00306 /* now end string with eol */ 00307 chCAP[4] = '\0'; 00308 return; 00309 } 00310 00311 /*uncaps convert input command line (through eol) to all lowercase */ 00312 void uncaps(char *chCard ) 00313 { 00314 long int i; 00315 00316 DEBUG_ENTRY( "caps()" ); 00317 00318 /* convert full character string in chCard to ALL CAPS */ 00319 i = 0; 00320 while( chCard[i]!= '\0' ) 00321 { 00322 chCard[i] = (char)tolower( chCard[i] ); 00323 ++i; 00324 } 00325 return; 00326 } 00327 00328 /*caps convert input command line (through eol) to ALL CAPS */ 00329 void caps(char *chCard ) 00330 { 00331 long int i; 00332 00333 DEBUG_ENTRY( "caps()" ); 00334 00335 /* convert full character string in chCard to ALL CAPS */ 00336 i = 0; 00337 while( chCard[i]!= '\0' ) 00338 { 00339 chCard[i] = (char)toupper( chCard[i] ); 00340 ++i; 00341 } 00342 return; 00343 } 00344 00345 /*e2 second exponential integral */ 00346 /*>>chng 07 jan 17, PvH discover that exp-t is not really 00347 * exp-t - this changed results in several tests */ 00348 double e2( 00349 /* the argument to E2 */ 00350 double t ) 00351 { 00352 /* use recurrence relation */ 00353 /* ignore exp_mt, it is *very* unreliable */ 00354 double hold = sexp(t) - t*ee1(t); 00355 DEBUG_ENTRY( "e2()" ); 00356 /* guard against negative results, this can happen for very large t */ 00357 return max( hold, 0. ); 00358 } 00359 00360 /*ee1 first exponential integral */ 00361 double ee1(double x) 00362 { 00363 double ans, 00364 bot, 00365 top; 00366 static double a[6]={-.57721566,.99999193,-.24991055,.05519968,-.00976004, 00367 .00107857}; 00368 static double b[4]={8.5733287401,18.0590169730,8.6347608925,.2677737343}; 00369 static double c[4]={9.5733223454,25.6329561486,21.0996530827,3.9584969228}; 00370 00371 DEBUG_ENTRY( "ee1()" ); 00372 00373 /* computes the exponential integral E1(x), 00374 * from Abramowitz and Stegun 00375 * stops with error condition for negative argument, 00376 * returns zero in large x limit 00377 * */ 00378 00379 /* error - does not accept negative arguments */ 00380 if( x <= 0 ) 00381 { 00382 fprintf( ioQQQ, " DISASTER negative argument in function ee1, x<0\n" ); 00383 cdEXIT(EXIT_FAILURE); 00384 } 00385 00386 /* branch for x less than 1 */ 00387 else if( x < 1. ) 00388 { 00389 /* abs. accuracy better than 2e-7 */ 00390 ans = ((((a[5]*x + a[4])*x + a[3])*x + a[2])*x + a[1])*x + a[0] - log(x); 00391 } 00392 00393 /* branch for x greater than, or equal to, one */ 00394 else 00395 { 00396 /* abs. accuracy better than 2e-8 */ 00397 top = (((x + b[0])*x + b[1])*x + b[2])*x + b[3]; 00398 bot = (((x + c[0])*x + c[1])*x + c[2])*x + c[3]; 00399 ans = top/bot/x*exp(-x); 00400 } 00401 return ans; 00402 } 00403 00404 /* same as ee1, except without factor of exp(x) in returned value */ 00405 double ee1_safe(double x) 00406 { 00407 double ans, 00408 bot, 00409 top; 00410 /*static double a[6]={-.57721566,.99999193,-.24991055,.05519968,-.00976004, 00411 .00107857};*/ 00412 static double b[4]={8.5733287401,18.0590169730,8.6347608925,.2677737343}; 00413 static double c[4]={9.5733223454,25.6329561486,21.0996530827,3.9584969228}; 00414 00415 DEBUG_ENTRY( "ee1_safe()" ); 00416 00417 ASSERT( x > 1. ); 00418 00419 /* abs. accuracy better than 2e-8 */ 00420 /* top = powi(x,4) + b[0]*powi(x,3) + b[1]*x*x + b[2]*x + b[3]; */ 00421 top = (((x + b[0])*x + b[1])*x + b[2])*x + b[3]; 00422 /* bot = powi(x,4) + c[0]*powi(x,3) + c[1]*x*x + c[2]*x + c[3]; */ 00423 bot = (((x + c[0])*x + c[1])*x + c[2])*x + c[3]; 00424 00425 ans = top/bot/x; 00426 return ans; 00427 } 00428 00429 /*FFmtRead scan input line for free format number */ 00430 double FFmtRead(const char *chCard, 00431 long int *ipnt, 00432 /* the contents of this array element is the last that will be read */ 00433 long int last, 00434 bool *lgEOL) 00435 { 00436 DEBUG_ENTRY( "FFmtRead()" ); 00437 00438 const int MAX_SIZE = 1001; 00439 char chr = '\0', chLocal[MAX_SIZE]; 00440 const char *eol_ptr = &chCard[last]; // eol_ptr points one beyond last valid char 00441 const char *ptr = min(&chCard[*ipnt-1],eol_ptr); // ipnt is on fortran scale 00442 00443 ASSERT( *ipnt > 0 && last - *ipnt + 2 <= MAX_SIZE ); 00444 00445 while( ptr < eol_ptr && ( chr = *ptr++ ) != '\0' ) 00446 { 00447 const char *lptr = ptr; 00448 char lchr = chr; 00449 if( lchr == '-' || lchr == '+' ) 00450 lchr = *lptr++; 00451 if( lchr == '.' ) 00452 lchr = *lptr; 00453 if( isdigit(lchr) ) 00454 break; 00455 } 00456 00457 if( ptr == eol_ptr || chr == '\0' ) 00458 { 00459 *ipnt = last+1; 00460 *lgEOL = true; 00461 return 0.; 00462 } 00463 00464 // stripping commas is deprecated, this loop should be deleted in due course 00465 char* ptr2 = chLocal; 00466 do 00467 { 00468 if( chr != ',' ) 00469 *ptr2++ = chr; 00470 # if 0 00471 else 00472 { 00473 if( ptr != eol_ptr ) 00474 { 00475 /* don't complain about comma if it appears after number, 00476 * as determined by space following comma */ 00477 char tmpChar = *ptr; 00478 } 00479 } 00480 # endif 00481 if( ptr == eol_ptr ) 00482 break; 00483 chr = *ptr++; 00484 } 00485 while( isdigit(chr) || chr == '.' || chr == '-' || chr == '+' || chr == ',' || chr == 'e' || chr == 'E' ); 00486 *ptr2 = '\0'; 00487 00488 /*if( lgCommaFound ) 00489 fprintf( ioQQQ, " PROBLEM - a comma was found embedded in a number, this is deprecated.\n" );*/ 00490 00491 double value = atof( chLocal ); 00492 00493 *ipnt = (long)(ptr - chCard); // ptr already points 1 beyond where next read should start 00494 *lgEOL = false; 00495 return value; 00496 } 00497 00498 /*nMatch determine whether match to a keyword occurs on command line, 00499 * return value is 0 if no match, and position of match within string if hit */ 00500 long nMatch(const char *chKey, 00501 const char *chCard) 00502 { 00503 const char *ptr; 00504 long Match_v; 00505 00506 DEBUG_ENTRY( "nMatch()" ); 00507 00508 ASSERT( strlen(chKey) > 0 ); 00509 00510 if( ( ptr = strstr( chCard, chKey ) ) == NULL ) 00511 { 00512 /* did not find match, return 0 */ 00513 Match_v = 0L; 00514 } 00515 else 00516 { 00517 /* return position within chCard (fortran scale) */ 00518 Match_v = (long)(ptr-chCard+1); 00519 } 00520 return Match_v; 00521 } 00522 00523 /* fudge enter fudge factors, or some arbitrary number, with fudge command 00524 * other sections of the code access these numbers by calling fudge 00525 * fudge(0) returns the first number that was entered 00526 * prototype for this function is in cddefines.h so that it can be used without 00527 * declarations 00528 * fudge(-1) queries the routine for the number of fudge parameters that were entered, 00529 * zero returned if none */ 00530 double fudge(long int ipnt) 00531 { 00532 double fudge_v; 00533 00534 DEBUG_ENTRY( "fudge()" ); 00535 00536 if( ipnt < 0 ) 00537 { 00538 /* this is special case, return number of arguments */ 00539 fudge_v = fudgec.nfudge; 00540 fudgec.lgFudgeUsed = true; 00541 } 00542 else if( ipnt >= fudgec.nfudge ) 00543 { 00544 fprintf( ioQQQ, " FUDGE factor not entered for array number %3ld\n", 00545 ipnt ); 00546 cdEXIT(EXIT_FAILURE); 00547 } 00548 else 00549 { 00550 fudge_v = fudgec.fudgea[ipnt]; 00551 fudgec.lgFudgeUsed = true; 00552 } 00553 return fudge_v; 00554 } 00555 00556 /*GetElem scans line image, finds element. returns atomic number j, 00557 * on C scale, -1 if no hit. chCARD_CAPS must be in CAPS to hit element */ 00558 long int GetElem(char *chCARD_CAPS ) 00559 { 00560 int i; 00561 00562 DEBUG_ENTRY( "GetElem()" ); 00563 00564 /* find which element */ 00565 00566 /* >>>chng 99 apr 17, lower limit to loop had been 1, so search started with helium, 00567 * change to 0 so we can pick up hydrogen. needed for parseasserts command */ 00568 /* find match with element name, start with helium */ 00569 for( i=0; i<(int)LIMELM; ++i ) 00570 { 00571 if( nMatch( elementnames.chElementNameShort[i], chCARD_CAPS ) ) 00572 { 00573 /* return value is in C counting, hydrogen would be 0*/ 00574 return i; 00575 } 00576 } 00577 /* fall through, did not hit, return -1 as error condition */ 00578 return (-1 ); 00579 } 00580 00581 /* GetQuote get any name between double quotes off command line 00582 * return string as chLabel, is null terminated 00583 * returns zero for success, 1 for did not find double quotes */ 00584 int GetQuote( 00585 /* we will generate a label and stuff it here */ 00586 char *chLabel, 00587 /* local capd line, we will blank this out */ 00588 char *chCard , 00589 /* if true then abort if no double quotes found, 00590 * if false then return null string in this case */ 00591 bool lgABORT ) 00592 { 00593 char *i0, /* pointer to first " */ 00594 *i1, /* pointer to second ", name is in between */ 00595 *iLoc; /* pointer to first " in local version of card in calling routine */ 00596 size_t len; 00597 00598 DEBUG_ENTRY( "GetQuote()" ); 00599 00600 /*label within quotes returned within chLabel 00601 *label in line image is set to blanks when complete */ 00602 00603 /* find start of string, string must begin and end with double quotes */ 00604 /* get pointer to first quote */ 00605 i0 = strchr( input.chOrgCard,'\"' ); 00606 00607 if( i0 != NULL ) 00608 { 00609 /* get pointer to next quote */ 00610 /*i1 = strrchr( input.chOrgCard,'\"' );*/ 00611 i1 = strchr( i0+1 ,'\"' ); 00612 } 00613 else 00614 { 00615 i1 = NULL; 00616 } 00617 00618 /* check that pointers were not NULL */ 00619 /* >>chng 00 apr 27, check for i0 and i1 equal not needed anymore, by PvH */ 00620 if( i0 == NULL || i1 == NULL ) 00621 { 00622 if( lgABORT ) 00623 { 00624 /* this case really do need this to work - it did not, so must abort */ 00625 fprintf( ioQQQ, 00626 " A filename or label must be specified within double quotes, but no quotes were encountered on this command.\n" ); 00627 fprintf( ioQQQ, " Name must be surrounded by exactly two double quotes, like \"name.txt\". \n" ); 00628 fprintf( ioQQQ, " The line image follows:\n" ); 00629 fprintf( ioQQQ, " %s\n", input.chOrgCard); 00630 fprintf( ioQQQ, " Sorry\n" ); 00631 cdEXIT(EXIT_FAILURE); 00632 } 00633 else 00634 { 00635 /* this branch, ok if not present, return null string in that case */ 00636 chLabel[0] = '\0'; 00637 /* return value of 1 indicates did not find double quotes */ 00638 return 1; 00639 } 00640 } 00641 00642 /* now copy the text in between quotes */ 00643 len = (size_t)(i1-i0-1); 00644 strncpy(chLabel,i0+1,len); 00645 /* strncpy doesn't terminate the label */ 00646 chLabel[len] = '\0'; 00647 00648 /* get pointer to first quote in local copy of line image in calling routine */ 00649 iLoc = strchr( chCard ,'\"' ); 00650 if( iLoc == NULL ) 00651 { 00652 fprintf( ioQQQ, " Insanity in GetQuote - line image follows:\n" ); 00653 fprintf( ioQQQ, " %s\n", input.chOrgCard); 00654 cdEXIT(EXIT_FAILURE); 00655 } 00656 00657 /* >>chng 97 jul 01, blank out label once finished, to not be picked up later */ 00658 /* >>chng 00 apr 27, erase quotes as well, so that we can find second label, by PvH */ 00659 while( i0 <= i1 ) 00660 { 00661 *i0 = ' '; 00662 *iLoc = ' '; 00663 ++i0; 00664 ++iLoc; 00665 } 00666 /* return condition of 0 indicates success */ 00667 return 0; 00668 } 00669 00670 /* powi.c - calc x^n, where n is an integer! */ 00671 00672 /* Very slightly modified version of power() from Computer Language, Sept. 86, 00673 pg 87, by Jon Snader (who extracted the binary algorithm from Donald Knuth, 00674 "The Art of Computer Programming", vol 2, 1969). 00675 powi() will only be called when an exponentiation with an integer 00676 exponent is performed, thus tests & code for fractional exponents were 00677 removed. 00678 */ 00679 00680 /* want to define this only if no native os support exists */ 00681 #if !HAVE_POWI 00682 00683 double powi( double x, long int n ) /* returns: x^n */ 00684 /* x; base */ 00685 /* n; exponent */ 00686 { 00687 double p; /* holds partial product */ 00688 00689 DEBUG_ENTRY( "powi()" ); 00690 00691 if( x == 0 ) 00692 return 0.; 00693 00694 /* test for negative exponent */ 00695 if( n < 0 ) 00696 { 00697 n = -n; 00698 x = 1/x; 00699 } 00700 00701 p = is_odd(n) ? x : 1; /* test & set zero power */ 00702 00703 /*lint -e704 shift right of signed quantity */ 00704 /*lint -e720 Boolean test of assignment */ 00705 while( n >>= 1 ) 00706 { /* now do the other powers */ 00707 x *= x; /* sq previous power of x */ 00708 if( is_odd(n) ) /* if low order bit set */ 00709 p *= x; /* then, multiply partial product by latest power of x */ 00710 } 00711 /*lint +e704 shift right of signed quantity */ 00712 /*lint +e720 Boolean test of assignment */ 00713 return p; 00714 } 00715 00716 #endif /* HAVE_POWI */ 00717 00718 long ipow( long m, long n ) /* returns: m^n */ 00719 /* m; base */ 00720 /* n; exponent */ 00721 { 00722 long p; /* holds partial product */ 00723 00724 DEBUG_ENTRY( "ipow()" ); 00725 00726 if( m == 0 || (n < 0 && m > 1) ) 00727 return 0L; 00728 /* NOTE: negative exponent always results in 0 for integers! 00729 * (except for the case when m==1 or -1) */ 00730 00731 if( n < 0 ) 00732 { /* test for negative exponent */ 00733 n = -n; 00734 m = 1/m; 00735 } 00736 00737 p = is_odd(n) ? m : 1; /* test & set zero power */ 00738 00739 /*lint -e704 shift right of signed quantity */ 00740 /*lint -e720 Boolean test of assignment */ 00741 while( n >>= 1 ) 00742 { /* now do the other powers */ 00743 m *= m; /* sq previous power of m */ 00744 if( is_odd(n) ) /* if low order bit set */ 00745 p *= m; /* then, multiply partial product by latest power of m */ 00746 } 00747 /*lint +e704 shift right of signed quantity */ 00748 /*lint +e720 Boolean test of assignment */ 00749 return p; 00750 } 00751 00752 /*PrintE82 - series of routines to mimic 1p, e8.2 fortran output */ 00753 /*********************************************************** 00754 * contains the following sets of routines to get around * 00755 * the MS C++ compilers unusual exponential output. * 00756 * PrintEfmt <= much faster, no overhead with unix * 00757 * PrintE93 * 00758 * PrintE82 * 00759 * PrintE71 * 00760 **********************************************************/ 00761 00762 /**********************************************************/ 00763 /* 00764 * Instead of printf'ing with %e or %.5e or whatever, call 00765 * efmt("%e", val) and print the result with %s. This lets 00766 * us work around bugs in VS C 6.0. 00767 */ 00768 char *PrintEfmt(const char *fmt, double val /* , char *buf */) 00769 { 00770 static char buf[30]; /* or pass it in */ 00771 00772 DEBUG_ENTRY( "PrintEfmt()" ); 00773 00774 /* create the string */ 00775 sprintf(buf, fmt, val); 00776 00777 /* we need to fix e in format if ms vs */ 00778 # ifdef _MSC_VER 00779 { 00780 /* code to fix incorrect ms v e format. works only for case where there is 00781 * a leading space in the format - for formats like 7.1, 8.2, 9.3, 10.4, etc 00782 * result will have 1 too many characters */ 00783 char *ep , buf2[30]; 00784 00785 /* msvc behaves badly in different ways for positive vs negative sign vals, 00786 * if val is positive must create a leading space */ 00787 if( val >= 0.) 00788 { 00789 strcpy(buf2 , " " ); 00790 strcat(buf2 , buf); 00791 strcpy( buf , buf2); 00792 } 00793 00794 /* allow for both e and E formats */ 00795 if((ep = strchr(buf, 'e')) == NULL) 00796 { 00797 ep = strchr(buf, 'E'); 00798 } 00799 00800 /* ep can still be NULL if val is Inf or NaN */ 00801 if(ep != NULL) 00802 { 00803 /* move pointer two char past the e, to pick up the e and sign */ 00804 ep += 2; 00805 00806 /* terminate buf where the e is, *ep points to this location */ 00807 *ep = '\0'; 00808 00809 /* skip next char, */ 00810 ++ep; 00811 00812 /* copy resulting string to return string */ 00813 strcat( buf, ep ); 00814 } 00815 } 00816 # endif 00817 return buf; 00818 } 00819 00820 /**********************************************************/ 00821 void PrintE82( FILE* ioOUT, double value ) 00822 { 00823 double frac , xlog , xfloor , tvalue; 00824 int iExp; 00825 00826 DEBUG_ENTRY( "PrintE82()" ); 00827 00828 if( value < 0 ) 00829 { 00830 fprintf(ioOUT,"********"); 00831 } 00832 else if( value == 0 ) 00833 { 00834 fprintf(ioOUT,"0.00E+00"); 00835 } 00836 else 00837 { 00838 /* round number off for 8.2 format (not needed since can't be negative) */ 00839 tvalue = value; 00840 xlog = log10( tvalue ); 00841 xfloor = floor(xlog); 00842 /* this is now the fractional part */ 00843 frac = tvalue*pow(10.,-xfloor); 00844 /*this is the possibly signed exponential part */ 00845 iExp = (int)xfloor; 00846 if( frac>9.9945 ) 00847 { 00848 frac /= 10.; 00849 iExp += 1; 00850 } 00851 /* print the fractional part*/ 00852 fprintf(ioOUT,"%.2f",frac); 00853 /* E for exponent */ 00854 fprintf(ioOUT,"E"); 00855 /* if positive throw a + sign*/ 00856 if(iExp>=0 ) 00857 { 00858 fprintf(ioOUT,"+"); 00859 } 00860 fprintf(ioOUT,"%.2d",iExp); 00861 } 00862 return; 00863 } 00864 /* 00865 *============================================================================== 00866 */ 00867 void PrintE71( FILE* ioOUT, double value ) 00868 { 00869 double frac , xlog , xfloor , tvalue; 00870 int iExp; 00871 00872 DEBUG_ENTRY( "PrintE71()" ); 00873 00874 if( value < 0 ) 00875 { 00876 fprintf(ioOUT,"*******"); 00877 } 00878 else if( value == 0 ) 00879 { 00880 fprintf(ioOUT,"0.0E+00"); 00881 } 00882 else 00883 { 00884 /* round number off for 8.2 format (not needed since can't be negative) */ 00885 tvalue = value; 00886 xlog = log10( tvalue ); 00887 xfloor = floor(xlog); 00888 /* this is now the fractional part */ 00889 frac = tvalue*pow(10.,-xfloor); 00890 /*this is the possibly signed exponential part */ 00891 iExp = (int)xfloor; 00892 if( frac>9.9945 ) 00893 { 00894 frac /= 10.; 00895 iExp += 1; 00896 } 00897 /* print the fractional part*/ 00898 fprintf(ioOUT,"%.1f",frac); 00899 /* E for exponent */ 00900 fprintf(ioOUT,"E"); 00901 /* if positive throw a + sign*/ 00902 if(iExp>=0 ) 00903 { 00904 fprintf(ioOUT,"+"); 00905 } 00906 fprintf(ioOUT,"%.2d",iExp); 00907 } 00908 return; 00909 } 00910 00911 /* 00912 *============================================================================== 00913 */ 00914 void PrintE93( FILE* ioOUT, double value ) 00915 { 00916 double frac , xlog , xfloor, tvalue; 00917 int iExp; 00918 00919 DEBUG_ENTRY( "PrintE93()" ); 00920 00921 if( value < 0 ) 00922 { 00923 fprintf(ioOUT,"*********"); 00924 } 00925 else if( value == 0 ) 00926 { 00927 fprintf(ioOUT,"0.000E+00"); 00928 } 00929 else 00930 { 00931 /* round number off for 9.3 format, neg numb not possible */ 00932 tvalue = value; 00933 xlog = log10( tvalue ); 00934 xfloor = floor(xlog); 00935 /* this is now the fractional part */ 00936 frac = tvalue*pow(10.,-xfloor); 00937 /*this is the possibly signed exponential part */ 00938 iExp = (int)xfloor; 00939 if( frac>9.99949 ) 00940 { 00941 frac /= 10.; 00942 iExp += 1; 00943 } 00944 /* print the fractional part*/ 00945 fprintf(ioOUT,"%5.3f",frac); 00946 /* E for exponent */ 00947 fprintf(ioOUT,"E"); 00948 /* if positive throw a + sign*/ 00949 if(iExp>=0 ) 00950 { 00951 fprintf(ioOUT,"+"); 00952 } 00953 fprintf(ioOUT,"%.2d",iExp); 00954 } 00955 return; 00956 } 00957 00958 /*TotalInsanity general error handler for something that cannot happen */ 00959 NORETURN void TotalInsanity(void) 00960 { 00961 00962 DEBUG_ENTRY( "TotalInsanity()" ); 00963 00964 /* something that cannot happen, happened, 00965 * if this message is triggered, simply place a breakpoint here 00966 * and debug the error */ 00967 fprintf( ioQQQ, " Something that cannot happen, has happened.\n" ); 00968 fprintf( ioQQQ, " This is TotalInsanity, I live in service.cpp.\n" ); 00969 ShowMe(); 00970 00971 cdEXIT(EXIT_FAILURE); 00972 } 00973 00974 00975 /*BadRead general error handler for trying to read data, but failing */ 00976 NORETURN void BadRead(void) 00977 { 00978 00979 DEBUG_ENTRY( "BadRead()" ); 00980 00981 /* read failed */ 00982 fprintf( ioQQQ, " A read of internal input data has failed.\n" ); 00983 fprintf( ioQQQ, " This is BadRead, I live in service.c.\n" ); 00984 ShowMe(); 00985 00986 cdEXIT(EXIT_FAILURE); 00987 } 00988 00989 /*BadOpen general error handler for trying to open file, but failing */ 00990 NORETURN void BadOpen(void) 00991 { 00992 00993 DEBUG_ENTRY( "BadOpen()" ); 00994 00995 /* read failed */ 00996 fprintf( ioQQQ, " An attempt at opening a files has failed.\n" ); 00997 fprintf( ioQQQ, " This is BadOpen, I live in service.c.\n" ); 00998 ShowMe(); 00999 01000 cdEXIT(EXIT_FAILURE); 01001 } 01002 01003 /*NoNumb general error handler for no numbers on input line */ 01004 NORETURN void NoNumb(char *chCard) 01005 { 01006 01007 DEBUG_ENTRY( "NoNumb()" ); 01008 01009 /* general catch-all for no number when there should have been */ 01010 fprintf( ioQQQ, " There is a problem on the following command line:\n" ); 01011 fprintf( ioQQQ, " %s\n", chCard ); 01012 fprintf( ioQQQ, " There should have been a number on this line. Sorry.\n" ); 01013 cdEXIT(EXIT_FAILURE); 01014 } 01015 01016 /*sexp safe exponential function */ 01017 sys_float sexp(sys_float x) 01018 { 01019 sys_float sexp_v; 01020 01021 DEBUG_ENTRY( "sexp()" ); 01022 01023 /* SEXP_LIMIT is 84 in cddefines.h */ 01024 if( x < SEXP_LIMIT ) 01025 { 01026 sexp_v = exp(-x); 01027 } 01028 else 01029 { 01030 sexp_v = 0.f; 01031 } 01032 return sexp_v; 01033 } 01034 01035 /*sexp safe exponential function */ 01036 double sexp(double x) 01037 { 01038 double sexp_v; 01039 01040 DEBUG_ENTRY( "sexp()" ); 01041 01042 /* SEXP_LIMIT is 84 in cddefines.h */ 01043 if( x < SEXP_LIMIT ) 01044 { 01045 sexp_v = exp(-x); 01046 } 01047 else 01048 { 01049 sexp_v = 0.; 01050 } 01051 return sexp_v; 01052 } 01053 01054 01055 /*dsexp safe exponential function for doubles */ 01056 double dsexp(double x) 01057 { 01058 double dsexp_v; 01059 01060 DEBUG_ENTRY( "dsexp()" ); 01061 01062 if( x < 680. ) 01063 { 01064 dsexp_v = exp(-x); 01065 } 01066 else 01067 { 01068 dsexp_v = 0.; 01069 } 01070 return dsexp_v; 01071 } 01072 01073 /*TestCode set flag saying that test code is in place 01074 * prototype in cddefines.h */ 01075 void TestCode(void) 01076 { 01077 01078 DEBUG_ENTRY( "TestCode( )" ); 01079 01080 /* called if test code is in place */ 01081 lgTestCodeCalled = true; 01082 return; 01083 } 01084 01085 /*broken set flag saying that the code is broken, */ 01086 void broken(void) 01087 { 01088 01089 DEBUG_ENTRY( "broken( )" ); 01090 01091 broke.lgBroke = true; 01092 return; 01093 } 01094 01095 /*fixit say that code needs to be fixed */ 01096 void fixit(void) 01097 { 01098 DEBUG_ENTRY( "fixit( )" ); 01099 01100 broke.lgFixit = true; 01101 return; 01102 } 01103 01104 /*CodeReview placed next to code that needs to be checked */ 01105 void CodeReview(void) 01106 { 01107 DEBUG_ENTRY( "CodeReview( )" ); 01108 01109 broke.lgCheckit = true; 01110 return; 01111 } 01112 01114 int dprintf(FILE *fp, const char *format, ...) 01115 { 01116 va_list ap; 01117 int i1, i2; 01118 01119 DEBUG_ENTRY( "dprintf()" ); 01120 va_start(ap,format); 01121 i1 = fprintf(fp,"DEBUG "); 01122 if (i1 >= 0) 01123 i2 = vfprintf(fp,format,ap); 01124 else 01125 i2 = 0; 01126 if (i2 < 0) 01127 i1 = 0; 01128 va_end(ap); 01129 01130 return i1+i2; 01131 } 01132 01133 /* dbg_printf is a debug print routine that was provided by Peter Teuben, 01134 * as a component from his NEMO package. It offers run-time specification 01135 * of the level of debugging */ 01136 int dbg_printf(int debug, const char *fmt, ...) 01137 { 01138 va_list ap; 01139 int i=0; 01140 01141 DEBUG_ENTRY( "dbg_printf()" ); 01142 01143 /* print this debug message? (debug_level not currently used)*/ 01144 if(debug <= trace.debug_level) 01145 { 01146 va_start(ap, fmt); 01147 01148 i = vfprintf(ioQQQ, fmt, ap); 01149 /* drain ioQQQ */ 01150 fflush(ioQQQ); 01151 va_end(ap); 01152 } 01153 return i; 01154 } 01155 01156 01157 /*qg32 32 point Gaussian quadrature, originally given to Gary F by Jim Lattimer */ 01158 double qg32( 01159 double xl, /*lower limit to integration range*/ 01160 double xu, /*upper limit to integration range*/ 01161 /*following is the pointer to the function that will be evaluated*/ 01162 double (*fct)(double) ) 01163 { 01164 double a, 01165 b, 01166 c, 01167 y; 01168 01169 DEBUG_ENTRY( "qg32()" ); 01170 01171 /******************************************************************************** 01172 * * 01173 * 32-point Gaussian quadrature * 01174 * xl : the lower limit of integration * 01175 * xu : the upper limit * 01176 * fct : the (external) function * 01177 * returns the value of the integral * 01178 * * 01179 * simple call to integrate sine from 0 to pi * 01180 * double agn = qg32( 0., 3.141592654 , sin ); * 01181 * * 01182 *******************************************************************************/ 01183 01184 a = .5*(xu + xl); 01185 b = xu - xl; 01186 c = .498631930924740780*b; 01187 y = .35093050047350483e-2*((*fct)(a+c) + (*fct)(a-c)); 01188 c = b*.49280575577263417; 01189 y += .8137197365452835e-2*((*fct)(a+c) + (*fct)(a-c)); 01190 c = b*.48238112779375322; 01191 y += .1269603265463103e-1*((*fct)(a+c) + (*fct)(a-c)); 01192 c = b*.46745303796886984; 01193 y += .17136931456510717e-1*((*fct)(a+c) + (*fct)(a-c)); 01194 c = b*.44816057788302606; 01195 y += .21417949011113340e-1*((*fct)(a+c) + (*fct)(a-c)); 01196 c = b*.42468380686628499; 01197 y += .25499029631188088e-1*((*fct)(a+c) + (*fct)(a-c)); 01198 c = b*.3972418979839712; 01199 y += .29342046739267774e-1*((*fct)(a+c) + (*fct)(a-c)); 01200 c = b*.36609105937014484; 01201 y += .32911111388180923e-1*((*fct)(a+c) + (*fct)(a-c)); 01202 c = b*.3315221334651076; 01203 y += .36172897054424253e-1*((*fct)(a+c) + (*fct)(a-c)); 01204 c = b*.29385787862038116; 01205 y += .39096947893535153e-1*((*fct)(a+c) + (*fct)(a-c)); 01206 c = b*.2534499544661147; 01207 y += .41655962113473378e-1*((*fct)(a+c) + (*fct)(a-c)); 01208 c = b*.21067563806531767; 01209 y += .43826046502201906e-1*((*fct)(a+c) + (*fct)(a-c)); 01210 c = b*.16593430114106382; 01211 y += .45586939347881942e-1*((*fct)(a+c) + (*fct)(a-c)); 01212 c = b*.11964368112606854; 01213 y += .46922199540402283e-1*((*fct)(a+c) + (*fct)(a-c)); 01214 c = b*.7223598079139825e-1; 01215 y += .47819360039637430e-1*((*fct)(a+c) + (*fct)(a-c)); 01216 c = b*.24153832843869158e-1; 01217 y = b*(y + .482700442573639e-1*((*fct)(a+c) + (*fct)(a-c))); 01218 /* the answer */ 01219 return y; 01220 } 01221 01222 /*spsort netlib routine to sort array returning sorted indices */ 01223 void spsort( 01224 /* input array to be sorted */ 01225 realnum x[], 01226 /* number of values in x */ 01227 long int n, 01228 /* permutation output array */ 01229 long int iperm[], 01230 /* flag saying what to do - 1 sorts into increasing order, not changing 01231 * the original vector, -1 sorts into decreasing order. 2, -2 change vector */ 01232 int kflag, 01233 /* error condition, should be 0 */ 01234 int *ier) 01235 { 01236 /* 01237 ****BEGIN PROLOGUE SPSORT 01238 ****PURPOSE Return the permutation vector generated by sorting a given 01239 * array and, optionally, rearrange the elements of the array. 01240 * The array may be sorted in increasing or decreasing order. 01241 * A slightly modified quicksort algorithm is used. 01242 ****LIBRARY SLATEC 01243 ****CATEGORY N6A1B, N6A2B 01244 ****TYPE SINGLE PRECISION (SPSORT-S, DPSORT-D, IPSORT-I, HPSORT-H) 01245 ****KEY WORDS NUMBER SORTING, PASSIVE SORTING, SINGLETON QUICKSORT, SORT 01246 ****AUTHOR Jones, R. E., (SNLA) 01247 * Rhoads, G. S., (NBS) 01248 * Wisniewski, J. A., (SNLA) 01249 ****DESCRIPTION 01250 * 01251 * SPSORT returns the permutation vector IPERM generated by sorting 01252 * the array X and, optionally, rearranges the values in X. X may 01253 * be sorted in increasing or decreasing order. A slightly modified 01254 * quicksort algorithm is used. 01255 * 01256 * IPERM is such that X(IPERM(I)) is the Ith value in the rearrangement 01257 * of X. IPERM may be applied to another array by calling IPPERM, 01258 * SPPERM, DPPERM or HPPERM. 01259 * 01260 * The main difference between SPSORT and its active sorting equivalent 01261 * SSORT is that the data are referenced indirectly rather than 01262 * directly. Therefore, SPSORT should require approximately twice as 01263 * long to execute as SSORT. However, SPSORT is more general. 01264 * 01265 * Description of Parameters 01266 * X - input/output -- real array of values to be sorted. 01267 * If ABS(KFLAG) = 2, then the values in X will be 01268 * rearranged on output; otherwise, they are unchanged. 01269 * N - input -- number of values in array X to be sorted. 01270 * IPERM - output -- permutation array such that IPERM(I) is the 01271 * index of the value in the original order of the 01272 * X array that is in the Ith location in the sorted 01273 * order. 01274 * KFLAG - input -- control parameter: 01275 * = 2 means return the permutation vector resulting from 01276 * sorting X in increasing order and sort X also. 01277 * = 1 means return the permutation vector resulting from 01278 * sorting X in increasing order and do not sort X. 01279 * = -1 means return the permutation vector resulting from 01280 * sorting X in decreasing order and do not sort X. 01281 * = -2 means return the permutation vector resulting from 01282 * sorting X in decreasing order and sort X also. 01283 * IER - output -- error indicator: 01284 * = 0 if no error, 01285 * = 1 if N is zero or negative, 01286 * = 2 if KFLAG is not 2, 1, -1, or -2. 01287 ****REFERENCES R. C. Singleton, Algorithm 347, An efficient algorithm 01288 * for sorting with minimal storage, Communications of 01289 * the ACM, 12, 3 (1969), pp. 185-187. 01290 ****ROUTINES CALLED XERMSG 01291 ****REVISION HISTORY (YYMMDD) 01292 * 761101 DATE WRITTEN 01293 * 761118 Modified by John A. Wisniewski to use the Singleton 01294 * quicksort algorithm. 01295 * 870423 Modified by Gregory S. Rhoads for passive sorting with the 01296 * option for the rearrangement of the original data. 01297 * 890620 Algorithm for rearranging the data vector corrected by R. 01298 * Boisvert. 01299 * 890622 Prologue upgraded to Version 4.0 style by D. Lozier. 01300 * 891128 Error when KFLAG.LT.0 and N=1 corrected by R. Boisvert. 01301 * 920507 Modified by M. McClain to revise prologue text. 01302 * 920818 Declarations section rebuilt and code restructured to use 01303 * IF-THEN-ELSE-ENDIF. (SMR, WRB) 01304 ****END PROLOGUE SPSORT 01305 * .. Scalar Arguments .. 01306 */ 01307 long int i, 01308 ij, 01309 il[21], 01310 indx, 01311 indx0, 01312 istrt, 01313 istrt_, 01314 iu[21], 01315 j, 01316 k, 01317 kk, 01318 l, 01319 lm, 01320 lmt, 01321 m, 01322 nn; 01323 realnum r, 01324 ttemp; 01325 01326 DEBUG_ENTRY( "spsort()" ); 01327 01328 /* .. Array Arguments .. */ 01329 /* .. Local Scalars .. */ 01330 /* .. Local Arrays .. */ 01331 /* .. External Subroutines .. */ 01332 /* .. Intrinsic Functions .. */ 01333 /****FIRST EXECUTABLE STATEMENT SPSORT */ 01334 *ier = 0; 01335 nn = n; 01336 if( nn < 1 ) 01337 { 01338 *ier = 1; 01339 return; 01340 } 01341 else 01342 { 01343 kk = labs(kflag); 01344 if( kk != 1 && kk != 2 ) 01345 { 01346 *ier = 2; 01347 return; 01348 } 01349 else 01350 { 01351 01352 /* Initialize permutation vector to index on f scale 01353 * */ 01354 for( i=0; i < nn; i++ ) 01355 { 01356 iperm[i] = i+1; 01357 } 01358 01359 /* Return if only one value is to be sorted */ 01360 if( nn == 1 ) 01361 { 01362 --iperm[0]; 01363 return; 01364 } 01365 01366 /* Alter array X to get decreasing order if needed */ 01367 if( kflag <= -1 ) 01368 { 01369 for( i=0; i < nn; i++ ) 01370 { 01371 x[i] = -x[i]; 01372 } 01373 } 01374 01375 /* Sort X only */ 01376 m = 1; 01377 i = 1; 01378 j = nn; 01379 r = .375e0; 01380 } 01381 } 01382 01383 while( true ) 01384 { 01385 if( i == j ) 01386 goto L_80; 01387 if( r <= 0.5898437e0 ) 01388 { 01389 r += 3.90625e-2; 01390 } 01391 else 01392 { 01393 r -= 0.21875e0; 01394 } 01395 01396 L_40: 01397 k = i; 01398 01399 /* Select a central element of the array and save it in location L 01400 * */ 01401 ij = i + (long)((j-i)*r); 01402 lm = iperm[ij-1]; 01403 01404 /* If first element of array is greater than LM, interchange with LM 01405 * */ 01406 if( x[iperm[i-1]-1] > x[lm-1] ) 01407 { 01408 iperm[ij-1] = iperm[i-1]; 01409 iperm[i-1] = lm; 01410 lm = iperm[ij-1]; 01411 } 01412 l = j; 01413 01414 /* If last element of array is less than LM, interchange with LM 01415 * */ 01416 if( x[iperm[j-1]-1] < x[lm-1] ) 01417 { 01418 iperm[ij-1] = iperm[j-1]; 01419 iperm[j-1] = lm; 01420 lm = iperm[ij-1]; 01421 01422 /* If first element of array is greater than LM, interchange 01423 * with LM 01424 * */ 01425 if( x[iperm[i-1]-1] > x[lm-1] ) 01426 { 01427 iperm[ij-1] = iperm[i-1]; 01428 iperm[i-1] = lm; 01429 lm = iperm[ij-1]; 01430 } 01431 } 01432 01433 /* Find an element in the second half of the array which is smaller 01434 * than LM */ 01435 while( true ) 01436 { 01437 l -= 1; 01438 if( x[iperm[l-1]-1] <= x[lm-1] ) 01439 { 01440 01441 /* Find an element in the first half of the array which is greater 01442 * than LM */ 01443 while( true ) 01444 { 01445 k += 1; 01446 if( x[iperm[k-1]-1] >= x[lm-1] ) 01447 break; 01448 } 01449 01450 /* Interchange these elements */ 01451 if( k > l ) 01452 break; 01453 lmt = iperm[l-1]; 01454 iperm[l-1] = iperm[k-1]; 01455 iperm[k-1] = lmt; 01456 } 01457 } 01458 01459 /* Save upper and lower subscripts of the array yet to be sorted */ 01460 if( l - i > j - k ) 01461 { 01462 il[m-1] = i; 01463 iu[m-1] = l; 01464 i = k; 01465 m += 1; 01466 } 01467 else 01468 { 01469 il[m-1] = k; 01470 iu[m-1] = j; 01471 j = l; 01472 m += 1; 01473 } 01474 01475 L_90: 01476 if( j - i >= 1 ) 01477 goto L_40; 01478 if( i == 1 ) 01479 continue; 01480 i -= 1; 01481 01482 while( true ) 01483 { 01484 i += 1; 01485 if( i == j ) 01486 break; 01487 lm = iperm[i]; 01488 if( x[iperm[i-1]-1] > x[lm-1] ) 01489 { 01490 k = i; 01491 01492 while( true ) 01493 { 01494 iperm[k] = iperm[k-1]; 01495 k -= 1; 01496 01497 if( x[lm-1] >= x[iperm[k-1]-1] ) 01498 break; 01499 } 01500 iperm[k] = lm; 01501 } 01502 } 01503 01504 /* Begin again on another portion of the unsorted array */ 01505 L_80: 01506 m -= 1; 01507 if( m == 0 ) 01508 break; 01509 /*lint -e644 not explicitly initialized */ 01510 i = il[m-1]; 01511 j = iu[m-1]; 01512 /*lint +e644 not explicitly initialized */ 01513 goto L_90; 01514 } 01515 01516 /* Clean up */ 01517 if( kflag <= -1 ) 01518 { 01519 for( i=0; i < nn; i++ ) 01520 { 01521 x[i] = -x[i]; 01522 } 01523 } 01524 01525 /* Rearrange the values of X if desired */ 01526 if( kk == 2 ) 01527 { 01528 01529 /* Use the IPERM vector as a flag. 01530 * If IPERM(I) < 0, then the I-th value is in correct location */ 01531 for( istrt=1; istrt <= nn; istrt++ ) 01532 { 01533 istrt_ = istrt - 1; 01534 if( iperm[istrt_] >= 0 ) 01535 { 01536 indx = istrt; 01537 indx0 = indx; 01538 ttemp = x[istrt_]; 01539 while( iperm[indx-1] > 0 ) 01540 { 01541 x[indx-1] = x[iperm[indx-1]-1]; 01542 indx0 = indx; 01543 iperm[indx-1] = -iperm[indx-1]; 01544 indx = labs(iperm[indx-1]); 01545 } 01546 x[indx0-1] = ttemp; 01547 } 01548 } 01549 01550 /* Revert the signs of the IPERM values */ 01551 for( i=0; i < nn; i++ ) 01552 { 01553 iperm[i] = -iperm[i]; 01554 } 01555 } 01556 01557 for( i=0; i < nn; i++ ) 01558 { 01559 --iperm[i]; 01560 } 01561 return; 01562 } 01563 01564 /*MyMalloc wrapper for malloc(). Returns a good pointer or dies. 01565 * memory is filled with NaN 01566 * >>chng 05 dec 14, do not set to NaN since tricks debugger 01567 * routines within code do not call this or malloc, but rather MALLOC 01568 * which is resolved into MyMalloc or malloc depending on whether 01569 * NDEBUG is set by the compiler to indicate "not debugging", 01570 * in typical negative C style */ 01571 void *MyMalloc( 01572 /*use same type as library function MALLOC*/ 01573 size_t size , 01574 const char *chFile, 01575 int line 01576 ) 01577 { 01578 void *ptr; 01579 01580 DEBUG_ENTRY( "MyMalloc()" ); 01581 01582 ASSERT( size > 0 ); 01583 01584 /* debug branch for printing malloc args */ 01585 { 01586 /*@-redef@*/ 01587 enum{DEBUG_LOC=false}; 01588 /*@+redef@*/ 01589 if( DEBUG_LOC) 01590 { 01591 static long int kount=0, nTot=0; 01592 nTot += (long)size; 01593 fprintf(ioQQQ,"%li\t%li\t%li\n", 01594 kount , 01595 (long)size , 01596 nTot ); 01597 ++kount; 01598 } 01599 } 01600 01601 if( ( ptr = malloc( size ) ) == NULL ) 01602 { 01603 fprintf(ioQQQ,"DISASTER MyMalloc could not allocate %lu bytes. Exit in MyMalloc.", 01604 (unsigned long)size ); 01605 fprintf(ioQQQ,"MyMalloc called from file %s at line %i.\n", 01606 chFile , line ); 01607 cdEXIT(EXIT_FAILURE); 01608 } 01609 01610 /* flag -DNOINIT will turn off this initialization which can fool valgrind/purify */ 01611 # if !defined(NDEBUG) && !defined(NOINIT) 01612 01613 size_t nFloat = size/4; 01614 size_t nDouble = size/8; 01615 sys_float *fptr = static_cast<sys_float*>(ptr); 01616 double *dptr = static_cast<double*>(ptr); 01617 01618 /* >>chng 04 feb 03, fill memory with invalid numbers, PvH */ 01619 /* on IA32/AMD64 processors this will generate NaN's for both float and double; 01620 * on most other (modern) architectures it is likely to do the same... */ 01621 /* >>chng 05 dec 14, change code to generate signaling NaN's for most cases (but not all!) */ 01622 if( size == nDouble*8 ) 01623 { 01624 /* this could be an array of doubles as well as floats -> we will hedge our bets 01625 * we will fill the array with a pattern that is interpreted as all signaling 01626 * NaN's for doubles, and alternating signaling and quiet NaN's for floats: 01627 * byte offset: 0 4 8 12 16 01628 * double | SNaN | SNan | 01629 * float | SNaN | QNaN | SNan | QNaN | (little-endian, e.g. Intel, AMD, alpha) 01630 * float | QNaN | SNaN | QNan | SNaN | (big-endian, e.g. Sparc, PowerPC, MIPS) */ 01631 set_NaN( dptr, (long)nDouble ); 01632 } 01633 else if( size == nFloat*4 ) 01634 { 01635 /* this could be an arrays of floats, but not doubles -> init to all float SNaN */ 01636 set_NaN( fptr, (long)nFloat ); 01637 } 01638 else 01639 { 01640 memset( ptr, 0xff, size ); 01641 } 01642 01643 # endif /* !defined(NDEBUG) && !defined(NOINIT) */ 01644 return ptr; 01645 } 01646 01647 01648 /* wrapper for calloc(). Returns a good pointer or dies. 01649 * routines within code do not call this or malloc, but rather CALLOC 01650 * which is resolved into MyCalloc or calloc depending on whether 01651 * NDEBUG is set in cddefines. \h */ 01652 void *MyCalloc( 01653 /*use same type as library function CALLOC*/ 01654 size_t num , 01655 size_t size ) 01656 { 01657 void *ptr; 01658 01659 DEBUG_ENTRY( "MyCalloc()" ); 01660 01661 ASSERT( size > 0 ); 01662 01663 /* debug branch for printing malloc args */ 01664 { 01665 /*@-redef@*/ 01666 enum{DEBUG_LOC=false}; 01667 /*@+redef@*/ 01668 if( DEBUG_LOC) 01669 { 01670 static long int kount=0; 01671 fprintf(ioQQQ,"%li\tcall\t%li\tbytes\n", kount, 01672 (long)size ); 01673 ++kount; 01674 } 01675 } 01676 01677 if( ( ptr = calloc( num , size ) ) == NULL ) 01678 { 01679 fprintf(ioQQQ,"MyCalloc could not allocate %lu bytes. Exit in MyCalloc.", 01680 (unsigned long)size ); 01681 cdEXIT(EXIT_FAILURE); 01682 } 01683 return ptr; 01684 } 01685 01686 /* wrapper for realloc(). Returns a good pointer or dies. 01687 * routines within code do not call this or malloc, but rather REALLOC 01688 * which is resolved into MyRealloc or realloc depending on whether 01689 * NDEBUG is set in cddefines.h */ 01690 void *MyRealloc( 01691 /*use same type as library function realloc */ 01692 void *p , 01693 size_t size ) 01694 { 01695 void *ptr; 01696 01697 DEBUG_ENTRY( "MyRealloc()" ); 01698 01699 ASSERT( size > 0 ); 01700 01701 /* debug branch for printing malloc args */ 01702 { 01703 /*@-redef@*/ 01704 enum{DEBUG_LOC=false}; 01705 /*@+redef@*/ 01706 if( DEBUG_LOC) 01707 { 01708 static long int kount=0; 01709 fprintf(ioQQQ,"%li\tcall\t%li\tbytes\n", kount, 01710 (long)size ); 01711 ++kount; 01712 } 01713 } 01714 01715 if( ( ptr = realloc( p , size ) ) == NULL ) 01716 { 01717 fprintf(ioQQQ,"MyRealloc could not allocate %lu bytes. Exit in MyRealloc.", 01718 (unsigned long)size ); 01719 cdEXIT(EXIT_FAILURE); 01720 } 01721 return ptr; 01722 } 01723 01724 /* function to facilitate addressing opacity array */ 01725 double csphot( 01726 /* INU is array index pointing to frequency where opacity is to be evaluated 01727 * on f not c scale */ 01728 long int inu, 01729 /* ITHR is pointer to threshold*/ 01730 long int ithr, 01731 /* IOFSET is offset as defined in opac0*/ 01732 long int iofset) 01733 { 01734 double csphot_v; 01735 01736 DEBUG_ENTRY( "csphot()" ); 01737 01738 csphot_v = opac.OpacStack[inu-ithr+iofset-1]; 01739 return csphot_v; 01740 } 01741 01742 /*RandGauss normal Gaussian random number generator 01743 * the user must call srand to set the seed before using this routine. 01744 01745 * the random numbers will have a mean of xMean and a standard deviation of s 01746 01747 * The convention is for srand to be called when the command setting 01748 * the noise is parsed 01749 01750 * for very small dispersion there are no issues, but when the dispersion becomes 01751 * large the routine will find negative values - so most often used in this case 01752 * to find dispersion in log space 01753 * this routine will return a normal Gaussian - must be careful in how this is 01754 * used when adding noise to physical quantity */ 01755 /* 01756 NB - following from Ryan Porter: 01757 I discovered that I unintentionally created an antisymmetric skew in my 01758 Monte Carlo. RandGauss is symmetric in log space, which means it is not 01759 symmetric in linear space. But to get the right standard deviation you 01760 have to take 10^x, where x is the return from RandGauss. The problem is 01761 10^x will happen less frequently than 10^-x, so without realizing it, the 01762 average "tweak" to every piece of atomic data in my Monte Carlo run was 01763 not 1.0 but something greater than 1.0, causing every single line to have 01764 an average Monte Carlo emissivity greater than the regular value. Any place 01765 that takes 10^RandGauss() needs to be adjusted if what is intended is +/- x. */ 01766 double RandGauss( 01767 /* mean value */ 01768 double xMean, 01769 /*standard deviation s */ 01770 double s ) 01771 { 01772 double x1, x2, w, yy1; 01773 static double yy2=-BIGDOUBLE; 01774 static int use_last = false; 01775 01776 DEBUG_ENTRY( "RandGauss()" ); 01777 01778 if( use_last ) 01779 { 01780 yy1 = yy2; 01781 use_last = false; 01782 } 01783 else 01784 { 01785 do { 01786 x1 = 2.*genrand_real3() - 1.; 01787 x2 = 2.*genrand_real3() - 1.; 01788 w = x1 * x1 + x2 * x2; 01789 } while ( w >= 1.0 ); 01790 01791 w = sqrt((-2.0*log(w))/w); 01792 yy1 = x1 * w; 01793 yy2 = x2 * w; 01794 use_last = true; 01795 } 01796 return xMean + yy1 * s; 01797 } 01798 01799 /* MyGaussRand takes as input a percent uncertainty less than 50% 01800 * (expressed as 0.5). The routine then assumes this input variable represents two 01801 * standard deviations about a mean of unity, and returns a random number within 01802 * that range. A hard cutoff is imposed at the two standard deviations, which 01803 * eliminates roughly 2% of the normal distribution. In other words, the routine 01804 * returns a number in a normal distribution with standard deviation equal to 01805 * half of the input, and returns a number between 1-2*stdev and 1+2*stdev. */ 01806 double MyGaussRand( double PctUncertainty ) 01807 { 01808 double StdDev; 01809 double result; 01810 01811 DEBUG_ENTRY( "MyGaussRand()" ); 01812 01813 ASSERT( PctUncertainty < 0.5 ); 01814 /* We want this "percent uncertainty to represent two standard deviations */ 01815 StdDev = 0.5 * PctUncertainty; 01816 01817 do 01818 { 01819 /*result = pow( 10., RandGauss( 0., logStdDev ) );*/ 01820 result = 1.+RandGauss( 0., StdDev ); 01821 } 01822 /* This will give us a result that is less than or equal to the 01823 * percent uncertainty about 98% of the time. */ 01824 while( (result < 1.-PctUncertainty) || (result > 1+PctUncertainty) ); 01825 01826 ASSERT( result>0. && result<2. ); 01827 return result; 01828 } 01829 01830 /*plankf evaluate Planck function for any cell at current electron temperature */ 01831 double plankf(long int ip) 01832 { 01833 double plankf_v; 01834 01835 DEBUG_ENTRY( "plankf()" ); 01836 01837 /* evaluate Planck function 01838 * argument is pointer to cell energy in ANU 01839 * return photon flux for cell IP */ 01840 if( rfield.ContBoltz[ip] <= 0. ) 01841 { 01842 plankf_v = 1e-35; 01843 } 01844 else 01845 { 01846 plankf_v = 6.991e-21*POW2(FR1RYD*rfield.anu[ip])/ 01847 (1./rfield.ContBoltz[ip] - 1.)*FR1RYD*4.; 01848 } 01849 return plankf_v; 01850 }