Show
Ignore:
Timestamp:
05/19/08 12:39:03 (8 months ago)
Author:
rjrw
Message:

Merged from trunk r2078:2085

Files:
1 modified

Legend:

Unmodified
Added
Removed

  • branches/newmole/source/radius_increment.cpp

    r2080 r2086  
    105105                                         * change is to only do this branch if abundance is large enough to be detected by 
    106106                                         * the ionization convergence solvers */ 
    107 #                                       define ERR_CHK  1.001 
    108                                         double err_chk = ERR_CHK; 
     107                                        const realnum err_tol = 1e-3; 
    109108                                        /* >>chng 05 sep 02, when low-T solver used solution is approximate, 
    110109                                         * and it must not matter (lot-T solver should not be used if it 
    111110                                         * does matter - so use more liberal expectation */ 
    112                                         if( iso.chTypeAtomUsed[ipISO][nelem][0]=='L' ) 
    113                                                 err_chk *= 10.; 
     111                                        //if( iso.chTypeAtomUsed[ipISO][nelem][0]=='L' ) 
     112                                        //      err_chk *= 10.; 
     113                                        realnum abund = 0.; 
     114                                        for (long n=0; n<iso.numLevels_local[ipISO][nelem]; ++n) 
     115                                        { 
     116                                                abund += StatesElem[ipISO][nelem][n].Pop; 
     117                                        } 
    114118                                        /* make sure that populations are valid, first check Pop2Ion  */ 
    115                                         if( StatesElem[ipISO][nelem][0].Pop > 
    116                                                 dense.xIonDense[nelem][nelem-ipISO]/(double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO])*err_chk && 
     119                                        if( iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15 && 
     120                                                        ! fp_equal_tol(abund*dense.xIonDense[nelem][nelem+1-ipISO],dense.xIonDense[nelem][nelem-ipISO], 
     121                                                                                                 err_tol*dense.xIonDense[nelem][nelem-ipISO]) 
    117122                                                /*>>chng 06 jul 22, only do check if pops are within tolerance of double */ 
    118                                                 iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15 ) 
     123                                                ) 
    119124                                        { 
    120125                                                fprintf(ioQQQ, 
    … …  
    122127                                                        fnzone, 
    123128                                                        ipISO , nelem , 
    124                                                         StatesElem[ipISO][nelem][0].Pop*dense.xIonDense[nelem][nelem+1-ipISO] , 
     129                                                        abund*dense.xIonDense[nelem][nelem+1-ipISO], 
    125130                                                        dense.xIonDense[nelem][nelem-ipISO] , 
    126131                                                        iso.chTypeAtomUsed[ipISO][nelem]); 
    … …  
    131136                                                        dense.xIonDense[nelem][nelem+1-ipISO]/SDIV( dense.xIonDense[nelem][nelem-ipISO] ) ); 
    132137                                                fprintf(ioQQQ, 
    133                                                         "simple %.3e Test shows Pop2Ion (%.6e) > xIonDense[nelem]/[nelem+1] (%.6e) \n",  
     138                                                        "simple %.3e Test shows Pop2Ion (%.6e) != xIonDense[nelem]/[nelem+1] (%.6e) \n",  
    134139                                                        iso.xIonSimple[ipISO][nelem], 
    135                                                         StatesElem[ipISO][nelem][0].Pop , 
     140                                                        abund , 
    136141                                                        dense.xIonDense[nelem][nelem-ipISO]/(double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO]) ); 
    137142                                        } 
    … …  
    142147                                         * if ok, rm this and all ref to lgSearch, perhaps conv.h header */ 
    143148                                        ASSERT( !conv.lgSearch ); 
    144                                         ASSERT( StatesElem[ipISO][nelem][0].Pop <= 
    145                                                 dense.xIonDense[nelem][nelem-ipISO]/ 
    146                                                 (double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO])*err_chk || 
    147                                                 /*>>chng 06 jul 22, only do check if pops are within tolerance of double */ 
    148                                                 iso.pop_ion_ov_neut[ipISO][nelem] < 1e-15); 
    149  
    150                                         /* make sure that populations are valid, first check Pop2Ion  */ 
    151                                         if( StatesElem[ipISO][nelem][0].Pop* 
    152                                                 dense.xIonDense[nelem][nelem+1-ipISO] > 
    153                                                 dense.xIonDense[nelem][nelem-ipISO]*err_chk && 
    154                                                 /*>>chng 06 jul 22, only do check if pops are within tolerance of double */ 
    155                                                 iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15) 
    156                                         { 
    157                                                 fprintf(ioQQQ," DISASTER PopLo fnz %.2f ipISO %li nelem %li pop_ion_ov_neut %.2e 2pops %.6e %.6e solver:%s\n",  
    158                                                         fnzone, 
    159                                                         ipISO , nelem , 
    160                                                         iso.pop_ion_ov_neut[ipISO][nelem] ,  
    161                                                         StatesElem[ipISO][nelem][0].Pop*dense.xIonDense[nelem][nelem+1-ipISO], 
    162                                                         dense.xIonDense[nelem][nelem-ipISO]*err_chk , 
    163                                                         iso.chTypeAtomUsed[ipISO][nelem]); 
    164                                         } 
    165                                         ASSERT(  
    166                                                 /*(conv.lgSearch || !conv.nPres2Ioniz) || */ 
    167                                                 (StatesElem[ipISO][nelem][0].Pop*dense.xIonDense[nelem][nelem+1-ipISO]<= 
    168                                                         dense.xIonDense[nelem][nelem-ipISO]*err_chk) || 
    169                                                 /*>>chng 06 jul 22, only do check if pops are within tolerance of double */ 
    170                                                 iso.pop_ion_ov_neut[ipISO][nelem] < 1e-15); 
    171 #                                       undef ERR_CHK 
     149                                        ASSERT( /*>>chng 06 jul 22, only do check if pops are within tolerance of double */ 
     150                                                iso.pop_ion_ov_neut[ipISO][nelem] < 1e-15 || 
     151                                                fp_equal_tol(abund*dense.xIonDense[nelem][nelem+1-ipISO],dense.xIonDense[nelem][nelem-ipISO], 
     152                                                                                                 err_tol*dense.xIonDense[nelem][nelem-ipISO]) ); 
    172153                                } 
    173154                        }