Changeset 2035

Show
Ignore:
Timestamp:
05/10/08 10:47:02 (2 months ago)
Author:
rporter
Message:

Make a few minor changes to iso error generation. This includes a change to MyGaussRand? such the passed parameter is now one standard deviation and the returned value is within three standard deviations. Fix and activate the recent issue in sanity_check.cpp (was broken in c07 and is commented out in release branch).

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • trunk/source/iso_collide.cpp

    r1969 r2035  
    3232                t_ADfA::Inst().coll_ion( nelem+1, 1+ipISO, phycon.te ); 
    3333 
     34        iso_put_error(ipISO,nelem,iso.numLevels_max[ipISO][nelem],0,IPCOLLIS,0.20f); 
     35         
    3436        for( long ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ ) 
    3537        { 
    …… 
    5456                if( ipHi < iso.numLevels_max[ipISO][nelem] - 1 ) 
    5557                        iso.ColIoniz[ipISO][nelem][ipHi] *= iso.lgColl_ionize[ipISO]; 
     58         
     59                iso_put_error(ipISO,nelem,iso.numLevels_max[ipISO][nelem],ipHi,IPCOLLIS,0.20f); 
    5660        } 
    5761         

  • trunk/source/iso_error.cpp

    r1732 r2035  
    6969                iso_put_error(ipISO,nelem,iso.numLevels_max[ipISO][nelem],ipHi,IPRAD,ErrorToPut); 
    7070        } 
     71 
    7172        /* Error for total recombination */ 
    7273        iso_put_error(ipISO,nelem,iso.numLevels_max[ipISO][nelem],iso.numLevels_max[ipISO][nelem],IPRAD,0.02f); 
     74        iso.ErrorFactor[ipISO][nelem][iso.numLevels_max[ipISO][nelem]][iso.numLevels_max[ipISO][nelem]][IPRAD] = 
     75                (realnum)MyGaussRand( iso.Error[ipISO][nelem][iso.numLevels_max[ipISO][nelem]][iso.numLevels_max[ipISO][nelem]][IPRAD] ); 
    7376 
    7477        for( ipHi=1; ipHi<= iso.numLevels_max[ipISO][nelem]; ipHi++ ) 
    7578        { 
    7679                /* >>chng 06 mar 15, the upper limit incorrectly went to numLevels_max */ 
    77                 for( ipLo=0; ipLo<= ipHi; ipLo++ ) 
     80                for( ipLo=0; ipLo< ipHi; ipLo++ ) 
    7881                { 
    7982                        for( typeOfRate=0; typeOfRate<=1; typeOfRate++ ) 

  • trunk/source/iso_radiative_recomb.cpp

    r1969 r2035  
    173173                        else 
    174174                        { 
    175                                 ASSERT( (TotRRThisN<TotRRLastN) || (ThisN<=2L) || (phycon.te>3E7) || (nelem!=ipHELIUM) ); 
     175                                ASSERT( iso.lgRandErrGen[ipISO] || (TotRRThisN<TotRRLastN) || (ThisN<=2L) || (phycon.te>3E7) || (nelem!=ipHELIUM) ); 
    176176                                LastN = ThisN; 
    177177                                ThisN = N_(ipLevel); 

  • trunk/source/sanity_check.cpp

    r1960 r2035  
    507507        } 
    508508 
    509         fixit();/* make this work */ 
    510 #if     0 
    511509        /* check photo cross sections for H */ 
    512510        long ipISO = ipH_LIKE; 
    513511        nelem = 0; 
    514         /* loop starts at 3, the first level with n = n and full l */ 
    515         for( n=3; n < MIN2(100, iso.n_HighestResolved_max[ipISO][nelem]); ++n ) 
    516         { 
    517                 realnum anu[1]={1.} , cs[1]; 
    518                 double energy; 
    519                 long index = iso.QuantumNumbers2Index[ipISO][nelem][n][0][2]; 
    520  
    521                 /* photon energy where cross section will be evaluated, 
    522                  * this is in Ryd */ 
    523                 energy = iso.xIsoLevNIonRyd[ipISO][nelem][index]; 
    524                 anu[0] = (realnum)(energy*1.01); 
    525  
    526                 H_photo_cs( photon , N_(n), L_(n), nelem+1 ); 
    527                 hypho( 
    528                         /* Z-Nelec, the residual charge, 1 for hydrogen, 2 for helium */ 
    529                         (double)(nelem+1),  
    530                         /* principal quantum number */ 
    531                   n,  
    532                   /* lowest angular momentum */ 
    533                   0,  
    534                   /* highest angular momentum */ 
    535                   n-1,  
    536                   /* scaled lowest photon energy,  
    537                    * => incorrect?? in units of zed^2/n^2,  
    538                    * at which the cs will be done */ 
    539                   energy,  
    540                   /* number of points to do */ 
    541                   1,  
    542                   /* array of energies (in units given above) */ 
    543                   anu ,  
    544                   /* calculated photoionization cross section in cm^-2 */ 
    545                   cs); 
    546  
    547                 error = fabs(cs[0] - opac.OpacStack[iso.ipOpac[ipISO][nelem][index]-1] )/ 
    548                         (        (cs[0] + opac.OpacStack[iso.ipOpac[ipISO][nelem][index]-1] )/2.); 
    549                 /*fprintf(ioQQQ,"z=%ld n=%ld error %g old %e new %e\n",nelem,n, error, 
    550                         opac.OpacStack[iso.ipOpac[ipISO][nelem][n]-1] ,cs[0] );*/ 
    551                 if( error > 0.05 ) 
    552                 { 
    553                         fprintf(ioQQQ , "SanityCheck finds insane H photo cs\n"); 
    554                         lgOK = false; 
    555                 } 
    556         } 
    557 #endif 
     512        for( long n=1; n <= iso.n_HighestResolved_max[ipISO][nelem]; ++n ) 
     513        { 
     514                double rel_photon_energy = 1. + FLT_EPSILON*2.; 
     515                for( long l=0; l < n; l++ ) 
     516                { 
     517                        double cs; 
     518                        long index = iso.QuantumNumbers2Index[ipISO][nelem][n][l][2]; 
     519 
     520                        cs = H_photo_cs( rel_photon_energy, n, l, nelem+1 ); 
     521 
     522                        error = fabs(cs - opac.OpacStack[iso.ipOpac[ipISO][nelem][index]-1] )/ 
     523                                ( (cs + opac.OpacStack[iso.ipOpac[ipISO][nelem][index]-1] )/2.); 
     524 
     525                        if( error > 0.05 ) 
     526                        { 
     527                                fprintf(ioQQQ , "SanityCheck finds insane H photo cs, n,l = %li, %li, expected, found = %e, %e\n"); 
     528                                lgOK = false; 
     529                        } 
     530                } 
     531        } 
    558532 
    559533        /********************************************************* 

  • trunk/source/service.cpp

    r1954 r2035  
    17981798 
    17991799/* MyGaussRand takes as input a percent uncertainty less than 50% 
    1800  * (expressed as 0.5). The routine then assumes this input variable represents two 
    1801  * standard deviations about a mean of unity, and returns a random number within 
    1802  * that range.  A hard cutoff is imposed at the two standard deviations, which  
    1803  * eliminates roughly 2% of the normal distribution.  In other words, the routine 
     1800 * (expressed as 0.5). The routine then assumes this input variable represents one 
     1801 * standard deviation about a mean of unity, and returns a random number within 
     1802 * that range.  A hard cutoff is imposed at two standard deviations, which  
     1803 * eliminates roughly 5% of the normal distribution.  In other words, the routine 
    18041804 * returns a number in a normal distribution with standard deviation equal to 
    1805  * half of the input, and returns a number between 1-2*stdev and 1+2*stdev. */ 
     1805 * the input.  The number will be between 1-3*stdev and 1+3*stdev. */ 
    18061806double MyGaussRand( double PctUncertainty ) 
    18071807{ 
    …… 
    18121812 
    18131813        ASSERT( PctUncertainty < 0.5 ); 
    1814         /* We want this "percent uncertainty to represent two standard deviations */ 
    1815         StdDev = 0.5 * PctUncertainty; 
     1814        /* We want this "percent uncertainty" to represent one standard deviation */ 
     1815        StdDev = PctUncertainty; 
    18161816 
    18171817        do 
    …… 
    18201820                result = 1.+RandGauss( 0., StdDev ); 
    18211821        } 
    1822         /* This will give us a result that is less than or equal to the 
    1823          * percent uncertainty about 98% of the time. */ 
    1824         while( (result < 1.-PctUncertainty) || (result > 1+PctUncertainty) ); 
     1822        /* only allow values that are within 3 standard deviations */ 
     1823        while( (result < 1.-3.*PctUncertainty) || (result > 1.+3.*PctUncertainty) ); 
    18251824 
    18261825        ASSERT( result>0. && result<2. );