Changeset 2023 for branches/newmole
- Timestamp:
- 05/05/08 05:09:17 (8 months ago)
- Location:
- branches/newmole/source
- Files:
-
- 97 modified
-
assertresults.h (modified) (1 diff)
-
assert_results.cpp (modified) (13 diffs)
-
atmdat.h (modified) (1 diff)
-
atmdat_char_tran.cpp (modified) (2 diffs)
-
atmdat_HS_caseb.cpp (modified) (2 diffs)
-
atmdat_lines_setup.cpp (modified) (1 diff)
-
atom_feii.cpp (modified) (2 diffs)
-
atom_pop5.cpp (modified) (1 diff)
-
atom_seq_beryllium.cpp (modified) (1 diff)
-
cdgetlinelist.cpp (modified) (2 diffs)
-
container_classes.h (modified) (1 diff)
-
cont_createmesh.cpp (modified) (2 diffs)
-
cont_ffun.cpp (modified) (2 diffs)
-
cont_gammas.cpp (modified) (2 diffs)
-
cont_setintensity.cpp (modified) (1 diff)
-
conv_base.cpp (modified) (2 diffs)
-
conv_init_solution.cpp (modified) (3 diffs)
-
cool_carb.cpp (modified) (1 diff)
-
cool_oxyg.cpp (modified) (1 diff)
-
cpu.cpp (modified) (4 diffs)
-
date.h (modified) (1 diff)
-
eden_sum.cpp (modified) (2 diffs)
-
gammas.h (modified) (2 diffs)
-
grains.cpp (modified) (2 diffs)
-
grains_qheat.cpp (modified) (2 diffs)
-
heat_punch.cpp (modified) (1 diff)
-
heat_sum.cpp (modified) (2 diffs)
-
helike_cs.cpp (modified) (1 diff)
-
helike_einsta.cpp (modified) (6 diffs)
-
helike_recom.cpp (modified) (2 diffs)
-
highen.cpp (modified) (1 diff)
-
hydroeinsta.cpp (modified) (3 diffs)
-
hydrolevel.cpp (modified) (2 diffs)
-
hydrot2low.cpp (modified) (2 diffs)
-
hydro_bauman.cpp (modified) (2 diffs)
-
hydro_recom.cpp (modified) (2 diffs)
-
hydro_vs_rates.cpp (modified) (1 diff)
-
init_coreload_postparse.cpp (modified) (2 diffs)
-
init_defaults_preparse.cpp (modified) (1 diff)
-
ion_carbo.cpp (modified) (2 diffs)
-
ion_nitro.cpp (modified) (1 diff)
-
ion_oxyge.cpp (modified) (2 diffs)
-
ion_recomb.cpp (modified) (2 diffs)
-
ion_silic.cpp (modified) (7 diffs)
-
ion_solver.cpp (modified) (1 diff)
-
ion_sulph.cpp (modified) (1 diff)
-
ion_trim.cpp (modified) (2 diffs)
-
iso.h (modified) (9 diffs)
-
iso_collide.cpp (modified) (2 diffs)
-
iso_continuum_lower.cpp (modified) (1 diff)
-
iso_cool.cpp (modified) (3 diffs)
-
iso_create.cpp (modified) (12 diffs)
-
iso_level.cpp (modified) (2 diffs)
-
iso_photo.cpp (modified) (2 diffs)
-
iso_radiative_recomb.cpp (modified) (7 diffs)
-
iso_solve.cpp (modified) (8 diffs)
-
iter_startend.cpp (modified) (1 diff)
-
lines.h (modified) (1 diff)
-
mole_h2.cpp (modified) (4 diffs)
-
mole_h2_io.cpp (modified) (3 diffs)
-
nemala2.cpp (modified) (1 diff)
-
opacity.h (modified) (1 diff)
-
optimize_subplx.cpp (modified) (2 diffs)
-
parse_abundances.cpp (modified) (1 diff)
-
parse_atom_iso.cpp (modified) (3 diffs)
-
parse_commands.cpp (modified) (9 diffs)
-
parse_compile.cpp (modified) (4 diffs)
-
parse_crashdo.cpp (modified) (2 diffs)
-
parse_dlaw.cpp (modified) (3 diffs)
-
parse_element.cpp (modified) (2 diffs)
-
parse_fluc.cpp (modified) (3 diffs)
-
parse_globule.cpp (modified) (1 diff)
-
parse_grain.cpp (modified) (1 diff)
-
parse_hden.cpp (modified) (3 diffs)
-
parse_punch.cpp (modified) (6 diffs)
-
pressure_total.cpp (modified) (1 diff)
-
prt_comment.cpp (modified) (2 diffs)
-
prt_final.cpp (modified) (1 diff)
-
prt_linepres.cpp (modified) (1 diff)
-
prt_lines_helium.cpp (modified) (5 diffs)
-
prt_lines_lv1_na_ar.cpp (modified) (1 diff)
-
prt_met.cpp (modified) (1 diff)
-
punch.h (modified) (1 diff)
-
punch_do.cpp (modified) (2 diffs)
-
punch_fits.cpp (modified) (3 diffs)
-
punch_line.cpp (modified) (1 diff)
-
radius_increment.cpp (modified) (1 diff)
-
rt_line_all.cpp (modified) (5 diffs)
-
rt_line_driving.cpp (modified) (1 diff)
-
rt_line_one.cpp (modified) (1 diff)
-
rt_stark.cpp (modified) (2 diffs)
-
rt_tau_inc.cpp (modified) (9 diffs)
-
sanity_check.cpp (modified) (1 diff)
-
service.cpp (modified) (1 diff)
-
sil.h (modified) (1 diff)
-
stars.cpp (modified) (4 diffs)
-
zero.cpp (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/newmole/source/assertresults.h
r1739 r2023 27 27 EXTERN bool lgAssertsOK , lgBigBotch, lgPrtSciNot; 28 28 29 EXTERN struct t_assertresults { 30 double SumErrorCaseAssert; 31 long int nSumErrorCaseAssert; 32 33 } assertresults; 34 29 35 #endif /* _ASSERTRESULTS_H_ */ -
branches/newmole/source/assert_results.cpp
r1942 r2023 48 48 static char **chAssertLineLabel; 49 49 50 /* element to assert H-like case B */51 static char chElemLabelCaseB[5];52 static long int nelemCaseB , nISOCaseB;53 54 50 /* this will be = for equal, < or > for limit */ 55 51 static char *chAssertLimit; … … 122 118 123 119 /* these are a pair of optional parameters */ 124 Param[i] = ((double *)MALLOC( 2*sizeof(double )));120 Param[i] = ((double *)MALLOC(5*sizeof(double ))); 125 121 } 126 122 … … 457 453 AssertError[nAsserts] = DEF_ERROR; 458 454 AssertQuantity[nAsserts] = 0; 459 strcpy( chAssertLineLabel[nAsserts] , "Ca B" );460 455 wavelength[nAsserts] = 0.; 461 456 462 /* assert case b for H or He */ 457 /* faint option - do not test line if relative intensity is less 458 * than entered value */ 459 if( (i=nMatch("FAINT",input.chCARDCAPS ))>0 ) 460 { 461 Param[nAsserts][4] = 462 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 463 if( lgEOL ) 464 { 465 /* did not get 2 numbers */ 466 fprintf(ioQQQ," The assert Case B faint option must have a number," 467 " the relative intensity of the fainest line to assert.\n"); 468 cdEXIT(EXIT_FAILURE); 469 } 470 /* number is log if <= 0 */ 471 if( Param[nAsserts][4]<=0. ) 472 Param[nAsserts][4] = pow(10., Param[nAsserts][4] ); 473 } 474 else 475 { 476 /* use default - include everything*/ 477 Param[nAsserts][4] = SMALLFLOAT; 478 } 479 480 /* range option - to limit check on a certain wavelength range */ 481 if( (i=nMatch("RANG",input.chCARDCAPS ))>0 ) 482 { 483 Param[nAsserts][2] = 484 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 485 Param[nAsserts][3] = 486 FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL); 487 if( lgEOL ) 488 { 489 /* did not get 2 numbers */ 490 fprintf(ioQQQ," The assert Case B range option must have two numbers," 491 " the lower and upper limit to the wavelengths in Angstroms.\n"); 492 fprintf(ioQQQ," There must be a total of three numbers on the line," 493 " the relative error followed by the lower and upper limits to the " 494 "wavelength in Angstroms.\n"); 495 cdEXIT(EXIT_FAILURE); 496 } 497 if( Param[nAsserts][2]>Param[nAsserts][3]) 498 { 499 /* make sure in increasing order */ 500 double sav = Param[nAsserts][3]; 501 Param[nAsserts][3] = Param[nAsserts][2]; 502 Param[nAsserts][2] = sav; 503 } 504 } 505 else 506 { 507 /* use default - include everything*/ 508 Param[nAsserts][2] = 0.; 509 Param[nAsserts][3] = 1e30; 510 } 511 /* assert case b for H - O checking against Hummer & Storey tables */ 463 512 if( nMatch("H-LI",input.chCARDCAPS ) ) 464 513 { … … 478 527 cdEXIT(EXIT_FAILURE); 479 528 } 480 nelemCaseB = nelem; 481 nISOCaseB = 0; 482 strcpy( chElemLabelCaseB , elementnames.chElementSym[nelem] ); 483 char chNumb[4]; 484 sprintf( chNumb , "%2li" , nelem+1 ); 485 strcat( chElemLabelCaseB , chNumb ); 529 Param[nAsserts][0] = ipH_LIKE; 530 Param[nAsserts][1] = nelem; 531 /* generate string to find simple prediction, as in "Ca B" */ 532 strcpy( chAssertLineLabel[nAsserts], "Ca "); 533 if( nMatch("CASE A ",input.chCARDCAPS ) ) 534 strcat( chAssertLineLabel[nAsserts] , "A"); 535 else 536 strcat( chAssertLineLabel[nAsserts] , "B"); 486 537 } 487 538 else if( nMatch("HE-L",input.chCARDCAPS ) ) 488 539 { 489 540 /* He-like - only helium itself */ 490 nelemCaseB = ipHELIUM; 491 nISOCaseB = 0; 492 strcpy( chElemLabelCaseB , "He 1" ); 493 /* beyond reach of tables */ 494 fprintf(ioQQQ, "assert He-like case b not implemented yet %s\n", 495 input.chCARDCAPS ); 496 cdEXIT(EXIT_FAILURE); 541 Param[nAsserts][0] = ipHE_LIKE; 542 Param[nAsserts][1] = ipHELIUM; 543 544 strcpy( chAssertLineLabel[nAsserts] , "+Col"); 497 545 } 498 546 } … … 1477 1525 strcpy( chAssertLineLabel[nAsserts], "noth" ); 1478 1526 wavelength[nAsserts] = 0; 1479 AssertQuantity[nAsserts] = 0.; ;1527 AssertQuantity[nAsserts] = 0.; 1480 1528 AssertError[nAsserts] = DEF_ERROR; 1481 1529 } … … 1885 1933 else if( strcmp( chAssertType[i] , "CB") == 0 ) 1886 1934 { 1935 long int nISOCaseB = (long)Param[i][0]; 1936 long int nelemCaseB = (long)Param[i][1]; 1937 char chElemLabelCaseB[5]; 1938 /* generate label to find element, as "H 1" */ 1939 strcpy( chElemLabelCaseB , elementnames.chElementSym[nelemCaseB] ); 1940 char chNumb[4]; 1941 sprintf( chNumb , "%2li" , nelemCaseB+1-nISOCaseB ); 1942 strcat( chElemLabelCaseB , chNumb ); 1943 1887 1944 int iCase = 1; 1888 1945 RelError[i] = 0.; 1889 fprintf(ioQQQ," Species nHi nLo Wl Computed Asserted error\n");1890 1946 long nHighestPrinted = StatesElem[nISOCaseB][nelemCaseB][iso.numPrintLevels[nISOCaseB][nelemCaseB]-1].n; 1891 for( long int ipLo=1+iCase; ipLo< MIN2(24,nHighestPrinted-1); ++ipLo ) 1892 { 1893 for( long int ipHi=ipLo+1; ipHi< MIN2(25,nHighestPrinted); ++ipHi ) 1947 if( nISOCaseB == ipH_LIKE ) 1948 { 1949 fprintf(ioASSERT," Species nHi nLo Wl Computed Asserted error\n"); 1950 /* limit of 10 is because that is all we printed and saved in prt_lines_hydro */ 1951 for( long int ipLo=1+iCase; ipLo< MIN2(10,nHighestPrinted-1); ++ipLo ) 1952 { 1953 for( long int ipHi=ipLo+1; ipHi< MIN2(25,nHighestPrinted); ++ipHi ) 1954 { 1955 /* assert the line */ 1956 realnum wl = atmdat.WaveLengthCaseB[nelemCaseB][ipHi][ipLo]; 1957 /* range option to restrict wavelength coverage */ 1958 if( wl < Param[i][2] || wl > Param[i][3] ) 1959 continue; 1960 1961 double relint , absint,CBrelint , CBabsint; 1962 /* find the predicted line intensity */ 1963 cdLine( chAssertLineLabel[i] , wl , &CBrelint , &CBabsint ); 1964 if( CBrelint < Param[i][4] ) 1965 continue; 1966 CBabsint = pow( 10., CBabsint ); 1967 double error; 1968 /* now find the Case B intensity - may not all be present */ 1969 if( cdLine( chElemLabelCaseB , wl , &relint , &absint ) >0) 1970 { 1971 absint = pow( 10., absint ); 1972 error = (CBabsint - absint)/MAX2(CBabsint , absint); 1973 double RelativeError = fabs(error) / AssertError[i]; 1974 /* start of line, flag problems */ 1975 if( RelativeError < 1. ) 1976 { 1977 if( RelativeError < 0.25 ) 1978 { 1979 fprintf( ioASSERT, " ChkAssert "); 1980 } 1981 else if( RelativeError < 0.50 ) 1982 { 1983 fprintf( ioASSERT, " ChkAssert - "); 1984 } 1985 else if( RelativeError < 0.75 ) 1986 { 1987 fprintf( ioASSERT, " ChkAssert -- "); 1988 } 1989 else if( RelativeError < 0.90 ) 1990 { 1991 fprintf( ioASSERT, " ChkAssert --- "); 1992 } 1993 else if( RelativeError < 0.95 ) 1994 { 1995 fprintf( ioASSERT, " ChkAssert ---- "); 1996 } 1997 else if( RelativeError < 0.98 ) 1998 { 1999 fprintf( ioASSERT, " ChkAssert ----- "); 2000 } 2001 else 2002 { 2003 fprintf( ioASSERT, " ChkAssert ------ "); 2004 } 2005 2006 } 2007 else 2008 { 2009 fprintf( ioASSERT, " ChkAssert botch>>"); 2010 } 2011 fprintf(ioASSERT," %s %3li %3li ", 2012 chElemLabelCaseB , ipHi , ipLo ); 2013 prt_wl(ioASSERT, wl ); 2014 fprintf(ioASSERT," %.2e %.2e %10.3f", 2015 absint , CBabsint , error ); 2016 } 2017 else 2018 TotalInsanity(); 2019 if( fabs(error) > AssertError[i] ) 2020 fprintf(ioASSERT , " botch \n"); 2021 else 2022 fprintf(ioASSERT , "\n"); 2023 2024 PredQuan[i] = 0; 2025 AssertQuantity[i] = 0.; 2026 RelError[i] = MAX2( RelError[i] , fabs(error) ); 2027 2028 /* save sum which we will report later */ 2029 assertresults.SumErrorCaseAssert += RelError[i]; 2030 ++assertresults.nSumErrorCaseAssert; 2031 2032 } 2033 } 2034 fprintf(ioASSERT,"\n"); 2035 } 2036 else if( nISOCaseB == ipHE_LIKE ) 2037 { 2038 if( !dense.lgElmtOn[ipHELIUM] ) 2039 { 2040 fprintf(ioQQQ,"PROBLEM assert case B for a He is requested but He is not " 2041 "included.\n"); 2042 fprintf(ioQQQ,"Do not turn off He if you want to assert its spectrum.\n"); 2043 cdEXIT(EXIT_FAILURE); 2044 } 2045 # if 0 2046 # define N_CASEB_HEI 11 2047 realnum CaseBWlHeI[N_CASEB_HEI]= 2048 { 10830.f, 3889.f, 3188.f, 5016.f, 3965.f, 7065.f, 5876.f, 4471.f, 2049 4026.f, 6678.f, 4922.f }; 2050 # endif 2051 /* do He I as special case */ 2052 fprintf(ioASSERT," Wl Computed Asserted error\n"); 2053 for( long int ipLine=0; ipLine< atmdat.nCaseBHeI ; ++ipLine ) 1894 2054 { 1895 2055 /* assert the line */ 1896 realnum wl = atmdat.WaveLengthCaseB[nelemCaseB][ipHi][ipLo]; 2056 realnum wl = atmdat.CaseBWlHeI[ipLine]; 2057 /* range option to restrict wavelength coverage */ 2058 if( wl < Param[i][2] || wl > Param[i][3] ) 2059 continue; 1897 2060 double relint , absint,CBrelint , CBabsint; 1898 cdLine( "Ca B" , wl , &CBrelint , &CBabsint ); 2061 cdLine( chAssertLineLabel[i] , wl , &CBrelint , &CBabsint ); 2062 if( CBrelint < Param[i][4] ) 2063 continue; 1899 2064 CBabsint = pow( 10., CBabsint ); 1900 2065 double error; … … 1903 2068 absint = pow( 10., absint ); 1904 2069 error = (CBabsint - absint)/MAX2(CBabsint , absint); 1905 fprintf(ioQQQ,"ChkAssert %s %3li %3li ", 1906 chElemLabelCaseB , ipHi , ipLo ); 1907 prt_wl(ioQQQ, wl ); 1908 fprintf(ioQQQ," %.2e %.2e %10.3f", 2070 double RelativeError = fabs(error) / AssertError[i]; 2071 /* start of line, flag problems */ 2072 if( RelativeError < 1. ) 2073 { 2074 if( RelativeError < 0.25 ) 2075 { 2076 fprintf( ioASSERT, " ChkAssert "); 2077 } 2078 else if( RelativeError < 0.50 ) 2079 { 2080 fprintf( ioASSERT, " ChkAssert - "); 2081 } 2082 else if( RelativeError < 0.75 ) 2083 { 2084 fprintf( ioASSERT, " ChkAssert -- "); 2085 } 2086 else if( RelativeError < 0.90 ) 2087 { 2088 fprintf( ioASSERT, " ChkAssert --- "); 2089 } 2090 else if( RelativeError < 0.95 ) 2091 { 2092 fprintf( ioASSERT, " ChkAssert ---- "); 2093 } 2094 else if( RelativeError < 0.98 ) 2095 { 2096 fprintf( ioASSERT, " ChkAssert ----- "); 2097 } 2098 else 2099 { 2100 fprintf( ioASSERT, " ChkAssert ------ "); 2101 } 2102 2103 } 2104 else 2105 { 2106 fprintf( ioASSERT, " ChkAssert botch>>"); 2107 } 2108 prt_wl(ioASSERT, wl ); 2109 fprintf(ioASSERT," %.2e %.2e %10.3f", 1909 2110 absint , CBabsint , error ); 1910 2111 } … … 1912 2113 TotalInsanity(); 1913 2114 if( fabs(error) > AssertError[i] ) 1914 fprintf(io QQQ , " BOTCH\n");2115