Navigation

← Previous Change
Next Change →

assert_results.cpp

Timestamp:

05/10/08 09:03:02 (2 months ago)

Author:

peter

Message:

Merging all changes from mainline upto r2033, except r1902.

Files:

branches/c08_branch/source/assert_results.cpp (modified) (14 diffs)

Legend:

: Unmodified
: Added
: Removed
: Modified
: Copied
: Moved

branches/c08_branch/source/assert_results.cpp

r1896	r2034
10	10	#include "iso.h"
11	11	#include "called.h"
	12	#include "atmdat.h"
12	13	#include "hcmap.h"
13	14	#include "thermal.h"
…	…
117	118
118	119	/* these are a pair of optional parameters */
119		Param[i] = ((double )MALLOC(2sizeof(double )));
	120	Param[i] = ((double )MALLOC(5sizeof(double )));
120	121	}
121	122
…	…
344	345
345	346	/* assert line "LINE" -- key is ine_ since linear option appears on some commands */
346		else if( nMatch("INE ",input.chCARDCAPS ) )
	347	else if( nMatch(" LINE ",input.chCARDCAPS ) )
347	348	{
348	349	if( nMatch("LUMI",input.chCARDCAPS ) \|\| nMatch("INTE",input.chCARDCAPS ))
…	…
440	441
441	442	/* assert departure coefficients */
	443	else if( nMatch("CASE",input.chCARDCAPS ) )
	444	{
	445	/* this is Case B for some element */
	446	strcpy( chAssertType[nAsserts] , "CB" );
	447	i = 5;
	448	/* this is relative error */
	449	AssertError[nAsserts] =
	450	FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
	451	if( lgEOL )
	452	/* default error was set in define above */
	453	AssertError[nAsserts] = DEF_ERROR;
	454	AssertQuantity[nAsserts] = 0;
	455	wavelength[nAsserts] = 0.;
	456
	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 */
	512	if( nMatch("H-LI",input.chCARDCAPS ) )
	513	{
	514	/* H-like - now get an element */
	515	if( (nelem = GetElem( input.chCARDCAPS )) < 0 )
	516	{
	517	/* no name found */
	518	fprintf(ioQQQ, "assert H-like case B did not find an element on this line, sorry %s\n",
	519	input.chCARDCAPS );
	520	cdEXIT(EXIT_FAILURE);
	521	}
	522	if( nelem>7 )
	523	{
	524	/* beyond reach of tables */
	525	fprintf(ioQQQ, "assert H-like cannot do elements heavier than O, sorry %s\n",
	526	input.chCARDCAPS );
	527	cdEXIT(EXIT_FAILURE);
	528	}
	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");
	537	}
	538	else if( nMatch("HE-L",input.chCARDCAPS ) )
	539	{
	540	/* He-like - only helium itself */
	541	Param[nAsserts][0] = ipHE_LIKE;
	542	Param[nAsserts][1] = ipHELIUM;
	543
	544	strcpy( chAssertLineLabel[nAsserts] , "+Col");
	545	}
	546	}
442	547	else if( nMatch("DEPA",input.chCARDCAPS ) )
443	548	{
…	…
447	552	FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
448	553	if( lgEOL )
449		{
450	554	NoNumb(input.chOrgCard);
451		}
452	555
453	556	/* this is relative error, max departure from unity of any level or std */
…	…
455	558	FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
456	559	if( lgEOL )
457		{
458	560	/* default error was set in define above */
459	561	AssertError[nAsserts] = DEF_ERROR;
460		}
461	562
462	563	/* H-like key means do one of the hydrogenic ions */
…	…
1424	1525	strcpy( chAssertLineLabel[nAsserts], "noth" );
1425	1526	wavelength[nAsserts] = 0;
1426		AssertQuantity[nAsserts] = 0.;;
	1527	AssertQuantity[nAsserts] = 0.;
1427	1528	AssertError[nAsserts] = DEF_ERROR;
1428	1529	}
…	…
1830	1931	}
1831	1932
	1933	else if( strcmp( chAssertType[i] , "CB") == 0 )
	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
	1944	int iCase = 1;
	1945	RelError[i] = 0.;
	1946	long nHighestPrinted = StatesElem[nISOCaseB][nelemCaseB][iso.numPrintLevels[nISOCaseB][nelemCaseB]-1].n;
	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 )
	2054	{
	2055	/* assert the line */
	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;
	2060	double relint , absint,CBrelint , CBabsint;
	2061	cdLine( chAssertLineLabel[i] , wl , &CBrelint , &CBabsint );
	2062	if( CBrelint < Param[i][4] )
	2063	continue;
	2064	CBabsint = pow( 10., CBabsint );
	2065	double error;
	2066	if( cdLine( chElemLabelCaseB , wl , &relint , &absint ) >0)
	2067	{
	2068	absint = pow( 10., absint );
	2069	error = (CBabsint - absint)/MAX2(CBabsint , absint);
	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",
	2110	absint , CBabsint , error );
	2111	}
	2112	else
	2113	TotalInsanity();
	2114	if( fabs(error) > AssertError[i] )
	2115	fprintf(ioASSERT , " botch \n");
	2116	else
	2117	fprintf(ioASSERT , "\n");
	2118
	2119	PredQuan[i] = 0;
	2120	AssertQuantity[i] = 0.;
	2121	RelError[i] = MAX2( RelError[i] , fabs(error) );
	2122
	2123	/* save sum which we will report later */
	2124	assertresults.SumErrorCaseAssert += RelError[i];
	2125	++assertresults.nSumErrorCaseAssert;
	2126	}
	2127	fprintf(ioASSERT,"\n");
	2128	}
	2129	else
	2130	TotalInsanity();
	2131	}
	2132
1832	2133	else if( strcmp( chAssertType[i] , "D ") == 0 )
1833	2134	{
…	…
1838	2139	if( !dense.lgElmtOn[ipZ] )
1839	2140	{
1840		fprintf(ioQQQ," asserted element %ld is not turned on!\n",ipZ);
	2141	fprintf(ioQQQ,"PROBLEM asserted element %ld is not turned on!\n",ipZ);
1841	2142	PredQuan[i] = 0;
1842	2143	RelError[i] = 0.;
…	…
1867	2168	if( !dense.lgElmtOn[ipZ] )
1868	2169	{
1869		fprintf(ioQQQ," asserted element %ld is not turned on!\n",ipZ);
	2170	fprintf(ioQQQ,"PROBLEM asserted element %ld is not turned on!\n",ipZ);
1870	2171	PredQuan[i] = 0.;
1871	2172	RelError[i] = 0.;
…	…
2307	2608	" Label line computed asserted Rel Err Set err\n");
2308	2609	/* now print a summary */
2309		~~/lint -e771 vars conceivably not initialized /~~
2310	2610	for( i=0; i<nAsserts; ++i )
2311	2611	{
…	…
2392	2692
2393	2693	{
2394		~~/@-redef@/~~
2395	2694	enum {DEBUG_LOC=false};
2396		~~/@+redef@/~~
2397	2695
2398	2696	if( DEBUG_LOC )
…	…
2452	2750	}
2453	2751	}
2454		~~/lint +e771 vars conceivably not initialized /~~
2455	2752	fprintf( ioASSERT , " \n");
2456	2753
…	…
2467	2764	}
2468	2765
	2766	if( assertresults.nSumErrorCaseAssert>0 )
	2767	{
	2768	fprintf(ioASSERT,"\n The mean of the %li assert Case A, B relative "
	2769	"residuals is %.2f\n\n" ,
	2770	assertresults.nSumErrorCaseAssert,
	2771	assertresults.SumErrorCaseAssert /assertresults.nSumErrorCaseAssert );
	2772	}
	2773
2469	2774	/* explain how we were compiled, but only if printing time */
2470	2775	if( prt.lgPrintTime )

Navigation

Changeset 2034 for branches/c08_branch/source/assert_results.cpp

Legend:

branches/c08_branch/source/assert_results.cpp

Download in other formats: