Changeset 2023

Show
Ignore:
Timestamp:
05/05/08 05:09:17 (1 week ago)
Author:
rjrw
Message:

Merged from trunk r1941:2022

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • branches/newmole/source/assertresults.h

    r1739 r2023  
    2727EXTERN bool lgAssertsOK , lgBigBotch, lgPrtSciNot; 
    2828 
     29EXTERN struct t_assertresults { 
     30        double SumErrorCaseAssert; 
     31        long int nSumErrorCaseAssert; 
     32 
     33}       assertresults; 
     34 
    2935#endif /* _ASSERTRESULTS_H_ */ 

  • branches/newmole/source/assert_results.cpp

    r1942 r2023  
    4848static char **chAssertLineLabel; 
    4949 
    50 /* element to assert H-like case B */ 
    51 static char chElemLabelCaseB[5]; 
    52 static long int nelemCaseB , nISOCaseB; 
    53  
    5450/* this will be = for equal, < or > for limit */ 
    5551static char *chAssertLimit; 
    …… 
    122118 
    123119                        /* these are a pair of optional parameters */ 
    124                         Param[i] = ((double *)MALLOC(2*sizeof(double ))); 
     120                        Param[i] = ((double *)MALLOC(5*sizeof(double ))); 
    125121                } 
    126122 
    …… 
    457453                        AssertError[nAsserts] = DEF_ERROR; 
    458454                AssertQuantity[nAsserts] = 0; 
    459                 strcpy( chAssertLineLabel[nAsserts] , "Ca B" );  
    460455                wavelength[nAsserts] = 0.; 
    461456 
    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 */ 
    463512                if( nMatch("H-LI",input.chCARDCAPS ) ) 
    464513                { 
    …… 
    478527                                cdEXIT(EXIT_FAILURE); 
    479528                        } 
    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"); 
    486537                } 
    487538                else if( nMatch("HE-L",input.chCARDCAPS ) ) 
    488539                { 
    489540                        /* 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"); 
    497545                } 
    498546        } 
    …… 
    14771525                strcpy( chAssertLineLabel[nAsserts], "noth" ); 
    14781526                wavelength[nAsserts] = 0; 
    1479                 AssertQuantity[nAsserts] = 0.;; 
     1527                AssertQuantity[nAsserts] = 0.; 
    14801528                AssertError[nAsserts] = DEF_ERROR; 
    14811529        } 
    …… 
    18851933                else if( strcmp( chAssertType[i] , "CB") == 0 ) 
    18861934                { 
     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 
    18871944                        int iCase = 1; 
    18881945                        RelError[i] = 0.; 
    1889                         fprintf(ioQQQ,"        Species  nHi nLo  Wl  Computed  Asserted       error\n"); 
    18901946                        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 ) 
    18942054                                { 
    18952055                                        /* 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; 
    18972060                                        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; 
    18992064                                        CBabsint = pow( 10., CBabsint ); 
    19002065                                        double error; 
    …… 
    19032068                                                absint = pow( 10., absint ); 
    19042069                                                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",  
    19092110                                                        absint , CBabsint , error ); 
    19102111                                        } 
    …… 
    19122113                                                TotalInsanity(); 
    19132114                                        if( fabs(error) > AssertError[i]  ) 
    1914                                                 fprintf(ioQQQ , " BOTCH \n"); 
     2115                                                fprintf(ioASSERT , " botch \n"); 
    19152116                                        else 
    1916                                                 fprintf(ioQQQ , "\n"); 
     2117                                                fprintf(ioASSERT , "\n"); 
    19172118 
    19182119                                        PredQuan[i] = 0; 
    19192120                                        AssertQuantity[i] = 0.; 
    19202121                                        RelError[i] = MAX2( RelError[i] , fabs(error) ); 
     2122 
     2123                                        /* save sum which we will report later */ 
     2124                                        assertresults.SumErrorCaseAssert += RelError[i]; 
     2125                                        ++assertresults.nSumErrorCaseAssert; 
    19212126                                } 
    1922                         } 
    1923                         fprintf(ioQQQ,"\n"); 
     2127                                fprintf(ioASSERT,"\n"); 
     2128                        } 
     2129                        else 
     2130                                TotalInsanity(); 
    19242131                } 
    19252132 
    …… 
    19322139                        if( !dense.lgElmtOn[ipZ] ) 
    19332140                        { 
    1934                                 fprintf(ioQQQ," asserted element %ld is not turned on!\n",ipZ); 
     2141                                fprintf(ioQQQ,"PROBLEM asserted element %ld is not turned on!\n",ipZ); 
    19352142                                PredQuan[i] = 0; 
    19362143                                RelError[i] = 0.; 
    …… 
    19612168                        if( !dense.lgElmtOn[ipZ] ) 
    19622169                        { 
    1963                                 fprintf(ioQQQ," asserted element %ld is not turned on!\n",ipZ); 
     2170                                fprintf(ioQQQ,"PROBLEM asserted element %ld is not turned on!\n",ipZ); 
    19642171                                PredQuan[i] = 0.; 
    19652172                                RelError[i] = 0.; 
    …… 
    24012608                        "                  Label line    computed     asserted Rel Err Set err\n"); 
    24022609                /* now print a summary */ 
    2403                 /*lint -e771 vars conceivably not initialized */ 
    24042610                for( i=0; i<nAsserts; ++i ) 
    24052611                { 
    …… 
    25442750                        } 
    25452751                } 
    2546                 /*lint +e771 vars conceivably not initialized */ 
    25472752                fprintf( ioASSERT , " \n"); 
    25482753 
    …… 
    25592764                } 
    25602765 
     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 
    25612774                /* explain how we were compiled, but only if printing time */ 
    25622775                if( prt.lgPrintTime ) 

  • branches/newmole/source/atmdat.h

    r1942 r2023  
    205205        /** wavelengths of Hummer & Storey case B lines for H - O */ 
    206206        realnum WaveLengthCaseB[8][25][24]; 
     207 
     208        /** wavelengths for HeI case b */ 
     209        realnum *CaseBWlHeI; 
     210        /** number of HeI case b lines included */ 
     211        long int nCaseBHeI; 
    207212 
    208213        }       atmdat; 

  • branches/newmole/source/atmdat_char_tran.cpp

    r1918 r2023  
    2020#include "conv.h" 
    2121#include "atmdat.h" 
    22 /*lint -e662 possible creation out of bounds pointer */ 
    2322 
    2423/*HCTion H charge transfer ionization, H+ + A => H + A+ */ 
    …… 
    19581957        return; 
    19591958} 
    1960 /*lint +e662 possible creation out of bounds pointer */ 

  • branches/newmole/source/atmdat_HS_caseb.cpp

    r1780 r2023  
    2020/* general utility to interpolate line emissivities from the Storey & Hummer tables 
    2121 of case B emissivities.   
    22  iHi, iLo, two long integers are the principal quantum 
     22 iHi<25, iLo, the principal quantum 
    2323 numbers, and are upper and lower levels in any order.   
    24  nelem element number on physcial scale, 1 or 2 have data 
     24 nelem element number on physicial scale, 1 or 2 have data 
    2525 TempIn = temperature, and must lie within the range of the table, which depends on 
    2626 the ion charge, and is 500 - 30,000K for hydrogen.   
    …… 
    177177        yy2 = log10( atmdat.ElecTemp[iCase][nelem][ipTempHi] ); 
    178178 
    179         /*now generate the pointer to the array, expression from storey code -1 for c*/ 
     179        /*now generate the index to the array, expression from Storey code -1 for c*/ 
    180180        ip = (int)((((atmdat.ncut[iCase][nelem]-ipUp)*(atmdat.ncut[iCase][nelem]+ipUp-1))/2)+ipLo - 1); 
    181181 

  • branches/newmole/source/atmdat_lines_setup.cpp

    r1780 r2023  
    12611261        for( i=1; i<=nLevel1; ++i ) 
    12621262        { 
    1263                 if( lev2set[i] == false ) 
     1263                if( !lev2set[i] ) 
    12641264                { 
    12651265                        fprintf(ioQQQ,"PROBLEM endFindLevLine warning; line %li not claimed\n",i); 

  • branches/newmole/source/atom_feii.cpp

    r1918 r2023  
    6161#include "iso.h" 
    6262#include "pressure.h" 
    63 /*lint -e668 passing null pointer */ 
    6463 
    6564/* FeIIOvrLap handle overlapping FeII lines */ 
    …… 
    32873286#pragma OPTIMIZE ON 
    32883287#endif 
    3289 /*lint +e668 passing null pointer */ 

  • branches/newmole/source/atom_pop5.cpp

    r1780 r2023  
    224224 
    225225        /* solve matrix */ 
    226         /*lint -e740 unusual pointer case */ 
    227226        getrf_wrapper(5,5,(double*)amat,5,ipiv,&ner); 
    228227        getrs_wrapper('N',5,1,(double*)amat,5,ipiv,bvec,5,&ner); 
    229         /*lint +e740 unusual pointer case */ 
    230228 
    231229        if( ner != 0 ) 

  • branches/newmole/source/atom_seq_beryllium.cpp

    r1830 r2023  
    216216        nerror = 0; 
    217217 
    218         /*lint -e740 unusual cast */ 
    219218        getrf_wrapper(4, 4, (double*)amat, 4, ipiv, &nerror); 
    220219        getrs_wrapper('N', 4, 1, (double*)amat, 4, ipiv, bvec, 4, &nerror); 
    221