Changeset 1942 for branches/newmole/source/assert_results.cpp

Timestamp:

04/12/08 09:30:35 (9 months ago)

Author:

rjrw

Message:

Merged from trunk r1917:1941

Files:

• branches/newmole/source/assert_results.cpp (modified) (8 diffs)

branches/newmole/source/assert_results.cpp

@@ -10 +10 @@
 #include "iso.h"
 #include "called.h"
+#include "atmdat.h"
 #include "hcmap.h"
 #include "thermal.h"
@@ -47 +48 @@
 static char **chAssertLineLabel;
 
+/* element to assert H-like case B */
+static char chElemLabelCaseB[5];
+static long int nelemCaseB , nISOCaseB;
+
 /* this will be = for equal, < or > for limit */
 static char *chAssertLimit;
@@ -344 +349 @@
 
         /* assert line "LINE"  --  key is ine_ since linear option appears on some commands */
-        else if( nMatch("INE ",input.chCARDCAPS ) )
+        else if( nMatch(" LINE ",input.chCARDCAPS ) )
         {
                 if(  nMatch("LUMI",input.chCARDCAPS ) || nMatch("INTE",input.chCARDCAPS ))
@@ -440 +445 @@
 
         /* assert departure coefficients */
+        else if( nMatch("CASE",input.chCARDCAPS ) )
+        {
+                /* this is Case B for some element */
+                strcpy( chAssertType[nAsserts] , "CB" );
+                i = 5;
+                /* this is relative error */
+                AssertError[nAsserts] =
+                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
+                if( lgEOL )
+                        /* default error was set in define above */
+                        AssertError[nAsserts] = DEF_ERROR;
+                AssertQuantity[nAsserts] = 0;
+                strcpy( chAssertLineLabel[nAsserts] , "Ca B" ); 
+                wavelength[nAsserts] = 0.;
+
+                /* assert case b for H or He */
+                if( nMatch("H-LI",input.chCARDCAPS ) )
+                {
+                        /* H-like - now get an element */
+                        if( (nelem = GetElem( input.chCARDCAPS )) < 0 )
+                        {
+                                /* no name found */
+                                fprintf(ioQQQ, "assert H-like case B did not find an element on this line, sorry %s\n",
+                                        input.chCARDCAPS );
+                                cdEXIT(EXIT_FAILURE);
+                        }
+                        if( nelem>7 )
+                        {
+                                /* beyond reach of tables */
+                                fprintf(ioQQQ, "assert H-like cannot do elements heavier than O, sorry %s\n",
+                                        input.chCARDCAPS );
+                                cdEXIT(EXIT_FAILURE);
+                        }
+                        nelemCaseB = nelem;
+                        nISOCaseB = 0;
+                        strcpy( chElemLabelCaseB , elementnames.chElementSym[nelem] );
+                        char chNumb[4];
+                        sprintf( chNumb , "%2li" , nelem+1 );
+                        strcat( chElemLabelCaseB , chNumb );
+                }
+                else if( nMatch("HE-L",input.chCARDCAPS ) )
+                {
+                        /* He-like - only helium itself */
+                        nelemCaseB = ipHELIUM;
+                        nISOCaseB = 0;
+                        strcpy( chElemLabelCaseB , "He 1" );
+                        /* beyond reach of tables */
+                        fprintf(ioQQQ, "assert He-like case b not implemented yet %s\n",
+                                input.chCARDCAPS );
+                        cdEXIT(EXIT_FAILURE);
+                }
+        }
         else if( nMatch("DEPA",input.chCARDCAPS ) )
         {
@@ -447 +504 @@
                         FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
                 if( lgEOL )
-                {
                         NoNumb(input.chOrgCard);
-                }
 
                 /* this is relative error, max departure from unity of any level or std */
@@ -455 +510 @@
                         FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
                 if( lgEOL )
-                {
                         /* default error was set in define above */
                         AssertError[nAsserts] = DEF_ERROR;
-                }
 
                 /* H-like key means do one of the hydrogenic ions */
@@ -1828 +1881 @@
                         /* uncertainty in difference is zero */
                         AssertError[i] = 0.;
+                }
+
+                else if( strcmp( chAssertType[i] , "CB") == 0 )
+                {
+                        int iCase = 1;
+                        RelError[i] = 0.;
+                        fprintf(ioQQQ,"        Species  nHi nLo  Wl  Computed  Asserted       error\n");
+                        long nHighestPrinted = StatesElem[nISOCaseB][nelemCaseB][iso.numPrintLevels[nISOCaseB][nelemCaseB]-1].n;
+                        for( long int ipLo=1+iCase; ipLo< MIN2(24,nHighestPrinted-1); ++ipLo )
+                        {
+                                for( long int ipHi=ipLo+1; ipHi< MIN2(25,nHighestPrinted); ++ipHi )
+                                {
+                                        /* assert the line */
+                                        realnum wl = atmdat.WaveLengthCaseB[nelemCaseB][ipHi][ipLo];
+                                        double relint , absint,CBrelint , CBabsint;
+                                        cdLine( "Ca B" , wl , &CBrelint , &CBabsint );
+                                        CBabsint = pow( 10., CBabsint );
+                                        double error;
+                                        if( cdLine( chElemLabelCaseB , wl , &relint , &absint ) >0)
+                                        {
+                                                absint = pow( 10., absint );
+                                                error = (CBabsint - absint)/MAX2(CBabsint , absint);
+                                                fprintf(ioQQQ,"ChkAssert %s %3li %3li ", 
+                                                        chElemLabelCaseB , ipHi , ipLo );
+                                                prt_wl(ioQQQ, wl );
+                                                fprintf(ioQQQ," %.2e %.2e %10.3f", 
+                                                        absint , CBabsint , error );
+                                        }
+                                        else
+                                                TotalInsanity();
+                                        if( fabs(error) > AssertError[i]  )
+                                                fprintf(ioQQQ , " BOTCH \n");
+                                        else
+                                                fprintf(ioQQQ , "\n");
+
+                                        PredQuan[i] = 0;
+                                        AssertQuantity[i] = 0.;
+                                        RelError[i] = MAX2( RelError[i] , fabs(error) );
+                                }
+                        }
+                        fprintf(ioQQQ,"\n");
                 }
 
@@ -2392 +2486 @@
 
                                 {
-                                        /*@-redef@*/
                                         enum {DEBUG_LOC=false};
-                                        /*@+redef@*/
 
                                         if( DEBUG_LOC )