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 /*InitAssertResults, this must be called first, done at startup of ParseCommands*/ 00004 /*ParseAssertResults - parse input stream */ 00005 /*lgCheckAsserts, checks asserts, last thing cloudy calls, returns true if all are ok, false if problems */ 00006 #include "cddefines.h" 00007 #include "input.h" 00008 #include "conv.h" 00009 #include "optimize.h" 00010 #include "iso.h" 00011 #include "called.h" 00012 #include "atmdat.h" 00013 #include "hcmap.h" 00014 #include "thermal.h" 00015 #include "pressure.h" 00016 #include "struc.h" 00017 #include "wind.h" 00018 #include "h2.h" 00019 #include "colden.h" 00020 #include "dense.h" 00021 #include "lines.h" 00022 #include "secondaries.h" 00023 #include "radius.h" 00024 #include "version.h" 00025 #include "hmi.h" 00026 #include "prt.h" 00027 #include "grainvar.h" 00028 #include "atomfeii.h" 00029 #include "cddrive.h" 00030 #include "elementnames.h" 00031 #include "assertresults.h" 00032 #include "taulines.h" 00033 #include "lines_service.h" 00034 00035 /* flag to remember that InitAssertResults was called */ 00036 static bool lgInitDone=false , 00037 /* will be set true when space for asserts is allocated */ 00038 lgSpaceAllocated=false; 00039 00040 /* number of asserts we can handle, used in dim of space */ 00041 #define NASSERTS 100 00042 00043 /* default relative error for asserted quantities */ 00044 #define DEF_ERROR 0.05 00045 00046 /* dim of 5 also appears in MALLOC below */ 00047 #define NCHAR 5 00048 static char **chAssertLineLabel; 00049 00050 /* this will be = for equal, < or > for limit */ 00051 static char *chAssertLimit; 00052 00053 /* this will be a two character label identifying which type of assert */ 00054 static char **chAssertType; 00055 00056 /* the values and error in the asserted quantity */ 00057 static double *AssertQuantity ,*AssertQuantity2 ,*AssertError, **Param; 00058 00059 /* this flag says where we print linear or log quantity */ 00060 static int *lgQuantityLog; 00061 static long nAsserts=0; 00062 static realnum *wavelength; 00063 00064 /*======================================================================*/ 00065 /*InitAssertResults, this must be called first, done at startup of ParseCommands*/ 00066 void InitAssertResults(void) 00067 { 00068 /* set flag that init was called, and set number of asserts to zero. 00069 * this is done by ParseComments for every model, even when no asserts will 00070 * be done, so do not allocate space at this time */ 00071 lgInitDone = true; 00072 nAsserts = 0; 00073 } 00074 00075 /*======================================================================*/ 00076 /*ParseAssertResults parse the assert command */ 00077 void ParseAssertResults(void) 00078 { 00079 long i , 00080 nelem, 00081 n2; 00082 bool lgEOL; 00083 char chLabel[INPUT_LINE_LENGTH]; 00084 00085 DEBUG_ENTRY( "ParseAssertResults()" ); 00086 00087 if( !lgInitDone ) 00088 { 00089 fprintf( ioQQQ, " ParseAssertResults called before InitAsserResults\n" ); 00090 cdEXIT(EXIT_FAILURE); 00091 } 00092 00093 /* has space been allocated yet? */ 00094 if( !lgSpaceAllocated ) 00095 { 00096 /* - no, we must allocate it */ 00097 /* remember that space has been allocated */ 00098 lgSpaceAllocated = true; 00099 00100 /* create space for the array of labels*/ 00101 chAssertLineLabel = ((char **)MALLOC(NASSERTS*sizeof(char *))); 00102 00103 /* the 2-character string saying what type of assert */ 00104 chAssertType = ((char **)MALLOC(NASSERTS*sizeof(char *))); 00105 00106 /* these are a pair of optional parameters */ 00107 Param = ((double **)MALLOC(NASSERTS*sizeof(double * ))); 00108 00109 /* now fill out the 2D arrays */ 00110 for( i=0; i<NASSERTS; ++i ) 00111 { 00112 chAssertLineLabel[i] = ((char *)MALLOC(NCHAR*sizeof(char ))); 00113 strcpy(chAssertLineLabel[i],"unkn" ); 00114 00115 /* two char plus eol */ 00116 chAssertType[i] = ((char *)MALLOC(3*sizeof(char ))); 00117 strcpy(chAssertType[i],"un" ); 00118 00119 /* these are a pair of optional parameters */ 00120 Param[i] = ((double *)MALLOC(5*sizeof(double ))); 00121 } 00122 00123 /* now make space for the asserted quantities */ 00124 AssertQuantity = ((double *)MALLOC(NASSERTS*sizeof(double ))); 00125 00126 AssertQuantity2 = ((double *)MALLOC(NASSERTS*sizeof(double ))); 00127 00128 /* now the errors */ 00129 AssertError = ((double *)MALLOC(NASSERTS*sizeof(double ))); 00130 00131 /* now the line wavelengths */ 00132 wavelength = ((realnum *)MALLOC(NASSERTS*sizeof(realnum ))); 00133 00134 /* now the flag saying whether should be log */ 00135 lgQuantityLog = ((int *)MALLOC(NASSERTS*sizeof(int ))); 00136 00137 /* the flag for upper, lower limit, or equal */ 00138 chAssertLimit = ((char *)MALLOC(NASSERTS*sizeof(char ))); 00139 } 00140 /* end space allocation - we are ready to roll */ 00141 00142 /* first say we do not know what job to do */ 00143 strcpy( chAssertType[nAsserts] , "un" ); 00144 00145 /* false means print linear quantity - will be set true if entered 00146 * quantity comes in as a log */ 00147 lgQuantityLog[nAsserts] = false; 00148 00149 /* all asserts have option for quantity to be a limit, or the quantity itself */ 00150 if( nMatch("<",input.chCARDCAPS ) ) 00151 { 00152 chAssertLimit[nAsserts] = '<'; 00153 } 00154 else if( nMatch(">",input.chCARDCAPS ) ) 00155 { 00156 chAssertLimit[nAsserts] = '>'; 00157 } 00158 else 00159 { 00160 chAssertLimit[nAsserts] = '='; 00161 } 00162 00163 /* which quantity will we check?, first is */ 00164 00165 /* assert mean ionization */ 00166 if( nMatch("IONI",input.chCARDCAPS ) ) 00167 { 00168 00169 /* say that this will be mean ionization fraction */ 00170 00171 /* f will indicate average over radius, F over volume - 00172 * check whether keyword radius or volume occurs, 00173 * default will be radius */ 00174 if( nMatch("VOLU",input.chCARDCAPS ) ) 00175 { 00176 strcpy( chAssertType[nAsserts] , "F " ); 00177 } 00178 else 00179 { 00180 /* this is default case, Fraction over radius */ 00181 strcpy( chAssertType[nAsserts] , "f " ); 00182 } 00183 00184 /* first get element label and make null terminated string*/ 00185 if( (nelem = GetElem(input.chCARDCAPS)) < 0 ) 00186 { 00187 fprintf( ioQQQ, 00188 " I could not identify an element name on this line.\n"); 00189 fprintf( ioQQQ, " Sorry.\n" ); 00190 cdEXIT(EXIT_FAILURE); 00191 } 00192 ASSERT( nelem>= 0); 00193 ASSERT( nAsserts>= 0); 00194 /* we now have element name, copy 4-char string (elementnames.chElementNameShort[nelem]) 00195 * into array that will be used to get ionization after calculation */ 00196 strcpy( chAssertLineLabel[nAsserts], elementnames.chElementNameShort[nelem] ); 00197 00198 /* now get ionization stage, which will be saved into wavelength */ 00199 i = 5; 00200 wavelength[nAsserts] = 00201 (realnum)FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00202 if( lgEOL ) 00203 { 00204 NoNumb(input.chOrgCard); 00205 } 00206 /* ionization stage must be 1 or greater, but not greater than nelem (c scale)+2 */ 00207 if( wavelength[nAsserts] < 1 || wavelength[nAsserts] > nelem+2 ) 00208 { 00209 fprintf( ioQQQ, 00210 " The ionization stage is inappropriate for this element.\n"); 00211 fprintf( ioQQQ, " Sorry.\n" ); 00212 cdEXIT(EXIT_FAILURE); 00213 } 00214 00215 /* now get ionization fraction, log if number is negative or == 0, 00216 * linear if positive but less than or equal to 1.*/ 00217 AssertQuantity[nAsserts] = 00218 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00219 if( lgEOL ) 00220 NoNumb(input.chOrgCard); 00221 00222 /* optional error, default available (cannot do before loop since we 00223 * do not know how many numbers are on line */ 00224 AssertError[nAsserts] = 00225 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00226 if( lgEOL ) 00227 /* default error was set in define above */ 00228 AssertError[nAsserts] = DEF_ERROR; 00229 00230 if( nMatch( "GRID" , input.chCARDCAPS ) ) 00231 { 00232 long int j; 00233 /* skip nOptimiz values */ 00234 for( j=0; j<optimize.nOptimiz; ++j ) 00235 { 00236 AssertQuantity[nAsserts] = 00237 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00238 } 00239 /*fprintf(ioQQQ,"DEBUG grid HIT %li val %.2f\n", 00240 optimize.nOptimiz , AssertQuantity[nAsserts] );*/ 00241 if( lgEOL ) 00242 { 00243 fprintf(ioQQQ,"PROBLEM this assert grid command does not have enough values. sorry\n"); 00244 } 00245 } 00246 00247 /* now make sure we end up with positive linear ionization fraction that 00248 * is greater than 0 but less than or equal to 1. */ 00249 if( AssertQuantity[nAsserts] <= 0. ) 00250 { 00251 /* log since negative or 0 */ 00252 AssertQuantity[nAsserts] = 00253 pow(10.,AssertQuantity[nAsserts] ); 00254 /* entered as a log, so print as a log too */ 00255 lgQuantityLog[nAsserts] = true; 00256 } 00257 else if( AssertQuantity[nAsserts] > 1. ) 00258 { 00259 fprintf( ioQQQ, 00260 " The ionization fraction must be less than one.\n"); 00261 fprintf( ioQQQ, " Sorry.\n" ); 00262 cdEXIT(EXIT_FAILURE); 00263 } 00264 00265 /* result cannot be zero */ 00266 if( fabs(AssertQuantity[nAsserts]) <= SMALLDOUBLE ) 00267 { 00268 fprintf( ioQQQ, 00269 " The ionization ionization fraction is too small, or zero. Check input\n"); 00270 fprintf( ioQQQ, " Sorry.\n" ); 00271 cdEXIT(EXIT_FAILURE); 00272 } 00273 } 00274 00275 /* molecular fraction averaged over model */ 00276 else if( nMatch("MOLE",input.chCARDCAPS )&&nMatch("FRAC",input.chCARDCAPS ) ) 00277 { 00278 00279 /* say that this will be mean molecular fraction */ 00280 00281 /* mf will indicate average over radius, MF over vol - 00282 * check whether keyword radius or volume occurs, 00283 * default will be radius */ 00285 if( nMatch("VOLU",input.chCARDCAPS ) ) 00286 { 00287 strcpy( chAssertType[nAsserts] , "MF" ); 00288 } 00289 else 00290 { 00291 /* this is default case, Fraction over radius */ 00292 strcpy( chAssertType[nAsserts] , "mf" ); 00293 } 00294 00295 i = nMatch(" H2 ",input.chCARDCAPS ); 00296 if( i ) 00297 { 00298 strcpy( chAssertLineLabel[nAsserts], "H2 " ); 00299 /* increment to get past the label */ 00300 i += 3; 00301 } 00302 else if( nMatch(" CO ",input.chCARDCAPS ) ) 00303 { 00304 strcpy( chAssertLineLabel[nAsserts], "CO " ); 00305 i = 5; 00306 } 00307 else 00308 { 00309 fprintf( ioQQQ, 00310 " I could not identify CO or H2 on this line.\n"); 00311 fprintf( ioQQQ, " Sorry.\n" ); 00312 cdEXIT(EXIT_FAILURE); 00313 } 00314 00315 /* not meaningful */ 00316 wavelength[nAsserts] = 0; 00317 00318 /* i was set above for start of scan */ 00319 /* now get log of molecular fraction */ 00320 AssertQuantity[nAsserts] = 00321 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00322 if( lgEOL ) 00323 { 00324 NoNumb(input.chOrgCard); 00325 } 00326 if( AssertQuantity[nAsserts] <= 0. ) 00327 { 00328 /* if negative then entered as log, but we will compare with linear */ 00329 AssertQuantity[nAsserts] = 00330 pow(10.,AssertQuantity[nAsserts] ); 00331 } 00332 00333 /* optional error, default available (cannot do before loop since we 00334 * do not know how many numbers are on line */ 00335 AssertError[nAsserts] = 00336 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00337 if( lgEOL ) 00338 { 00339 /* default error was set in define above */ 00340 AssertError[nAsserts] = DEF_ERROR; 00341 } 00342 /* print results as logs */ 00343 lgQuantityLog[nAsserts] = true; 00344 } 00345 00346 /* assert line "LINE" -- key is ine_ since linear option appears on some commands */ 00347 else if( nMatch(" LINE ",input.chCARDCAPS ) ) 00348 { 00349 if( nMatch("LUMI",input.chCARDCAPS ) || nMatch("INTE",input.chCARDCAPS )) 00350 { 00351 /* say that this is a line luminosity or intensity*/ 00352 strcpy( chAssertType[nAsserts] , "Ll" ); 00353 00354 /* entered as a log, so print as a log too */ 00355 lgQuantityLog[nAsserts] = true; 00356 } 00357 else 00358 { 00359 /* say that this is line relative to norm line - this is the default */ 00360 strcpy( chAssertType[nAsserts] , "Lr" ); 00361 00362 /* entered linear quantity, so print as linear too */ 00363 lgQuantityLog[nAsserts] = false; 00364 } 00365 00366 /* this will check a line intensity, first get labels within quotes, 00367 * will be null terminated */ 00368 GetQuote( chLabel , input.chCARDCAPS , true ); 00369 00370 /* check that there were exactly 4 characters in the label*/ 00371 if( chLabel[4] != '\0' ) 00372 { 00373 fprintf( ioQQQ, 00374 " The label must be exactly 4 char long, between double quotes.\n"); 00375 fprintf( ioQQQ, " Sorry.\n" ); 00376 cdEXIT(EXIT_FAILURE); 00377 } 00378 00379 /* copy string into array */ 00380 strcpy( chAssertLineLabel[nAsserts], chLabel ); 00381 00382 /* now get line wavelength */ 00383 i = 5; 00384 wavelength[nAsserts] = 00385 (realnum)FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00386 if( lgEOL ) 00387 { 00388 NoNumb(input.chOrgCard); 00389 } 00390 /* check for optional micron or cm units, else interpret as Angstroms */ 00391 if( input.chCARDCAPS[i-1] == 'M' ) 00392 { 00393 /* microns */ 00394 wavelength[nAsserts] *= 1e4; 00395 } 00396 else if( input.chCARDCAPS[i-1] == 'C' ) 00397 { 00398 /* centimeters */ 00399 wavelength[nAsserts] *= 1e8; 00400 } 00401 00402 /* now get intensity or luminosity - 00403 * rel intensity is linear and intensity or luminosity are log */ 00404 AssertQuantity[nAsserts] = 00405 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00406 if( lgEOL ) 00407 { 00408 NoNumb(input.chOrgCard); 00409 } 00410 /* luminosity was entered as a log */ 00411 if( lgQuantityLog[nAsserts] ) 00412 { 00413 if( AssertQuantity[nAsserts] > DBL_MAX_10_EXP || 00414 AssertQuantity[nAsserts] < -DBL_MAX_10_EXP ) 00415 { 00416 fprintf( ioQQQ, 00417 " The asserted quantity is a log, but is too large or small, value is %e.\n",AssertQuantity[nAsserts] ); 00418 fprintf( ioQQQ, " I would crash if I tried to use it.\n Sorry.\n" ); 00419 cdEXIT(EXIT_FAILURE); 00420 } 00421 AssertQuantity[nAsserts] = 00422 pow(10.,AssertQuantity[nAsserts] ); 00423 } 00424 if( AssertQuantity[nAsserts]<= 0. ) 00425 { 00426 fprintf( ioQQQ, 00427 " The relative intensity must be positive, and was %e.\n",AssertQuantity[nAsserts] ); 00428 fprintf( ioQQQ, " Sorry.\n" ); 00429 cdEXIT(EXIT_FAILURE); 00430 } 00431 00432 /* optional error, default available */ 00433 AssertError[nAsserts] = 00434 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00435 if( lgEOL ) 00436 { 00437 /* default error was set in define above */ 00438 AssertError[nAsserts] = DEF_ERROR; 00439 } 00440 } 00441 00442 /* assert departure coefficients */ 00443 else if( nMatch("CASE",input.chCARDCAPS ) ) 00444 { 00445 /* this is Case B for some element */ 00446 strcpy( chAssertType[nAsserts] , "CB" ); 00447 i = 5; 00448 /* this is relative error */ 00449 AssertError[nAsserts] = 00450 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00451 if( lgEOL ) 00452 /* default error was set in define above */ 00453 AssertError[nAsserts] = DEF_ERROR; 00454 AssertQuantity[nAsserts] = 0; 00455 wavelength[nAsserts] = 0.; 00456 00457 /* faint option - do not test line if relative intensity is less 00458 * than entered value */ 00459 if( (i=nMatch("FAINT",input.chCARDCAPS ))>0 ) 00460 { 00461 Param[nAsserts][4] = 00462 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00463 if( lgEOL ) 00464 { 00465 /* did not get 2 numbers */ 00466 fprintf(ioQQQ," The assert Case B faint option must have a number," 00467 " the relative intensity of the fainest line to assert.\n"); 00468 cdEXIT(EXIT_FAILURE); 00469 } 00470 /* number is log if <= 0 */ 00471 if( Param[nAsserts][4]<=0. ) 00472 Param[nAsserts][4] = pow(10., Param[nAsserts][4] ); 00473 } 00474 else 00475 { 00476 /* use default - include everything*/ 00477 Param[nAsserts][4] = SMALLFLOAT; 00478 } 00479 00480 /* range option - to limit check on a certain wavelength range */ 00481 if( (i=nMatch("RANG",input.chCARDCAPS ))>0 ) 00482 { 00483 Param[nAsserts][2] = 00484 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00485 Param[nAsserts][3] = 00486 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00487 if( lgEOL ) 00488 { 00489 /* did not get 2 numbers */ 00490 fprintf(ioQQQ," The assert Case B range option must have two numbers," 00491 " the lower and upper limit to the wavelengths in Angstroms.\n"); 00492 fprintf(ioQQQ," There must be a total of three numbers on the line," 00493 " the relative error followed by the lower and upper limits to the " 00494 "wavelength in Angstroms.\n"); 00495 cdEXIT(EXIT_FAILURE); 00496 } 00497 if( Param[nAsserts][2]>Param[nAsserts][3]) 00498 { 00499 /* make sure in increasing order */ 00500 double sav = Param[nAsserts][3]; 00501 Param[nAsserts][3] = Param[nAsserts][2]; 00502 Param[nAsserts][2] = sav; 00503 } 00504 } 00505 else 00506 { 00507 /* use default - include everything*/ 00508 Param[nAsserts][2] = 0.; 00509 Param[nAsserts][3] = 1e30; 00510 } 00511 /* assert case b for H - O checking against Hummer & Storey tables */ 00512 if( nMatch("H-LI",input.chCARDCAPS ) ) 00513 { 00514 /* H-like - now get an element */ 00515 if( (nelem = GetElem( input.chCARDCAPS )) < 0 ) 00516 { 00517 /* no name found */ 00518 fprintf(ioQQQ, "assert H-like case B did not find an element on this line, sorry %s\n", 00519 input.chCARDCAPS ); 00520 cdEXIT(EXIT_FAILURE); 00521 } 00522 if( nelem>7 ) 00523 { 00524 /* beyond reach of tables */ 00525 fprintf(ioQQQ, "assert H-like cannot do elements heavier than O, sorry %s\n", 00526 input.chCARDCAPS ); 00527 cdEXIT(EXIT_FAILURE); 00528 } 00529 Param[nAsserts][0] = ipH_LIKE; 00530 Param[nAsserts][1] = nelem; 00531 /* generate string to find simple prediction, as in "Ca B" */ 00532 strcpy( chAssertLineLabel[nAsserts], "Ca "); 00533 if( nMatch("CASE A ",input.chCARDCAPS ) ) 00534 strcat( chAssertLineLabel[nAsserts] , "A"); 00535 else 00536 strcat( chAssertLineLabel[nAsserts] , "B"); 00537 } 00538 else if( nMatch("HE-L",input.chCARDCAPS ) ) 00539 { 00540 /* He-like - only helium itself */ 00541 Param[nAsserts][0] = ipHE_LIKE; 00542 Param[nAsserts][1] = ipHELIUM; 00543 00544 strcpy( chAssertLineLabel[nAsserts] , "+Col"); 00545 } 00546 } 00547 else if( nMatch("DEPA",input.chCARDCAPS ) ) 00548 { 00549 i = 5; 00550 /* get expected average departure coefficient, almost certainly 1 */ 00551 AssertQuantity[nAsserts] = 00552 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00553 if( lgEOL ) 00554 NoNumb(input.chOrgCard); 00555 00556 /* this is relative error, max departure from unity of any level or std */ 00557 AssertError[nAsserts] = 00558 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00559 if( lgEOL ) 00560 /* default error was set in define above */ 00561 AssertError[nAsserts] = DEF_ERROR; 00562 00563 /* H-like key means do one of the hydrogenic ions */ 00564 if( nMatch("H-LI",input.chCARDCAPS ) ) 00565 { 00566 /* label is dHII */ 00567 strcpy( chAssertLineLabel[nAsserts] , "dHyd" ); 00568 /* remember this is departure coefficient for some element */ 00569 strcpy( chAssertType[nAsserts] , "D " ); 00570 /* now get element number for h ion from element name on card */ 00571 if( (wavelength[nAsserts] = (realnum)GetElem(input.chCARDCAPS)) < 0 ) 00572 { 00573 fprintf( ioQQQ, 00574 " I could not identify an element name on this line.\n"); 00575 fprintf( ioQQQ, " Sorry.\n" ); 00576 cdEXIT(EXIT_FAILURE); 00577 } 00578 } 00579 00580 /* He-like key means do one of the helike ions */ 00581 else if( nMatch("HE-L",input.chCARDCAPS ) ) 00582 { 00583 /* label is dHII */ 00584 strcpy( chAssertLineLabel[nAsserts] , "d He" ); 00585 /* remember this is departure coefficient for some element */ 00586 strcpy( chAssertType[nAsserts] , "De" ); 00587 /* now get element number for h ion from element name on card */ 00588 if( (wavelength[nAsserts] = (realnum)GetElem(input.chCARDCAPS)) < 0 ) 00589 { 00590 fprintf( ioQQQ, 00591 " I could not identify an element name on this line.\n"); 00592 fprintf( ioQQQ, " Sorry.\n" ); 00593 cdEXIT(EXIT_FAILURE); 00594 } 00595 } 00596 00597 /* this is the large FeII ion */ 00598 else if( nMatch("FEII",input.chCARDCAPS ) || nMatch("FE II",input.chCARDCAPS ) ) 00599 { 00600 /* label */ 00601 strcpy( chAssertLineLabel[nAsserts] , "d Fe" ); 00602 /* remember this is departure coefficient for feii */ 00603 strcpy( chAssertType[nAsserts] , "d " ); 00604 /* the wavelength is meaningless, but put in 2 since FeII */ 00605 wavelength[nAsserts] = 2; 00606 } 00607 00608 /* this is H- h minus */ 00609 else if( nMatch("HMIN",input.chCARDCAPS ) ) 00610 { 00611 /* label */ 00612 strcpy( chAssertLineLabel[nAsserts] , "d H-" ); 00613 /* remember this is departure coefficient for H- */ 00614 strcpy( chAssertType[nAsserts] , "d-" ); 00615 /* the wavelength is meaningless */ 00616 wavelength[nAsserts] = -1; 00617 } 00618 else 00619 { 00620 fprintf( ioQQQ, 00621 " There must be a second key: FEII, H-LIke, HMINus, or HE-Like followed by element.\n"); 00622 fprintf( ioQQQ, " Sorry.\n" ); 00623 cdEXIT(EXIT_FAILURE); 00624 } 00625 00626 /* last check for key "excited" - which means to skip the ground state */ 00627 if( nMatch("EXCI",input.chCARDCAPS ) ) 00628 { 00629 /* this is lowest level - do not do 0 */ 00630 AssertQuantity2[nAsserts] = 1.; 00631 } 00632 else 00633 { 00634 /* do the ground state */ 00635 AssertQuantity2[nAsserts] = 0.; 00636 } 00637 } 00638 00639 /* assert some results from map */ 00640 else if( nMatch(" MAP",input.chCARDCAPS ) ) 00641 { 00642 00643 /* must have heating or cooling, since will check one or the other */ 00644 /* check heating cooling results from map at some temperature */ 00645 if( nMatch("HEAT",input.chCARDCAPS ) ) 00646 { 00647 strcpy( chAssertType[nAsserts] , "mh" ); 00648 } 00649 else if( nMatch("COOL",input.chCARDCAPS ) ) 00650 { 00651 strcpy( chAssertType[nAsserts] , "mc" ); 00652 } 00653 else 00654 { 00655 fprintf( ioQQQ, 00656 " There must be a second key, HEATing or COOLing.\n"); 00657 fprintf( ioQQQ, " Sorry.\n" ); 00658 cdEXIT(EXIT_FAILURE); 00659 } 00660 00661 /* i was set above for start of scan */ 00662 /* now get temperature for AssertQuantity2 array*/ 00663 i = 5; 00664 AssertQuantity2[nAsserts] = 00665 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00666 if( lgEOL ) 00667 { 00668 NoNumb(input.chOrgCard); 00669 } 00670 00671 if( AssertQuantity2[nAsserts] <= 10. ) 00672 { 00673 /* entered as log, but we will compare with linear */ 00674 AssertQuantity2[nAsserts] = 00675 pow(10.,AssertQuantity2[nAsserts] ); 00676 } 00677 00678 /* print the temperature in the wavelength column */ 00679 wavelength[nAsserts] = (realnum)AssertQuantity2[nAsserts]; 00680 00681 /* heating or cooling, both log, put into error */ 00682 AssertQuantity[nAsserts] = 00683 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00684 if( lgEOL ) 00685 { 00686 NoNumb(input.chOrgCard); 00687 } 00688 00689 /* AssertQuantity array will have heating or cooling */ 00690 AssertQuantity[nAsserts] = pow(10., AssertQuantity[nAsserts]); 00691 00692 /* optional error, default available (cannot do before loop since we 00693 * do not know how many numbers are on line */ 00694 AssertError[nAsserts] = 00695 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00696 if( lgEOL ) 00697 { 00698 /* default error was set in define above */ 00699 AssertError[nAsserts] = DEF_ERROR; 00700 } 00701 00702 /* entered as a log, so print as a log too */ 00703 lgQuantityLog[nAsserts] = true; 00704 } 00705 00706 /* assert column density of something */ 00707 else if( nMatch("COLU",input.chCARDCAPS ) ) 00708 { 00709 /* this is column density */ 00710 strcpy( chAssertType[nAsserts] , "cd" ); 00711 00712 /* this says to look for molecular column density, also could be ion stage */ 00713 wavelength[nAsserts] = 0; 00714 00715 /* we want to remember where the match occurred within the string 00716 * since do not want to count the 2 as the first number */ 00717 /* NB - can't put this expression in the if since many compilers warn against this */ 00718 if( (i = nMatch(" H2 ",input.chCARDCAPS )) != 0 ) 00719 { 00720 strcpy( chAssertLineLabel[nAsserts], "H2 " ); 00721 /* increment to get past the 2 in the label */ 00722 i += 3; 00723 if( nMatch( "LEVE" , input.chCARDCAPS ) ) 00724 { 00725 /* this is option for level-specific column density, 00726 * next two numbers must be v then J */ 00727 Param[nAsserts][0] = 00728 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00729 Param[nAsserts][1] = 00730 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00731 if( lgEOL ) 00732 NoNumb( input.chCARDCAPS ); 00733 /* wavelength will be 10. * vib + rot */ 00734 wavelength[nAsserts] = (realnum)(100.*Param[nAsserts][0] + Param[nAsserts][1]); 00735 } 00736 else 00737 { 00738 /* these are flags saying not to do state specific column densities */ 00739 Param[nAsserts][0] = -1.; 00740 Param[nAsserts][1] = -1.; 00741 } 00742 } 00743 else if( (i = nMatch("H3+ ",input.chCARDCAPS )) != 0 ) 00744 { 00745 strcpy( chAssertLineLabel[nAsserts], "H3+ " ); 00746 /* increment to get past the 2 in the label */ 00747 i += 3; 00748 } 00749 else if( (i = nMatch("H2+ ",input.chCARDCAPS )) != 0 ) 00750 { 00751 strcpy( chAssertLineLabel[nAsserts], "H2+ " ); 00752 /* increment to get past the 2 in the label */ 00753 i += 3; 00754 } 00755 else if( (i = nMatch(" H- ",input.chCARDCAPS )) != 0 ) 00756 { 00757 strcpy( chAssertLineLabel[nAsserts], "H- " ); 00758 /* increment to get past the 2 in the label */ 00759 i += 3; 00760 } 00761 else if( (i = nMatch("H2G ",input.chCARDCAPS )) != 0 ) 00762 { 00763 strcpy( chAssertLineLabel[nAsserts], "H2g " ); 00764 /* increment to get past the 2 in the label */ 00765 i += 3; 00766 } 00767 else if( (i = nMatch("H2* ",input.chCARDCAPS )) != 0 ) 00768 { 00769 strcpy( chAssertLineLabel[nAsserts], "H2* " ); 00770 /* increment to get past the 2 in the label */ 00771 i += 3; 00772 } 00773 else if( (i = nMatch("HEH+",input.chCARDCAPS )) != 0 ) 00774 { 00775 strcpy( chAssertLineLabel[nAsserts], "HeH+" ); 00776 /* increment to get past the 2 in the label */ 00777 i += 3; 00778 } 00779 else if( (i = nMatch(" O2 ",input.chCARDCAPS )) != 0 ) 00780 { 00781 strcpy( chAssertLineLabel[nAsserts], "O2 " ); 00782 /* increment to get past the 2 in the label */ 00783 i += 3; 00784 } 00785 else if( (i = nMatch("H2O ",input.chCARDCAPS )) != 0 ) 00786 { 00787 strcpy( chAssertLineLabel[nAsserts], "H2O " ); 00788 /* increment to get past the 2 in the label */ 00789 i += 3; 00790 } 00791 else if( (i = nMatch(" C2 ",input.chCARDCAPS ) ) !=0 ) 00792 { 00793 strcpy( chAssertLineLabel[nAsserts], "C2 " ); 00794 /* increment to get past the 2 in the label */ 00795 i += 3; 00796 } 00797 else if( (i = nMatch(" C3 ",input.chCARDCAPS ) ) !=0 ) 00798 { 00799 strcpy( chAssertLineLabel[nAsserts], "C3 " ); 00800 /* increment to get past the 3 in the label */ 00801 i += 3; 00802 } 00803 else if( nMatch(" CO ",input.chCARDCAPS ) ) 00804 { 00805 strcpy( chAssertLineLabel[nAsserts], "CO " ); 00806 i = 5; 00807 } 00808 else if( nMatch("SIO ",input.chCARDCAPS ) ) 00809 { 00810 strcpy( chAssertLineLabel[nAsserts], "SiO " ); 00811 i = 5; 00812 } 00813 else if( nMatch(" OH ",input.chCARDCAPS ) ) 00814 { 00815 strcpy( chAssertLineLabel[nAsserts], "OH " ); 00816 i = 5; 00817 } 00818 else if( nMatch(" CN ",input.chCARDCAPS ) ) 00819 { 00820 strcpy( chAssertLineLabel[nAsserts], "CN " ); 00821 i = 5; 00822 } 00823 else if( nMatch(" CH ",input.chCARDCAPS ) ) 00824 { 00825 strcpy( chAssertLineLabel[nAsserts], "CH " ); 00826 i = 5; 00827 } 00828 else if( nMatch(" CH+",input.chCARDCAPS ) ) 00829 { 00830 strcpy( chAssertLineLabel[nAsserts], "CH+ " ); 00831 i = 5; 00832 } 00833 else 00834 { 00835 fprintf( ioQQQ, 00836 " I could not identify H2, H3+, H2+, H2g, H2*, H2H+, CO, O2, SiO, OH, C2 or C3 or on this line.\n"); 00837 fprintf( ioQQQ, " Sorry.\n" ); 00838 cdEXIT(EXIT_FAILURE); 00839 } 00840 00841 /* i was set above for start of scan */ 00842 /* now get log of column density */ 00843 AssertQuantity[nAsserts] = 00844 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00845 if( lgEOL ) 00846 { 00847 NoNumb(input.chOrgCard); 00848 } 00849 /* entered as log, but we will compare with linear */ 00850 AssertQuantity[nAsserts] = 00851 pow(10.,AssertQuantity[nAsserts] ); 00852 00853 /* optional error, default available (cannot do before loop since we 00854 * do not know how many numbers are on line */ 00855 AssertError[nAsserts] = 00856 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00857 if( lgEOL ) 00858 { 00859 /* default error was set in define above */ 00860 AssertError[nAsserts] = DEF_ERROR; 00861 } 00862 /* the keyword log is special for this case, since H2 and CO column densities can 00863 * be so very unstable. look for work log, in which case the error is log not linear. 00864 * main column is always a log */ 00865 if( nMatch( " LOG" , input.chCARDCAPS ) ) 00866 { 00867 AssertError[nAsserts] = pow(10. , AssertError[nAsserts] ); 00868 } 00869 00870 /* entered as a log, so print as a log too although asserted quantity is now linear */ 00871 lgQuantityLog[nAsserts] = true; 00872 } 00873 /* assert rate H2 forms on grain surfaces */ 00874 else if( nMatch("GRAI",input.chCARDCAPS ) && nMatch(" H2 ",input.chCARDCAPS ) ) 00875 { 00876 /* this flag will mean h2 form on grains */ 00877 strcpy( chAssertType[nAsserts] , "g2" ); 00878 /* a label */ 00879 strcpy( chAssertLineLabel[nAsserts], "R H2" ); 00880 /* now get the first number on the line, which must be the 2 in H2 */ 00881 i = 5; 00882 nelem = 00883 (long int)FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00884 if( nelem!=2 ) 00885 { 00886 fprintf( ioQQQ, 00887 " I did not find a 2, as in H2, as the first number on this line.\n"); 00888 fprintf( ioQQQ, 00889 " The rate should be the second number.\n"); 00890 fprintf( ioQQQ, " Sorry.\n" ); 00891 cdEXIT(EXIT_FAILURE); 00892 } 00893 /* now past the 2 in h2, get the real number */ 00894 AssertQuantity[nAsserts] = 00895 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00896 if( lgEOL ) 00897 { 00898 NoNumb(input.chOrgCard); 00899 } 00900 /* if negative (almost certainly the case) then the log of the rate coefficient */ 00901 if( AssertQuantity[nAsserts] <=0. ) 00902 AssertQuantity[nAsserts] = pow(10.,AssertQuantity[nAsserts] ); 00903 /* will not use this */ 00904 wavelength[nAsserts] = 0; 00905 00906 /* optional error, default available (cannot do before loop since we 00907 * do not know how many numbers are on line */ 00908 AssertError[nAsserts] = 00909 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00910 if( lgEOL ) 00911 { 00912 /* default error was set in define above */ 00913 AssertError[nAsserts] = DEF_ERROR; 00914 /* want to print as a log since so small */ 00915 lgQuantityLog[nAsserts] = true; 00916 } 00917 } 00918 00919 /* assert grain potential */ 00920 else if( nMatch( "GRAI", input.chCARDCAPS ) && nMatch( "POTE", input.chCARDCAPS ) ) 00921 { 00922 /* this flag will mean grain potential */ 00923 strcpy( chAssertType[nAsserts] , "gp" ); 00924 /* a label */ 00925 strcpy( chAssertLineLabel[nAsserts], "GPot" ); 00926 /* now get the first number on the line */ 00927 i = 5; 00928 /* grain bin number */ 00929 wavelength[nAsserts] = (realnum)FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00930 /* the potential itself, in volt, always linear */ 00931 AssertQuantity[nAsserts] = FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00932 00933 if( lgEOL ) 00934 { 00935 NoNumb(input.chOrgCard); 00936 } 00937 00938 /* optional error, default available (cannot do before loop since we 00939 * do not know how many numbers are on line */ 00940 AssertError[nAsserts] = FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 00941 00942 if( lgEOL ) 00943 { 00944 /* default error was set in define above */ 00945 AssertError[nAsserts] = DEF_ERROR; 00946 } 00947 } 00948 00949 /* assert mean temperature, assert temperature hydrogen 2 8000 */ 00950 else if( nMatch("TEMP",input.chCARDCAPS ) ) 00951 { 00952 /* say that this will be mean temperature, electron or grain */ 00953 00954 /* t will indicate temperature average over radius, T over volume - 00955 * check whether keyword radius or volume occurs, 00956 * default will be radius */ 00957 if( nMatch("VOLU",input.chCARDCAPS ) ) 00958 { 00959 strcpy( chAssertType[nAsserts] , "T "); 00960 } 00961 else 00962 { 00963 /* this is default case, Fraction over radius */ 00964 strcpy( chAssertType[nAsserts] , "t "); 00965 } 00966 00967 /* first look for keyword Grains, since label silicate may be on it, 00968 * and this would trigger the element search */ 00969 if( nMatch("GRAI",input.chCARDCAPS ) ) 00970 { 00971 /* grains, copy 4-char string "grai" 00972 * into array that will be used to get temperature after calculation */ 00973 strcpy( chAssertLineLabel[nAsserts], "GTem" ); 00974 /* this is to make sure that pointer to grain type is valid, we check 00975 * that it is less than this below */ 00976 nelem = NDUST-2; 00977 /* the first numerical arg on the temper grain command is the grain 00978 * type as defined in Hazy, counting from 1. When it is used to 00979 * to get mean temps below we will subtract 1 from this to get onto 00980 * the c scale. but must leave on physical scale here to pass sanity 00981 * checks that occur below */ 00982 } 00983 00984 /* face is temperature at illuminated face of cloud */ 00985 else if( nMatch("FACE",input.chCARDCAPS ) ) 00986 { 00987 strcpy( chAssertLineLabel[nAsserts], "face" ); 00988 /* not used so zero out */ 00989 nelem = 0; 00990 } 00991 00992 /* mean H2 temperature */ 00993 else if( nMatch(" H2 ",input.chCARDCAPS ) ) 00994 { 00995 strcpy( chAssertLineLabel[nAsserts], "H2 " ); 00996 /* not used so zero out */ 00997 nelem = 0; 00998 } 00999 01000 /* look for element label within quotes, 01001 * will be null terminated */ 01002 else if( (nelem = GetElem(input.chCARDCAPS)) >= 0 ) 01003 { 01004 /* we now have element name, copy 4-char string (elementnames.chElementNameShort[nelem]) 01005 * into array that will be used to get ionization after calculation */ 01006 strcpy( chAssertLineLabel[nAsserts], elementnames.chElementNameShort[nelem] ); 01007 } 01008 else 01009 { 01010 fprintf( ioQQQ, 01011 " I could not identify an element name on this line.\n"); 01012 fprintf( ioQQQ, " Sorry.\n" ); 01013 cdEXIT(EXIT_FAILURE); 01014 } 01015 01016 /* now get ionization stage (or grain type), which will be saved into wavelength */ 01017 i = 5; 01018 if( strcmp( chAssertLineLabel[nAsserts], "face" )== 0) 01019 { 01020 /* this option is temperature at illuminated face so no element */ 01021 wavelength[nAsserts] = 0; 01022 } 01023 else if( strcmp( chAssertLineLabel[nAsserts], "H2 " )== 0) 01024 { 01025 int j; 01026 /* this option is temperature of H2 so element is zero */ 01027 wavelength[nAsserts] = 0; 01028 /* first find the 2 in H2 */ 01029 i = 5; 01030 j = (int)FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01031 if( j != 2 ) 01032 TotalInsanity(); 01033 ++i; 01034 } 01035 else 01036 { 01037 wavelength[nAsserts] = 01038 (realnum)FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01039 if( lgEOL ) 01040 { 01041 NoNumb(input.chOrgCard); 01042 } 01043 /* ionization stage must be 1 or greater, but not greater than nelem (c scale)+2 */ 01044 if( wavelength[nAsserts] < 1 || wavelength[nAsserts] > nelem+2 ) 01045 { 01046 fprintf( ioQQQ, 01047 " The ionization stage is inappropriate for this element.\n"); 01048 fprintf( ioQQQ, " Sorry.\n" ); 01049 cdEXIT(EXIT_FAILURE); 01050 } 01051 } 01052 01053 /* now get temperature, log if number is <= 10, else linear */ 01054 AssertQuantity[nAsserts] = 01055 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01056 if( lgEOL ) 01057 { 01058 NoNumb(input.chOrgCard); 01059 } 01060 01061 /* optional error, default available (cannot do before loop since we 01062 * do not know how many numbers are on line */ 01063 AssertError[nAsserts] = 01064 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01065 if( lgEOL ) 01066 /* default error was set in define above */ 01067 AssertError[nAsserts] = DEF_ERROR; 01068 01069 if( nMatch( "GRID" , input.chCARDCAPS ) ) 01070 { 01071 long int j; 01072 /* skip nOptimiz values */ 01073 for( j=0; j<optimize.nOptimiz; ++j ) 01074 { 01075 AssertQuantity[nAsserts] = 01076 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01077 } 01078 /*fprintf(ioQQQ,"DEBUG grid HIT %li val %.2f\n", 01079 optimize.nOptimiz , AssertQuantity[nAsserts] );*/ 01080 if( lgEOL ) 01081 { 01082 fprintf(ioQQQ,"PROBLEM this assert grid command does not have enough values. sorry\n"); 01083 } 01084 } 01085 01086 /* now make sure we end up with positive linear temperature 01087 * number is log if <=10 unless linear keyword appears */ 01088 if( AssertQuantity[nAsserts] <= 10. && !nMatch( "LINE" ,input.chCARDCAPS ) ) 01089 { 01090 /* log since negative or 0 */ 01091 AssertQuantity[nAsserts] = 01092 pow(10.,AssertQuantity[nAsserts] ); 01093 /* entered as a log, so print as a log too */ 01094 lgQuantityLog[nAsserts] = true; 01095 } 01096 } 01097 01098 /* assert log of helium hydrogen ionization correction factor */ 01099 else if( nMatch("HHEI",input.chCARDCAPS ) ) 01100 { 01101 /* this flag will mean H-He icf */ 01102 strcpy( chAssertType[nAsserts] , "hh" ); 01103 /* say that this is zone numbering */ 01104 strcpy( chAssertLineLabel[nAsserts], "HHei" ); 01105 01106 /* now get the ionization correction factor, it is always the linear 01107 * quantity itself, since can be positive or negative*/ 01108 i = 5; 01109 AssertQuantity[nAsserts] = 01110 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01111 if( lgEOL ) 01112 { 01113 NoNumb(input.chOrgCard); 01114 } 01115 /* will not use this */ 01116 wavelength[nAsserts] = 0; 01117 01118 /* optional error, default available (cannot do before loop since we 01119 * do not know how many numbers are on line */ 01120 AssertError[nAsserts] = 01121 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01122 if( lgEOL ) 01123 { 01124 /* default error was set in define above */ 01125 AssertError[nAsserts] = DEF_ERROR; 01126 } 01127 } 01128 01129 /* this large H2 molecule */ 01130 else if( nMatch(" H2 ",input.chCARDCAPS ) ) 01131 { 01132 /* ortho to para ratio for last computed zone */ 01133 if( nMatch("ORTH",input.chCARDCAPS ) ) 01134 { 01135 /* this flag will mean ortho to para ratio */ 01136 strcpy( chAssertType[nAsserts] , "or" ); 01137 /* say that this is ortho to para density ratio */ 01138 strcpy( chAssertLineLabel[nAsserts], "orth" ); 01139 } 01140 else 01141 { 01142 fprintf( ioQQQ, 01143 " I could not identify a second keyword on this line.\n"); 01144 fprintf( ioQQQ, " Sorry.\n" ); 01145 cdEXIT(EXIT_FAILURE); 01146 } 01147 01148 /* now get the first number, which better be the 2 in H2 */ 01149 i = 5; 01150 n2 = (long)FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01151 if( lgEOL || n2 != 2 ) 01152 { 01153 NoNumb(input.chOrgCard); 01154 } 01155 AssertQuantity[nAsserts] = 01156 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01157 /* will not use this */ 01158 wavelength[nAsserts] = 0; 01159 01160 /* optional error, default available (cannot do before loop since we 01161 * do not know how many numbers are on line */ 01162 AssertError[nAsserts] = 01163 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01164 if( lgEOL ) 01165 { 01166 /* default error was set in define above */ 01167 AssertError[nAsserts] = DEF_ERROR; 01168 } 01169 } 01170 01171 /* assert (probably upper limit to) number of zones */ 01172 else if( nMatch("NZON",input.chCARDCAPS ) ) 01173 { 01174 /* this flag will mean number of zones */ 01175 strcpy( chAssertType[nAsserts] , "z " ); 01176 /* say that this is zone numbering */ 01177 strcpy( chAssertLineLabel[nAsserts], "zone" ); 01178 01179 /* now get number of zones, which will be saved into AssertQuantity */ 01180 i = 5; 01181 wavelength[nAsserts] = 01182 (realnum)FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01183 if( lgEOL ) 01184 NoNumb(input.chOrgCard); 01185 AssertQuantity[nAsserts] = (double)wavelength[nAsserts]; 01186 01187 /* optional error */ 01188 AssertError[nAsserts] = 01189 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01190 if( lgEOL ) 01191 /* default error was set in define above */ 01192 AssertError[nAsserts] = DEF_ERROR; 01193 } 01194 01195 /* assert (probably upper limit to) standard deviation of pressure across model */ 01196 else if( nMatch("PRES",input.chCARDCAPS ) ) 01197 { 01198 /* this flag indicates ratio of standard deviation to the mean pressure */ 01199 strcpy( chAssertType[nAsserts] , "pr" ); 01200 /* say that this is error in pressure */ 01201 strcpy( chAssertLineLabel[nAsserts], "pres" ); 01202 01203 /* now get the pressure error, which will be saved into wavelength 01204 * in nearly all cases this is limit to error */ 01205 i = 5; 01206 wavelength[nAsserts] = 0; 01207 AssertQuantity[nAsserts] = (double)FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01208 if( lgEOL ) 01209 { 01210 NoNumb(input.chOrgCard); 01211 } 01212 else if( AssertQuantity[nAsserts] <= 0.) 01213 { 01214 /* number <= 0 is log of error */ 01215 AssertQuantity[nAsserts] = pow(10. , AssertQuantity[nAsserts]); 01216 } 01217 01218 /* optional error, default available (cannot do before loop since we 01219 * do not know how many numbers are on line */ 01220 AssertError[nAsserts] = 01221 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01222 if( lgEOL ) 01223 { 01224 /* default error was set in define above */ 01225 AssertError[nAsserts] = DEF_ERROR; 01226 } 01227 } 01228 01229 else if( nMatch("PRAD",input.chCARDCAPS ) && nMatch("DMAX",input.chCARDCAPS ) ) 01230 { 01231 /* assert pradmax - max ratio of rad to gas pressure */ 01232 /* this flag indicates ratio of rad to gas pressure */ 01233 strcpy( chAssertType[nAsserts] , "RM" ); 01234 /* say that this is error in pressure */ 01235 strcpy( chAssertLineLabel[nAsserts], "Prad" ); 01236 01237 /* now get the pressure error, which will be saved into wavelength 01238 * in nearly all cases this is limit to error */ 01239 i = 5; 01240 wavelength[nAsserts] = 0; 01241 AssertQuantity[nAsserts] = (double)FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01242 if( lgEOL ) 01243 { 01244 NoNumb(input.chOrgCard); 01245 } 01246 else if( AssertQuantity[nAsserts] <= 0.) 01247 { 01248 /* number <= 0 is log of error */ 01249 AssertQuantity[nAsserts] = pow(10. , AssertQuantity[nAsserts]); 01250 } 01251 01252 /* optional error, default available (cannot do before loop since we 01253 * do not know how many numbers are on line */ 01254 AssertError[nAsserts] = 01255 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01256 if( lgEOL ) 01257 { 01258 /* default error was set in define above */ 01259 AssertError[nAsserts] = DEF_ERROR; 01260 } 01261 } 01262 01263 /* assert secondary ionization rate, csupra */ 01264 else if( nMatch("CSUP",input.chCARDCAPS ) ) 01265 { 01266 /* this flag will mean secondary ionization, entered as log */ 01267 strcpy( chAssertType[nAsserts] , "sc" ); 01268 /* say that this is sec ioniz */ 01269 strcpy( chAssertLineLabel[nAsserts], "sion" ); 01270 01271 /* now get rate, saved into assert quantity */ 01272 i = 5; 01273 AssertQuantity[nAsserts] = 01274 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01275 if( lgEOL ) 01276 { 01277 NoNumb(input.chOrgCard); 01278 } 01279 /* entered as log, make linear */ 01280 AssertQuantity[nAsserts] = pow(10., AssertQuantity[nAsserts] ); 01281 01282 /* no wavelength */ 01283 wavelength[nAsserts] = 0; 01284 01285 /* optional error, default available (cannot do before loop since we 01286 * do not know how many numbers are on line */ 01287 AssertError[nAsserts] = 01288 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01289 if( lgEOL ) 01290 { 01291 /* default error was set in define above */ 01292 AssertError[nAsserts] = DEF_ERROR; 01293 } 01294 01295 /* we want to print the log of eden, not linear value */ 01296 lgQuantityLog[nAsserts] = true; 01297 } 01298 01299 /* assert heating rate, erg/cm3/s, htot */ 01300 else if( nMatch("HTOT",input.chCARDCAPS ) ) 01301 { 01302 /* this flag will mean heating, entered as log */ 01303 strcpy( chAssertType[nAsserts] , "ht" ); 01304 01305 /* say that this is heating rate */ 01306 strcpy( chAssertLineLabel[nAsserts], "heat" ); 01307 01308 /* now get rate, saved into assert quantity */ 01309 i = 5; 01310 AssertQuantity[nAsserts] = 01311 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01312 if( lgEOL ) 01313 { 01314 NoNumb(input.chOrgCard); 01315 } 01316 /* entered as log, make linear */ 01317 AssertQuantity[nAsserts] = pow(10., AssertQuantity[nAsserts] ); 01318 01319 /* no wavelength */ 01320 wavelength[nAsserts] = 0; 01321 01322 /* optional error, default available (cannot do before loop since we 01323 * do not know how many numbers are on line */ 01324 AssertError[nAsserts] = 01325 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01326 if( lgEOL ) 01327 { 01328 /* default error was set in define above */ 01329 AssertError[nAsserts] = DEF_ERROR; 01330 } 01331 01332 /* we want to print the log of the heating, not linear value */ 01333 lgQuantityLog[nAsserts] = true; 01334 } 01335 01336 /* assert number of iterations per zone, a test of convergence */ 01337 else if( nMatch("ITRZ",input.chCARDCAPS ) ) 01338 { 01339 /* this flag will mean number of iterations per zone */ 01340 strcpy( chAssertType[nAsserts] , "iz" ); 01341 /* say that this is iterations per zone */ 01342 strcpy( chAssertLineLabel[nAsserts], "itrz" ); 01343 01344 /* now get quantity */ 01345 i = 5; 01346 AssertQuantity[nAsserts] = 01347 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01348 if( lgEOL ) 01349 { 01350 NoNumb(input.chOrgCard); 01351 } 01352 /* wavelength is meaningless */ 01353 wavelength[nAsserts] = 0; 01354 01355 /* optional error, default available (cannot do before loop since we 01356 * do not know how many numbers are on line */ 01357 AssertError[nAsserts] = 01358 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01359 if( lgEOL ) 01360 { 01361 /* default error was set in define above */ 01362 /* >>chng 04 mar 05, from default 5% to very small number, 01363 * since we will usually be doing upper limit */ 01364 AssertError[nAsserts] = 0.001; 01365 } 01366 } 01367 01368 /* assert electron density of the last zone */ 01369 else if( nMatch("EDEN",input.chCARDCAPS ) ) 01370 { 01371 /* this flag will mean electron density of the last zone */ 01372 strcpy( chAssertType[nAsserts] , "e " ); 01373 /* say that this is electron density */ 01374 strcpy( chAssertLineLabel[nAsserts], "eden" ); 01375 01376 /* now get electron density, which is a log */ 01377 i = 5; 01378 AssertQuantity[nAsserts] = 01379 pow(10., FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL) ); 01380 if( lgEOL ) 01381 { 01382 NoNumb(input.chOrgCard); 01383 } 01384 01385 /* optional error, default available (cannot do before loop since we 01386 * do not know how many numbers are on line */ 01387 AssertError[nAsserts] = 01388 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01389 if( lgEOL ) 01390 { 01391 /* default error was set in define above */ 01392 AssertError[nAsserts] = DEF_ERROR; 01393 } 01394 wavelength[nAsserts] = 0; 01395 01396 /* we want to print the log of eden, not linear value */ 01397 lgQuantityLog[nAsserts] = true; 01398 } 01399 01400 /* assert thickness or depth of model */ 01401 else if( nMatch("THIC",input.chCARDCAPS ) || nMatch("DEPT",input.chCARDCAPS ) ) 01402 { 01403 /* this flag will mean thickness or depth */ 01404 strcpy( chAssertType[nAsserts] , "th" ); 01405 /* say that this is thickness */ 01406 strcpy( chAssertLineLabel[nAsserts], "thic" ); 01407 01408 /* now get thickness, which is a log */ 01409 i = 5; 01410 AssertQuantity[nAsserts] = 01411 pow(10., FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL) ); 01412 if( lgEOL ) 01413 { 01414 NoNumb(input.chOrgCard); 01415 } 01416 01417 /* optional error, default available (cannot do before loop since we 01418 * do not know how many numbers are on line */ 01419 AssertError[nAsserts] = 01420 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01421 if( lgEOL ) 01422 { 01423 /* default error was set in define above */ 01424 AssertError[nAsserts] = DEF_ERROR; 01425 } 01426 wavelength[nAsserts] = 0; 01427 01428 /* we want to print the log of eden, not linear value */ 01429 lgQuantityLog[nAsserts] = true; 01430 } 01431 01432 /* assert outer radius of model */ 01433 else if( nMatch("RADI",input.chCARDCAPS ) ) 01434 { 01435 /* this flag will mean radius */ 01436 strcpy( chAssertType[nAsserts] , "ra" ); 01437 /* say that this is radius */ 01438 strcpy( chAssertLineLabel[nAsserts], "radi" ); 01439 01440 /* now get radius, which is a log */ 01441 i = 5; 01442 AssertQuantity[nAsserts] = 01443 pow(10., FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL) ); 01444 if( lgEOL ) 01445 { 01446 NoNumb(input.chOrgCard); 01447 } 01448 01449 /* optional error, default available (cannot do before loop since we 01450 * do not know how many numbers are on line */ 01451 AssertError[nAsserts] = 01452 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01453 if( lgEOL ) 01454 { 01455 /* default error was set in define above */ 01456 AssertError[nAsserts] = DEF_ERROR; 01457 } 01458 wavelength[nAsserts] = 0; 01459 01460 /* we want to print the log of radius, not linear value */ 01461 lgQuantityLog[nAsserts] = true; 01462 } 01463 01464 /* assert (probably upper limit to) number of iterations */ 01465 else if( nMatch("NITE",input.chCARDCAPS ) ) 01466 { 01467 /* this flag will mean number of iterations */ 01468 strcpy( chAssertType[nAsserts] , "Z " ); 01469 /* say that this is iteration numbering */ 01470 strcpy( chAssertLineLabel[nAsserts], "iter" ); 01471 01472 /* now get number of iterations, which will be saved into wavelength */ 01473 i = 5; 01474 wavelength[nAsserts] = 01475 (realnum)FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01476 if( lgEOL ) 01477 { 01478 NoNumb(input.chOrgCard); 01479 } 01480 /* copy it here although we will not use it */ 01481 AssertQuantity[nAsserts] = (double)wavelength[nAsserts]; 01482 01483 /* optional error, default available (cannot do before loop since we 01484 * do not know how many numbers are on line */ 01485 AssertError[nAsserts] = 01486 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01487 if( lgEOL ) 01488 { 01489 /* default error was set in define above */ 01490 AssertError[nAsserts] = DEF_ERROR; 01491 } 01492 } 01493 01494 /* assert terminal velocity, at end of calculation 01495 * input in km/s and saved that way, even though code uses cm/s */ 01496 else if( nMatch("VELO",input.chCARDCAPS ) ) 01497 { 01498 /* this flag will mean velocity */ 01499 strcpy( chAssertType[nAsserts] , "Ve" ); 01500 /* say that this is velocity */ 01501 strcpy( chAssertLineLabel[nAsserts], "vel " ); 01502 wavelength[nAsserts] = 0; 01503 01504 /* now get velocity */ 01505 i = 5; 01506 AssertQuantity[nAsserts] = 01507 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01508 if( lgEOL ) 01509 NoNumb(input.chOrgCard); 01510 01511 /* optional error, default available (cannot do before loop since we 01512 * do not know how many numbers are on line */ 01513 AssertError[nAsserts] = 01514 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 01515 if( lgEOL ) 01516 { 01517 /* default error was set in define above */ 01518 AssertError[nAsserts] = DEF_ERROR; 01519 } 01520 } 01521 /* assert nothing - a pacifier */ 01522 else if( nMatch("NOTH",input.chCARDCAPS ) ) 01523 { 01524 strcpy( chAssertType[nAsserts] , "NO" ); 01525 strcpy( chAssertLineLabel[nAsserts], "noth" ); 01526 wavelength[nAsserts] = 0; 01527 AssertQuantity[nAsserts] = 0.; 01528 AssertError[nAsserts] = DEF_ERROR; 01529 } 01530 else 01531 { 01532 /* did not recognize a command */ 01533 fprintf( ioQQQ, 01534 " Unrecognized command. The line image was==%s==\n", 01535 input.chOrgCard ); 01536 fprintf( ioQQQ, 01537 " The options I know about are: ionization, line, departure coefficient, map, column, " 01538 "temperature, nzone, csupre, htot, itrz, eden, thickness, niter, \n"); 01539 fprintf( ioQQQ, " Sorry.\n" ); 01540 cdEXIT(EXIT_FAILURE); 01541 } 01542 01543 /* increment number of asserts and confirm that limit not exceeded */ 01544 ++nAsserts; 01545 if( nAsserts >= NASSERTS ) 01546 { 01547 fprintf(ioQQQ, 01548 " ParseAssertResults: too many asserts, limit is NASSERT=%d\n", 01549 NASSERTS ); 01550 cdEXIT(EXIT_FAILURE); 01551 } 01552 return; 01553 } 01554 01555 /*============================================================================*/ 01556 /*lgCheckAsserts, checks asserts, last thing cloudy calls, returns true if all are ok, false if problems */ 01557 bool lgCheckAsserts( 01558 /* this is the file we will write this to, usually standard output, 01559 * but can be punch */ 01560 FILE * ioASSERT ) 01561 { 01562 double PredQuan[NASSERTS] , RelError[NASSERTS]; 01563 /* this will be true if the quantity was found, and false if not. Used to prevent 01564 * big botch flag when quantity not found (as in removed old he atom) */ 01565 bool lgFound[NASSERTS]; 01566 double relint , absint; 01567 bool lg1OK[NASSERTS]; 01568 long i,j; 01569 /* This structure is for reporting another close match for asserts of line 01570 * intensities only. The zeroth, first, and second elements for each assert are, 01571 * respectively, the first, second, and third matches the code finds, if any. 01572 * A negative number means no match. A positive number indicates the pointer 01573 * in the line stack of that match. */ 01574 long ipDisambiguate[NASSERTS][3]; 01575 long lgDisambiguate = false; 01576 char chLabelLoc[INPUT_LINE_LENGTH]; 01577 char chCaps[5], chFind[5]; 01578 realnum errorwave; 01579 01580 DEBUG_ENTRY( "lgCheckAsserts()" ); 01581 01582 /* this is a global variable in assertresults.h, and can be checked by 01583 * other routines to see if asserts are ok - (most runs will not use asserts) */ 01584 lgAssertsOK = true; 01585 01586 /* will be used if there were big botched asserts */ 01587 lgBigBotch = false; 01588 01589 /* the optimize*.in and stars_oppim*.in tests in the test suite include 01590 * asserts while optimizing. We do not want to test the asserts during 01591 * the optimization process, since we will find blown asserts and report 01592 * major problems. We do want to test asserts on the final model however. 01593 * Printout will usually be turned off in all except the final model, 01594 * so do not to the assert tests if we are optimizing but not printing */ 01595 if( !called.lgTalk && optimize.lgOptimize ) 01596 { 01597 /* just return */ 01598 return true; 01599 } 01600 01601 /*fprintf(ioQQQ , "DEBUG grid %li\n", optimize.nOptimiz );*/ 01602 01603 /* this will usually just return, but with table lines will check 01604 * existence of some lines */ 01605 if( lines_table() ) 01606 { 01607 lgBigBotch = true; 01608 lgAssertsOK = false; 01609 } 01610 01611 /* first check all asserts, there probably are none */ 01612 for(i=0; i<nAsserts; ++i ) 01613 { 01614 lg1OK[i] = true; 01615 PredQuan[i] = 0.; 01616 RelError[i] = 0.; 01617 ipDisambiguate[i][0] = -1; 01618 ipDisambiguate[i][1] = -1; 01619 ipDisambiguate[i][2] = -1; 01620 01621 /* this flag is set false if we don't find the requested quantity */ 01622 lgFound[i] = true; 01623 01624 /* which type of assert? */ 01625 /* is it intensity? */ 01626 if( strcmp( chAssertType[i] , "Lr")==0 ) 01627 { 01628 /* this is an intensity, get the line, returns false if could not find it */ 01629 ipDisambiguate[i][0] = cdLine( chAssertLineLabel[i] , wavelength[i] , &relint,&absint ); 01630 if( ipDisambiguate[i][0] <= 0 ) 01631 { 01632 fprintf( ioASSERT, 01633 " assert error: lgCheckAsserts could not find a line with label %s ", 01634 chAssertLineLabel[i] ); 01635 prt_wl( ioASSERT, wavelength[i] ); 01636 fprintf( ioASSERT, "\n" ); 01637 01638 fprintf( ioASSERT, 01639 " assert error: The \"punch line labels\" command is a good way to get a list of line labels.\n\n"); 01640 /* go to next line */ 01641 lg1OK[i] = false; 01642 RelError[i] = 100000.; 01643 PredQuan[i] = 0; 01644 lg1OK[i] = false; 01645 lgAssertsOK = false; 01646 lgFound[i] = false; 01647 continue; 01648 } 01649 else 01650 { 01651 /********* LINE DISAMBIGUATION *************/ 01652 /* Here we look for lines with same label and small wavelength 01653 * differences so that we can disambiguate below */ 01654 01655 /* change chLabel to all caps */ 01656 strcpy( chLabelLoc , chAssertLineLabel[i] ); 01657 /*cap4(chFind,chLabel);*/ 01658 cap4(chFind,chLabelLoc); 01659 01660 /* get the error associated with 4 significant figures */ 01661 errorwave = WavlenErrorGet( wavelength[i] ); 01662 01663 /* go through rest of line stack to look for close matches */ 01664 for( j=1; j < LineSave.nsum; j++ ) 01665 { 01666 /* don't bother with this one, we've already identified it. */ 01667 if( j==ipDisambiguate[i][0] ) 01668 continue; 01669 01670 /* change chLabel to all caps to be like input chALab */ 01671 cap4(chCaps , (char*)LineSv[j].chALab); 01672 01673 /* look for wavelengths within 3 error bars. 01674 * For example, for a line entered in Angstroms with 01675 * four significant figures, the error bar is 0.5 Ang. 01676 * So here we will find any lines within 1.5 Angstroms 01677 * of the 01678 * asserted wavelength. check wavelength and chLabel for a match */ 01679 if( fabs(LineSv[j].wavelength-wavelength[i]) < 3.f*errorwave ) 01680 { 01681 /* now see if labels agree */ 01682 if( strcmp(chCaps,chFind) == 0 ) 01683 { 01684 double relint1, absint1, current_error; 01685 01686 cdLine_ip( j, &relint1, &absint1 ); 01687 current_error = fabs(1. - relint1/AssertQuantity[i]); 01688 01689 if( current_error < 2.*AssertError[i] || 01690 current_error < 2.*fabs(RelError[i]) ) 01691 { 01692 lgDisambiguate = true; 01693 /* if second match (element 1) is already set, 01694 * this is third match (element 2). Set and break out. */ 01695 if( ipDisambiguate[i][1] > 0 ) 01696 { 01697 ipDisambiguate[i][2] = j; 01698 break; 01699 } 01700 else 01701 { 01702 ipDisambiguate[i][1] = j; 01703 } 01704 } 01705 } 01706 } 01707 } 01708 } 01709 01710 PredQuan[i] = relint; 01711 RelError[i] = 1. - PredQuan[i]/AssertQuantity[i]; 01712 } 01713 01714 /*this is line luminosity */ 01715 else if( strcmp( chAssertType[i] , "Ll")==0 ) 01716 { 01717 /* this is a luminosity, get the line, returns false if could not find it */ 01718 if( cdLine( chAssertLineLabel[i] , wavelength[i] , &relint,&absint )<=0 ) 01719 { 01720 fprintf( ioASSERT, 01721 " assert error: lgCheckAsserts could not find a line with label %s ", 01722 chAssertLineLabel[i] ); 01723 prt_wl( ioASSERT, wavelength[i] ); 01724 fprintf( ioASSERT, "\n" ); 01725 01726 fprintf( ioASSERT, 01727 " assert error: The \"punch line labels\" command is a good way to get a list of line labels.\n\n"); 01728 /* go to next line */ 01729 lg1OK[i] = false; 01730 RelError[i] = 10000000.; 01731 PredQuan[i] = 0; 01732 lg1OK[i] = false; 01733 lgFound[i] = false; 01734 lgAssertsOK = false; 01735 continue; 01736 } 01737 /* absint was returned as the log */ 01738 PredQuan[i] = pow(10.,absint); 01739 RelError[i] = 1. - PredQuan[i]/AssertQuantity[i]; 01740 01741 } 01742 else if( strcmp( chAssertType[i] , "hh" ) == 0) 01743 { 01744 double hfrac , hefrac; 01745 /* get H ionization fraction, returns false if could not find it */ 01746 if( cdIonFrac( 01747 /* four char string, null terminated, giving the element name */ 01748 "hydr", 01749 /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number */ 01750 1, 01751 /* will be fractional ionization */ 01752 &hfrac, 01753 /* how to weight the average, must be "VOLUME" or "RADIUS" */ 01754 "VOLUME" , 01755 /* do not want extra factor of density */ 01756 false) ) 01757 { 01758 fprintf( ioASSERT, 01759 " assert error: lgCheckAsserts could not find h ionization fraction \n"); 01760 lg1OK[i] = false; 01761 RelError[i] = 0; 01762 PredQuan[i] = 0; 01763 lgFound[i] = false; 01764 continue; 01765 } 01766 if( cdIonFrac( 01767 /* four char string, null terminated, giving the element name */ 01768 "heli", 01769 /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number */ 01770 1, 01771 /* will be fractional ionization */ 01772 &hefrac, 01773 /* how to weight the average, must be "VOLUME" or "RADIUS" */ 01774 "VOLUME" , 01775 /* do not want extra factor of density */ 01776 false) ) 01777 { 01778 fprintf( ioASSERT, 01779 " assert error: lgCheckAsserts could not find h ionization fraction \n"); 01780 lg1OK[i] = false; 01781 RelError[i] = 0; 01782 PredQuan[i] = 0; 01783 lgFound[i] = false; 01784 continue; 01785 } 01786 /* the helium hydrogen ionization correction factor */ 01787 PredQuan[i] = hefrac-hfrac; 01788 /* two icf's in difference, no need to div by mean since already on scale with unity */ 01789 RelError[i] = fabs(AssertQuantity[i] - (hefrac-hfrac) ); 01790 } 01791 01792 else if( strcmp( chAssertType[i] , "z " ) == 0) 01793 { 01794 /* this is (probably an upper limit) to the number of zones */ 01795 PredQuan[i] = (double)nzone; 01796 /* two integers in difference */ 01797 RelError[i] = 1. - PredQuan[i]/AssertQuantity[i]; 01798 } 01799 01800 else if( strcmp( chAssertType[i] , "or" ) == 0) 01801 { 01802 /* ortho to para ratio for large H2 molecule in last zone */ 01803 PredQuan[i] = h2.ortho_density / SDIV( h2.para_density ); 01804 01805 /* this is relative error */ 01806 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i]; 01807 } 01808 01809 else if( strcmp( chAssertType[i] , "g2" ) == 0) 01810 { 01811 /* check Jura rate, rate per vol that H2 forms on grain surfaces */ 01812 PredQuan[i] = gv.rate_h2_form_grains_used_total; 01813 /* this is relative error */ 01814 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i]; 01815 } 01816 01817 else if( strcmp( chAssertType[i] , "RM" ) == 0) 01818 { 01819 /* check Jura rate, rate per vol that H2 forms on grain surfaces */ 01820 PredQuan[i] = pressure.RadBetaMax; 01821 /* this is relative error */ 01822 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i]; 01823 } 01824 01825 else if( strcmp( chAssertType[i] , "pr" ) == 0) 01826 { 01827 /* standard deviation of the pressure */ 01828 double sumx=0., sumx2=0., average; 01829 long int n; 01830 /* do sums to form standard deviation */ 01831 for( n=0; n<nzone; n++ ) 01832 { 01833 sumx += struc.pressure[n]; 01834 sumx2 += POW2(struc.pressure[n]); 01835 } 01836 if( nzone>1 ) 01837 { 01838 /* this is average */ 01839 average = sumx/nzone; 01840 /* this is abs std */ 01841 sumx = sqrt( (sumx2-POW2(sumx)/nzone)/(nzone-1) ); 01842 /* save the relative std */ 01843 PredQuan[i] = sumx / average; 01844 } 01845 else 01846 { 01847 PredQuan[i] = 0.; 01848 } 01849 01850 /* two integers in difference */ 01851 RelError[i] = 0.; 01852 } 01853 01854 else if( strcmp( chAssertType[i] , "iz") == 0 ) 01855 { 01856 /* this is number of iterations per zone, a test of convergence properties */ 01857 if( nzone > 0 ) 01858 PredQuan[i] = (double)(conv.nTotalIoniz-conv.nTotalIoniz_start)/(double)(nzone); 01859 else 01860 /* something big so assert will botch. */ 01861 PredQuan[i] = 1e10; 01862 01863 /* this is relative error */ 01864 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i]; 01865 } 01866 01867 else if( strcmp( chAssertType[i] , "e ") == 0 ) 01868 { 01869 /* this is electron density of the last zone */ 01870 PredQuan[i] = dense.eden; 01871 /* this is relative error */ 01872 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i]; 01873 } 01874 01875 else if( strcmp( chAssertType[i] , "th") == 0 ) 01876 { 01877 /* this is thickness */ 01878 PredQuan[i] = radius.depth; 01879 /* this is relative error */ 01880 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i]; 01881 } 01882 01883 else if( strcmp( chAssertType[i] , "ra") == 0 ) 01884 { 01885 /* this is thickness */ 01886 PredQuan[i] = radius.Radius; 01887 /* this is relative error */ 01888 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i]; 01889 } 01890 01891 else if( strcmp( chAssertType[i] , "Ve") == 0 ) 01892 { 01893 /* this is final velocity of wind in km/s (code uses cm/s) */ 01894 PredQuan[i] = wind.windv/1e5; 01895 /* this is relative error */ 01896 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i]; 01897 } 01898 01899 else if( strcmp( chAssertType[i] , "NO") == 0 ) 01900 { 01901 /* this is nothing */ 01902 PredQuan[i] = 0; 01903 /* this is relative error */ 01904 RelError[i] = 0.; 01905 } 01906 01907 else if( strcmp( chAssertType[i] , "sc") == 0 ) 01908 { 01909 /* this is secondary ionization rate */ 01910 PredQuan[i] = secondaries.csupra[ipHYDROGEN][0]; 01911 /* this is relative error */ 01912 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i]; 01913 } 01914 01915 else if( strcmp( chAssertType[i] , "ht") == 0 ) 01916 { 01917 /* this is heating rate */ 01918 PredQuan[i] = thermal.htot; 01919 /* this is relative error */ 01920 RelError[i] = 1. - PredQuan[i] / AssertQuantity[i]; 01921 } 01922 01923 else if( strcmp( chAssertType[i] , "Z ") == 0 ) 01924 { 01925 /* this is (probably an upper limit) to the number of iterations */ 01926 PredQuan[i] = (double)iteration; 01927 /* two integers in difference */ 01928 RelError[i] = (double)(wavelength[i] - iteration); 01929 /* uncertainty in difference is zero */ 01930 AssertError[i] = 0.; 01931 } 01932 01933 else if( strcmp( chAssertType[i] , "CB") == 0 ) 01934 { 01935 long int nISOCaseB = (long)Param[i][0]; 01936 long int nelemCaseB = (long)Param[i][1]; 01937 char chElemLabelCaseB[5]; 01938 /* generate label to find element, as "H 1" */ 01939 strcpy( chElemLabelCaseB , elementnames.chElementSym[nelemCaseB] ); 01940 char chNumb[4]; 01941 sprintf( chNumb , "%2li" , nelemCaseB+1-nISOCaseB ); 01942 strcat( chElemLabelCaseB , chNumb ); 01943 01944 int iCase = 1; 01945 RelError[i] = 0.; 01946 long nHighestPrinted = StatesElem[nISOCaseB][nelemCaseB][iso.numPrintLevels[nISOCaseB][nelemCaseB]-1].n; 01947 if( nISOCaseB == ipH_LIKE ) 01948 { 01949 fprintf(ioASSERT," Species nHi nLo Wl Computed Asserted error\n"); 01950 /* limit of 10 is because that is all we printed and saved in prt_lines_hydro */ 01951 for( long int ipLo=1+iCase; ipLo< MIN2(10,nHighestPrinted-1); ++ipLo ) 01952 { 01953 for( long int ipHi=ipLo+1; ipHi< MIN2(25,nHighestPrinted); ++ipHi ) 01954 { 01955 /* assert the line */ 01956 realnum wl = atmdat.WaveLengthCaseB[nelemCaseB][ipHi][ipLo]; 01957 /* range option to restrict wavelength coverage */ 01958 if( wl < Param[i][2] || wl > Param[i][3] ) 01959 continue; 01960 01961 double relint , absint,CBrelint , CBabsint; 01962 /* find the predicted line intensity */ 01963 cdLine( chAssertLineLabel[i] , wl , &CBrelint , &CBabsint ); 01964 if( CBrelint < Param[i][4] ) 01965 continue; 01966 CBabsint = pow( 10., CBabsint ); 01967 double error; 01968 /* now find the Case B intensity - may not all be present */ 01969 if( cdLine( chElemLabelCaseB , wl , &relint , &absint ) >0) 01970 { 01971 absint = pow( 10., absint ); 01972 error = (CBabsint - absint)/MAX2(CBabsint , absint); 01973 double RelativeError = fabs(error) / AssertError[i]; 01974 /* start of line, flag problems */ 01975 if( RelativeError < 1. ) 01976 { 01977 if( RelativeError < 0.25 ) 01978 { 01979 fprintf( ioASSERT, " ChkAssert "); 01980 } 01981 else if( RelativeError < 0.50 ) 01982 { 01983 fprintf( ioASSERT, " ChkAssert - "); 01984 } 01985 else if( RelativeError < 0.75 ) 01986 { 01987 fprintf( ioASSERT, " ChkAssert -- "); 01988 } 01989 else if( RelativeError < 0.90 ) 01990 { 01991 fprintf( ioASSERT, " ChkAssert --- "); 01992 } 01993 else if( RelativeError < 0.95 ) 01994 { 01995 fprintf( ioASSERT, " ChkAssert ---- "); 01996 } 01997 else if( RelativeError < 0.98 ) 01998 { 01999 fprintf( ioASSERT, " ChkAssert ----- "); 02000 } 02001 else 02002 { 02003 fprintf( ioASSERT, " ChkAssert ------ "); 02004 } 02005 02006 } 02007 else 02008 { 02009 fprintf( ioASSERT, " ChkAssert botch>>"); 02010 } 02011 fprintf(ioASSERT," %s %3li %3li ", 02012 chElemLabelCaseB , ipHi , ipLo ); 02013 prt_wl(ioASSERT, wl ); 02014 fprintf(ioASSERT," %.2e %.2e %10.3f", 02015 absint , CBabsint , error ); 02016 } 02017 else 02018 TotalInsanity(); 02019 if( fabs(error) > AssertError[i] ) 02020 fprintf(ioASSERT , " botch \n"); 02021 else 02022 fprintf(ioASSERT , "\n"); 02023 02024 PredQuan[i] = 0; 02025 AssertQuantity[i] = 0.; 02026 RelError[i] = MAX2( RelError[i] , fabs(error) ); 02027 02028 /* save sum which we will report later */ 02029 assertresults.SumErrorCaseAssert += RelError[i]; 02030 ++assertresults.nSumErrorCaseAssert; 02031 02032 } 02033 } 02034 fprintf(ioASSERT,"\n"); 02035 } 02036 else if( nISOCaseB == ipHE_LIKE ) 02037 { 02038 if( !dense.lgElmtOn[ipHELIUM] ) 02039 { 02040 fprintf(ioQQQ,"PROBLEM assert case B for a He is requested but He is not " 02041 "included.\n"); 02042 fprintf(ioQQQ,"Do not turn off He if you want to assert its spectrum.\n"); 02043 cdEXIT(EXIT_FAILURE); 02044 } 02045 # if 0 02046 # define N_CASEB_HEI 11 02047 realnum CaseBWlHeI[N_CASEB_HEI]= 02048 { 10830.f, 3889.f, 3188.f, 5016.f, 3965.f, 7065.f, 5876.f, 4471.f, 02049 4026.f, 6678.f, 4922.f }; 02050 # endif 02051 /* do He I as special case */ 02052 fprintf(ioASSERT," Wl Computed Asserted error\n"); 02053 for( long int ipLine=0; ipLine< atmdat.nCaseBHeI ; ++ipLine ) 02054 { 02055 /* assert the line */ 02056 realnum wl = atmdat.CaseBWlHeI[ipLine]; 02057 /* range option to restrict wavelength coverage */ 02058 if( wl < Param[i][2] || wl > Param[i][3] ) 02059 continue; 02060 double relint , absint,CBrelint , CBabsint; 02061 cdLine( chAssertLineLabel[i] , wl , &CBrelint , &CBabsint ); 02062 if( CBrelint < Param[i][4] ) 02063 continue; 02064 CBabsint = pow( 10., CBabsint ); 02065 double error; 02066 if( cdLine( chElemLabelCaseB , wl , &relint , &absint ) >0) 02067 { 02068 absint = pow( 10., absint ); 02069 error = (CBabsint - absint)/MAX2(CBabsint , absint); 02070 double RelativeError = fabs(error) / AssertError[i]; 02071 /* start of line, flag problems */ 02072 if( RelativeError < 1. ) 02073 { 02074 if( RelativeError < 0.25 ) 02075 { 02076 fprintf( ioASSERT, " ChkAssert "); 02077 } 02078 else if( RelativeError < 0.50 ) 02079 { 02080 fprintf( ioASSERT, " ChkAssert - "); 02081 } 02082 else if( RelativeError < 0.75 ) 02083 { 02084 fprintf( ioASSERT, " ChkAssert -- "); 02085 } 02086 else if( RelativeError < 0.90 ) 02087 { 02088 fprintf( ioASSERT, " ChkAssert --- "); 02089 } 02090 else if( RelativeError < 0.95 ) 02091 { 02092 fprintf( ioASSERT, " ChkAssert ---- "); 02093 } 02094 else if( RelativeError < 0.98 ) 02095 { 02096 fprintf( ioASSERT, " ChkAssert ----- "); 02097 } 02098 else 02099 { 02100 fprintf( ioASSERT, " ChkAssert ------ "); 02101 } 02102 02103 } 02104 else 02105 { 02106 fprintf( ioASSERT, " ChkAssert botch>>"); 02107 } 02108 prt_wl(ioASSERT, wl ); 02109 fprintf(ioASSERT," %.2e %.2e %10.3f", 02110 absint , CBabsint , error ); 02111 } 02112 else 02113 TotalInsanity(); 02114 if( fabs(error) > AssertError[i] ) 02115 fprintf(ioASSERT , " botch \n"); 02116 else 02117 fprintf(ioASSERT , "\n"); 02118 02119 PredQuan[i] = 0; 02120 AssertQuantity[i] = 0.; 02121 RelError[i] = MAX2( RelError[i] , fabs(error) ); 02122 02123 /* save sum which we will report later */ 02124 assertresults.SumErrorCaseAssert += RelError[i]; 02125 ++assertresults.nSumErrorCaseAssert; 02126 } 02127 fprintf(ioASSERT,"\n"); 02128 } 02129 else 02130 TotalInsanity(); 02131 } 02132 02133 else if( strcmp( chAssertType[i] , "D ") == 0 ) 02134 { 02135 long ipZ , n; 02136 /* this is departure coefficient for H-like ion given by wavelength */ 02137 /* stored number was element number on C scale */ 02138 ipZ = (long)wavelength[i]; 02139 if( !dense.lgElmtOn[ipZ] ) 02140 { 02141 fprintf(ioQQQ,"PROBLEM asserted element %ld is not turned on!\n",ipZ); 02142 PredQuan[i] = 0; 02143 RelError[i] = 0.; 02144 } 02145 else 02146 { 02147 RelError[i] = 0.; 02148 PredQuan[i] = 0.; 02149 for( n=(long)AssertQuantity2[i]; n<iso.numPrintLevels[ipH_LIKE][ipZ]; ++n ) 02150 { 02151 PredQuan[i] += iso.DepartCoef[ipH_LIKE][ipZ][n]; 02152 /* relerror is the largest deviation from unity for any single level*/ 02153 RelError[i] = MAX2( RelError[i] , iso.DepartCoef[ipH_LIKE][ipZ][n] -1.); 02154 RelError[i] = MAX2( RelError[i] , 1. - iso.DepartCoef[ipH_LIKE][ipZ][n] ); 02155 } 02156 PredQuan[i] /= (double)(iso.numPrintLevels[ipH_LIKE][ipZ]); 02157 RelError[i] /= (double)(iso.numPrintLevels[ipH_LIKE][ipZ]); 02158 } 02159 } 02160 02161 /* departure coefficients for something on he-like seq */ 02162 else if( strcmp( chAssertType[i] , "De") == 0 ) 02163 { 02164 long ipZ , n; 02165 /* this is departure coefficient for He-like ion given by wavelength */ 02166 /* stored number was element number on C scale */ 02167 ipZ = (long)wavelength[i]; 02168 if( !dense.lgElmtOn[ipZ] ) 02169 { 02170 fprintf(ioQQQ,"PROBLEM asserted element %ld is not turned on!\n",ipZ); 02171 PredQuan[i] = 0.; 02172 RelError[i] = 0.; 02173 } 02174 else 02175 { 02176 RelError[i] = 0.; 02177 PredQuan[i] = 0.; 02178 for( n=(long)AssertQuantity2[i]; n<iso.numPrintLevels[ipHE_LIKE][ipZ]; ++n ) 02179 { 02180 double relerror; 02181 relerror = 0.; 02182 PredQuan[i] += iso.DepartCoef[ipHE_LIKE][ipZ][n]; 02183 /* relerror is the largest deviation from unity for any single level*/ 02184 relerror = iso.DepartCoef[ipHE_LIKE][ipZ][n] -1.; 02185 relerror = MAX2( relerror , 1. - iso.DepartCoef[ipHE_LIKE][ipZ][n] ); 02186 RelError[i] = MAX2( RelError[i] , relerror / AssertQuantity[i] ); 02187 } 02188 PredQuan[i] /= (double)(iso.numPrintLevels[ipHE_LIKE][ipZ]); 02189 /* >>chng 03 mar 13, had not div by num levels */ 02190 RelError[i] /= (double)(iso.numPrintLevels[ipHE_LIKE][ipZ]); 02191 } 02192 } 02193 02194 /* this is FeII departure coefficient */ 02195 else if( strcmp( chAssertType[i] , "d ") == 0 ) 02196 { 02197 double BigError , StdDev; 02198 /* this is departure coefficient for FeII */ 02199 AssertFeIIDep( &PredQuan[i] , &BigError , &StdDev ); 02200 /* there are many missing A's in the FeII ion (no atomic data) so only the 02201 * average is relevant now, RefError returned above is single largest 02202 * excursion from unity */ 02203 RelError[i] = StdDev; 02204 } 02205 02206 /* this is H- departure coefficient */ 02207 else if( strcmp( chAssertType[i] , "d-") == 0 ) 02208 { 02209 PredQuan[i] = hmi.hmidep; 02210 RelError[i] = fabs( hmi.hmidep - 1.); 02211 } 02212 02213 /* this would be ionization fraction */ 02214 else if( (strcmp( chAssertType[i] , "f ") == 0) || 02215 (strcmp( chAssertType[i] , "F ") == 0) ) 02216 { 02217 char chWeight[7]; 02218 if( strcmp( chAssertType[i] , "F ") == 0 ) 02219 { 02220 strcpy( chWeight , "VOLUME" ); 02221 } 02222 else 02223 { 02224 /* this is default case, Fraction over radius */ 02225 strcpy( chWeight , "RADIUS" ); 02226 } 02227 /* get ionization fraction, returns false if could not find it */ 02228 if( cdIonFrac( 02229 /* four char string, null terminated, giving the element name */ 02230 chAssertLineLabel[i], 02231 /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number */ 02232 (long)wavelength[i], 02233 /* will be fractional ionization */ 02234 &relint, 02235 /* how to weight the average, must be "VOLUME" or "RADIUS" */ 02236 chWeight , 02237 /* do not want extra factor of density */ 02238 false) ) 02239 { 02240 fprintf( ioASSERT, 02241 " assert error: lgCheckAsserts could not find a line with label %s %f \n", 02242 chAssertLineLabel[i] , wavelength[i] ); 02243 /* go to next line */ 02244 lg1OK[i] = false; 02245 RelError[i] = 0; 02246 PredQuan[i] = 0; 02247 lgFound[i] = false; 02248 continue; 02249 } 02250 /* this is ionization fraction */ 02251 PredQuan[i] = relint; 02252 RelError[i] = 1. - PredQuan[i]/AssertQuantity[i]; 02253 } 02254 02255 /* this would be column density of several molecules */ 02256 else if( strcmp( chAssertType[i] , "cd") == 0) 02257 { 02258 /* for H2 column density - total or for a state? */ 02259 if( (strcmp( chAssertLineLabel[i], "H2 " ) == 0) && (Param[i][0] >= 0.) ) 02260 { 02261 /* this branch get state specific column density */ 02262 /* get total H2 column density */ 02263 if( (relint = cdH2_colden( (long)Param[i][0] , (long)Param[i][1] ) ) < 0. ) 02264 { 02265 fprintf(ioQQQ," PROBLEM lgCheckAsserts did not find v=%li, J=%li for H2 column density.\n", 02266 (long)Param[i][0] , (long)Param[i][1] ); 02267 lg1OK[i] = false; 02268 RelError[i] = 0; 02269 PredQuan[i] = 0; 02270 lgFound[i] = false; 02271 continue; 02272 } 02273 } 02274 else 02275 { 02276 /* get ionization fraction, returns 0 if all ok */ 02277 if( cdColm( 02278 /* four char string, null terminated, giving the element name */ 02279 chAssertLineLabel[i], 02280 /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number, 02281 * zero for molecule*/ 02282 (long)wavelength[i], 02283 /* will be fractional ionization */ 02284 &relint) ) 02285 { 02286 fprintf( ioASSERT, 02287 " assert error: lgCheckAsserts could not find a molecule with label %s %f \n", 02288 chAssertLineLabel[i] , wavelength[i] ); 02289 /* go to next line */ 02290 lg1OK[i] = false; 02291 RelError[i] = 0; 02292 PredQuan[i] = 0; 02293 lgFound[i] = false; 02294 continue; 02295 } 02296 } 02297 /* this is ionization fraction */ 02298 PredQuan[i] = relint; 02299 RelError[i] = 1. - PredQuan[i]/AssertQuantity[i]; 02300 } 02301 02302 /* this would be molecular fraction of CO or H2 */ 02303 else if( (strcmp( chAssertType[i] , "MF") == 0) || (strcmp( chAssertType[i] , "mf") == 0) ) 02304 { 02305 /* get molecular fraction, returns 0 if all ok */ 02306 relint = 0.; 02307 if( strcmp( chAssertLineLabel[i], "H2 ") ==0) 02308 { 02309 /* get total H2 column density */ 02310 if( cdColm("H2 " , 0, 02311 /* will be fractional ionization */ 02312 &relint) ) 02313 TotalInsanity(); 02314 02315 relint = relint / colden.colden[ipCOL_HTOT]; 02316 } 02317 else 02318 { 02319 fprintf( ioASSERT, 02320 " assert error: lgCheckAsserts could not find a molecule with label %s %f \n", 02321 chAssertLineLabel[i] , wavelength[i] ); 02322 /* go to next line */ 02323 lg1OK[i] = false; 02324 RelError[i] = 0; 02325 PredQuan[i] = 0; 02326 continue; 02327 } 02328 /* this is ionization fraction */ 02329 PredQuan[i] = relint; 02330 RelError[i] = 1. - PredQuan[i]/AssertQuantity[i]; 02331 } 02332 02333 /* check heating/cooling at some temperature in a thermal map */ 02334 else if( strcmp( chAssertType[i] , "mh") == 0 || strcmp( chAssertType[i] , "mc") == 0) 02335 { 02336 /* check heating or cooling (stored in error array) at temperature in assert results */ 02337 /* check that map was done, and arrays have nmap elements */ 02338 if( hcmap.nMapAlloc == 0 ) 02339 { 02340 /* this happens if map not done and space for h/c not allocated */ 02341 fprintf( ioASSERT, 02342 " assert error: lgCheckAsserts cannot check map since map not done.\n"); 02343 /* go to next line */ 02344 lg1OK[i] = false; 02345 RelError[i] = 0; 02346 PredQuan[i] = 0; 02347 continue; 02348 } 02349 /* now check that requested temperature is within the range of the map we computed */ 02350 if( AssertQuantity2[i]< hcmap.temap[0] || AssertQuantity2[i]> hcmap.temap[hcmap.nmap-1] ) 02351 { 02352 fprintf( ioASSERT, 02353 " assert error: lgCheckAsserts cannot check map since temperature not within range.\n"); 02354 /* go to next line */ 02355 lg1OK[i] = false; 02356 RelError[i] = 0; 02357 PredQuan[i] = 0; 02358 continue; 02359 } 02360 02361 /* we should have valid data - find closest temperature >- requested temperature */ 02362 j = 0; 02363 while( AssertQuantity2[i]>hcmap.temap[j]*1.001 && j < hcmap.nmap ) 02364 { 02365 ++j; 02366 } 02367 02368 /* j points to correct cell in heating cooling array */ 02369 /* we will not interpolate, just use this value, and clobber te to prove it*/ 02370 if( strcmp( chAssertType[i] , "mh") == 0 ) 02371 { 02372 /* heating */ 02373 PredQuan[i] = hcmap.hmap[j]; 02374 strcpy(chAssertLineLabel[i],"MapH" ); 02375 } 02376 else if( strcmp( chAssertType[i] , "mc") == 0) 02377 { 02378 /* cooling */ 02379 PredQuan[i] = hcmap.cmap[j]; 02380 strcpy(chAssertLineLabel[i],"MapC" ); 02381 } 02382 RelError[i] = 1. - PredQuan[i]/AssertQuantity[i]; 02383 } 02384 02385 /* this will be an average temperature */ 02386 else if( (strcmp( chAssertType[i] , "t ") == 0) || 02387 (strcmp( chAssertType[i] , "T ") == 0) ) 02388 { 02389 char chWeight[7]; 02390 if( strcmp( chAssertType[i] , "T ") == 0 ) 02391 { 02392 strcpy( chWeight , "VOLUME" ); 02393 } 02394 else 02395 { 02396 /* this is default case, Fraction over radius */ 02397 strcpy( chWeight , "RADIUS" ); 02398 } 02399 02400 /* options are average Te for ion, temp at ill face, or temp for grain */ 02401 if( strcmp( chAssertLineLabel[i], "GTem" ) == 0 ) 02402 { 02403 long nd; 02404 /* the minus one is because the grain types are counted from one, 02405 * but stuffed into the c array, that counts from zero */ 02406 nd = (long)(wavelength[i]-1); 02407 if( nd < 0 || nd >= gv.nBin ) { 02408 fprintf( ioQQQ, "Illegal grain number found: %f\n" , wavelength[i] ); 02409 fprintf( ioQQQ, "Use 1 for first grain that is turned on, " ); 02410 fprintf( ioQQQ, "2 for second, etc....\n" ); 02411 fprintf( ioQQQ, "Old style grain numbers are not valid anymore !!\n" ); 02412 cdEXIT(EXIT_FAILURE); 02413 } 02414 ASSERT( nd>= 0 ); 02415 relint = gv.bin[nd]->avdust/radius.depth_x_fillfac; 02416 } 02417 else if( strcmp( chAssertLineLabel[i], "face" ) == 0 ) 02418 { 02419 /* this is the temperature at the illuminated face */ 02420 relint = struc.testr[0]; 02421 } 02422 else 02423 { 02424 /* get temperature, returns false if could not find it */ 02425 if( cdTemp( 02426 /* four char string, null terminated, giving the element name */ 02427 chAssertLineLabel[i], 02428 /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number */ 02429 (long)wavelength[i], 02430 /* will be mean temperatue */ 02431 &relint, 02432 /* how to weight the average, must be "VOLUME" or "RADIUS" */ 02433 chWeight ) ) 02434 { 02435 fprintf( ioASSERT, 02436 " assert error: lgCheckAsserts could not find an ion with label %s ion %li \n", 02437 chAssertLineLabel[i] , (long)wavelength[i] ); 02438 /* go to next line */ 02439 lg1OK[i] = false; 02440 RelError[i] = 0; 02441 PredQuan[i] = 0; 02442 lgFound[i] = false; 02443 continue; 02444 } 02445 } 02446 /* this is the temperature */ 02447 PredQuan[i] = relint; 02448 RelError[i] = 1. - PredQuan[i]/AssertQuantity[i]; 02449 } 02450 02451 /* this would be grain potential in volt */ 02452 else if( strcmp( chAssertType[i], "gp" ) == 0 ) 02453 { 02454 /* the minus one is because the grain types are counted from one, 02455 * but stuffed into the c array, that counts from zero */ 02456 long nd = (long)(wavelength[i]-1); 02457 if( nd < 0 || nd >= gv.nBin ) { 02458 fprintf( ioQQQ, "Illegal grain number found: %g\n" , wavelength[i] ); 02459 fprintf( ioQQQ, "Use 1 for first grain that is turned on, " ); 02460 fprintf( ioQQQ, "2 for second, etc....\n" ); 02461 fprintf( ioQQQ, "Old style grain numbers are not valid anymore !!\n" ); 02462 cdEXIT(EXIT_FAILURE); 02463 } 02464 02465 /* get average grain potential in volt, always averaged over radius */ 02466 PredQuan[i] = gv.bin[nd]->avdpot/radius.depth_x_fillfac; 02467 /* actually absolute error, potential can be zero! */ 02468 RelError[i] = AssertQuantity[i] - PredQuan[i]; 02469 } 02470 02471 else 02472 { 02473 fprintf( ioASSERT, 02474 " assert error: lgCheckAsserts received an insane chAssertType=%s, impossible\n", 02475 chAssertType[i] ); 02476 ShowMe(); 02477 cdEXIT(EXIT_FAILURE); 02478 } 02479 02480 if( chAssertLimit[i] == '=' ) 02481 { 02482 /* predicted quantity should be within error of expected */ 02483 if( fabs(RelError[i]) > AssertError[i] ) 02484 { 02485 lg1OK[i] = false; 02486 lgAssertsOK = false; 02487 } 02488 } 02489 else if( chAssertLimit[i] == '<' ) 02490 { 02491 /* expected is an upper limit, so PredQuan/AssertQuantity should 02492 * be less than one, and so RelError should be positive 02493 * >>chng 04 feb 14,from <= to < - limit is really < not <=, 02494 * in case of iterations or zones, iter < iterations would not 02495 * trigger false when iter == iterations */ 02496 if( RelError[i] <= 0.-AssertError[i]) 02497 { 02498 lg1OK[i] = false; 02499 lgAssertsOK = false; 02500 } 02501 } 02502 else if( chAssertLimit[i] == '>' ) 02503 { 02504 /* expected is a lower limit, so PredQuan/AssertQuantity should 02505 * be greater than one, and so RelError should be negative */ 02506 if( RelError[i] >= 0.+AssertError[i]) 02507 { 02508 lg1OK[i] = false; 02509 lgAssertsOK = false; 02510 } 02511 } 02512 } 02513 02514 /* only print summary if we are talking */ 02515 if( called.lgTalk && nAsserts>0 ) 02516 { 02517 time_t now; 02518 02519 /* First disambiguate any line identifications */ 02520 if( lgDisambiguate ) 02521 { 02522 /* change significant figures of WL for this printout */ 02523 long sigfigsav = LineSave.sig_figs; 02524 double relint1, relint2, absint1; 02525 02526 LineSave.sig_figs = 6; 02527 02528 fprintf( ioASSERT, "=============Line Disambiguation=======================================================\n" ); 02529 fprintf( ioASSERT, " Wavelengths || Intensities \n" ); 02530 fprintf( ioASSERT, "Label line match1 match2 match3 || asserted match1 match2 match3\n" ); 02531 02532 for( i=0; i<nAsserts; ++i ) 02533 { 02534 if( ipDisambiguate[i][1] > 0 ) 02535 { 02536 fprintf( ioASSERT , "%4s ", chAssertLineLabel[i] ); 02537 prt_wl( ioASSERT , wavelength[i] ); 02538 fprintf( ioASSERT , " " ); 02539 prt_wl( ioASSERT , LineSv[ipDisambiguate[i][0]].wavelength ); 02540 fprintf( ioASSERT , " " ); 02541 prt_wl( ioASSERT , LineSv[ipDisambiguate[i][1]].wavelength ); 02542 fprintf( ioASSERT , " " ); 02543 if( ipDisambiguate[i][2] > 0 ) 02544 { 02545 prt_wl( ioASSERT , LineSv[ipDisambiguate[i][2]].wavelength ); 02546 cdLine_ip( ipDisambiguate[i][2], &relint2, &absint1 ); 02547 } 02548 else 02549 { 02550 fprintf( ioASSERT , "--------" ); 02551 relint2 = 0.0; 02552 } 02553 fprintf( ioASSERT , " ||" ); 02554 02555 cdLine_ip( ipDisambiguate[i][1], &relint1, &absint1 ); 02556 02557 if( lgPrtSciNot ) 02558 { 02559 fprintf( ioASSERT , " %10.3e %10.3e %10.3e %10.3e\n", 02560 AssertQuantity[i], 02561 PredQuan[i] , 02562 relint1, 02563 relint2 ); 02564 } 02565 else 02566 { 02567 fprintf( ioASSERT , " %10.4f %10.4f %10.4f %10.4f\n", 02568 AssertQuantity[i], 02569 PredQuan[i] , 02570 relint1, 02571 relint2 ); 02572 } 02573 } 02574 } 02575 fprintf( ioASSERT, "\n" ); 02576 02577 /* revert to original significant figures */ 02578 LineSave.sig_figs = sigfigsav; 02579 } 02580 02581 /* write start of title and version number of code */ 02582 fprintf( ioASSERT, "=============Results of asserts: Cloudy %s ", t_version::Inst().chVersion ); 02583 02584 /* usually print date and time info - do not if "no times" command entered, 02585 * which set this flag false */ 02586 if( prt.lgPrintTime ) 02587 { 02588 /* now add date of this run */ 02589 now = time(NULL); 02590 /* now print this time at the end of the string. the system put cr at the end of the string */ 02591 fprintf(ioASSERT,"%s", ctime(&now) ); 02592 } 02593 else 02594 { 02595 fprintf(ioASSERT,"\n" ); 02596 } 02597 02598 if( lgAssertsOK ) 02599 { 02600 fprintf( ioASSERT, " No errors were found. Summary follows.\n"); 02601 } 02602 else 02603 { 02604 fprintf( ioASSERT, " Errors were found. Summary follows.\n"); 02605 } 02606 02607 fprintf( ioASSERT, 02608 " Label line computed asserted Rel Err Set err\n"); 02609 /* now print a summary */ 02610 for( i=0; i<nAsserts; ++i ) 02611 { 02612 double prtPredQuan, prtAssertQuantity; 02613 /* this is option to print log of quantity rather than linear. default is 02614 * linear, and log only for a few such as ionization to see small numbers */ 02615 if( lgQuantityLog[i] ) 02616 { 02617 prtPredQuan = log10( MAX2( SMALLDOUBLE , PredQuan[i] ) ); 02618 prtAssertQuantity = log10( MAX2( SMALLDOUBLE , AssertQuantity[i] ) ); 02619 } 02620 else 02621 { 02622 prtPredQuan = PredQuan[i]; 02623 prtAssertQuantity = AssertQuantity[i]; 02624 } 02625 /* start of line, flag problems */ 02626 if( lg1OK[i] ) 02627 { 02628 /* the ChkAssert is a unique label so that we can grep on all output 02629 * and see what happened, without picking up input stream */ 02630 double relative = fabs(RelError[i] / SDIV( AssertError[i])); 02631 02632 if( relative < 0.25 || chAssertLimit[i] != '=' ) 02633 { 02634 fprintf( ioASSERT, " ChkAssert "); 02635 } 02636 else if( relative < 0.50 ) 02637 { 02638 fprintf( ioASSERT, " ChkAssert - "); 02639 } 02640 else if( relative < 0.75 ) 02641 { 02642 fprintf( ioASSERT, " ChkAssert -- "); 02643 } 02644 else if( relative < 0.90 ) 02645 { 02646 fprintf( ioASSERT, " ChkAssert --- "); 02647 } 02648 else if( relative < 0.95 ) 02649 { 02650 fprintf( ioASSERT, " ChkAssert ---- "); 02651 } 02652 else if( relative < 0.98 ) 02653 { 02654 fprintf( ioASSERT, " ChkAssert ----- "); 02655 } 02656 else 02657 { 02658 fprintf( ioASSERT, " ChkAssert ------ "); 02659 } 02660 02661 } 02662 else 02663 { 02664 fprintf( ioASSERT, " ChkAssert botch>>"); 02665 } 02666 if( strcmp( chAssertType[i] , "Ll")==0 || ( strcmp( chAssertType[i] , "Lr")==0 ) ) 02667 { 02668 /* special formatting for the emission lines */ 02669 fprintf( ioASSERT , "%4s ", 02670 chAssertLineLabel[i] ); 02671 02672 prt_wl( ioASSERT , wavelength[i] ); 02673 02674 if( lgPrtSciNot ) 02675 { 02676 fprintf( ioASSERT , " %10.3e %c %10.3e %7.3f %7.3f ", 02677 prtPredQuan , 02678 chAssertLimit[i] , 02679 prtAssertQuantity , 02680 RelError[i] , 02681 AssertError[i]); 02682 } 02683 else 02684 { 02685 fprintf( ioASSERT , " %10.4f %c %10.4f %7.3f %7.3f ", 02686 prtPredQuan , 02687 chAssertLimit[i] , 02688 prtAssertQuantity , 02689 RelError[i] , 02690 AssertError[i]); 02691 } 02692 02693 { 02694 enum {DEBUG_LOC=false}; 02695 02696 if( DEBUG_LOC ) 02697 { 02698 long ipHi, ipLo; 02699 errorwave = WavlenErrorGet( wavelength[i] ); 02700 02701 for( ipLo = ipHe1s1S; ipLo <= iso.numLevels_local[ipHE_LIKE][ipHELIUM] - 02702 iso.nCollapsed_local[ipHE_LIKE][ipHELIUM]; ipLo ++ ) 02703 { 02704 for( ipHi = ipLo+1; ipHi < iso.numLevels_local[ipHE_LIKE][ipHELIUM] - 02705 iso.nCollapsed_local[ipHE_LIKE][ipHELIUM]; ipHi++ ) 02706 { 02707 if( fabs(Transitions[ipHE_LIKE][ipHELIUM][ipHi][ipLo].WLAng - wavelength[i]) 02708 < errorwave && ipLo!=ipHe2p3P0 && ipLo!=ipHe2p3P1 ) 02709 { 02710 fprintf( ioQQQ, "\t%li %li %li %li %li %li", 02711 StatesElem[ipHE_LIKE][ipHELIUM][ipHi].n, 02712 StatesElem[ipHE_LIKE][ipHELIUM][ipHi].l, 02713 StatesElem[ipHE_LIKE][ipHELIUM][ipHi].S, 02714 StatesElem[ipHE_LIKE][ipHELIUM][ipLo].n, 02715 StatesElem[ipHE_LIKE][ipHELIUM][ipLo].l, 02716 StatesElem[ipHE_LIKE][ipHELIUM][ipLo].S ); 02717 break; 02718 } 02719 } 02720 } 02721 } 02722 } 02723 02724 } 02725 else 02726 { 02727 fprintf( ioASSERT , "%4s %6i %10.4f %c %10.4f %7.3f %7.3f ", 02728 chAssertLineLabel[i] , 02729 (int)wavelength[i] , 02730 prtPredQuan , 02731 chAssertLimit[i] , 02732 prtAssertQuantity , 02733 RelError[i] , 02734 AssertError[i]); 02735 } 02736 /* if botched and the botch is > 3 sigma, say BIG BOTCH, 02737 * the lg1OK is needed since some tests (number of zones, etc) 02738 * are limits, not the quantity, and if limit is large the 02739 * miss will be big too */ 02740 02741 /* >>chng 02 nov 27, added lgFound so don't get big botch when line simply missing */ 02742 if( !lg1OK[i] && (fabs(RelError[i]) > 3.*AssertError[i]) && lgFound[i] ) 02743 { 02744 fprintf( ioASSERT , " <<BIG BOTCH!!\n"); 02745 lgBigBotch = true; 02746 } 02747 else 02748 { 02749 fprintf( ioASSERT , "\n"); 02750 } 02751 } 02752 fprintf( ioASSERT , " \n"); 02753 02754 /* NB - in following, perl scripts detect these strings - be careful if they 02755 * are changed - the scripts may no longer detect major errors */ 02756 if( !lgAssertsOK && lgBigBotch ) 02757 { 02758 /* there were big botches */ 02759 fprintf( ioASSERT, " BIG BOTCHED ASSERTS!!! Big Botched Asserts!!! \n"); 02760 } 02761 else if( !lgAssertsOK ) 02762 { 02763 fprintf( ioASSERT, " BOTCHED ASSERTS!!! Botched Asserts!!! \n"); 02764 } 02765 02766 if( assertresults.nSumErrorCaseAssert>0 ) 02767 { 02768 fprintf(ioASSERT,"\n The mean of the %li assert Case A, B relative " 02769 "residuals is %.2f\n\n" , 02770 assertresults.nSumErrorCaseAssert, 02771 assertresults.SumErrorCaseAssert /assertresults.nSumErrorCaseAssert ); 02772 } 02773 02774 /* explain how we were compiled, but only if printing time */ 02775 if( prt.lgPrintTime ) 02776 { 02777 fprintf( ioQQQ, " %s\n\n", t_version::Inst().chInfo ); 02778 } 02779 } 02780 return lgAssertsOK; 02781 }