Show
Ignore:
Timestamp:
05/10/08 09:03:02 (2 months ago)
Author:
peter
Message:

Merging all changes from mainline upto r2033, except r1902.

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • branches/c08_branch/source/iso_level.cpp

    r1916 r2034  
    2424#include "trace.h" 
    2525 
    26 STATIC void PrtHydroTrace2( long ipISO, long int nelem, const multi_arr<double,2,C_TYPE>& SaveZ ); 
    27  
    28 /*lint -e662 Possible access of out-of-bounds pointer*/ 
    29 /*===================================================================================*/ 
    3026/* solve for level populations  */ 
    3127void iso_level( const long int ipISO, const long int nelem) 
    3228{ 
    33         long int n, 
    34                 ipHi, 
     29        long int ipHi, 
    3530                ipLo, 
    3631                i, 
    37                 j, 
    3832                level, 
    3933                level_error; 
    …… 
    4943                TooSmall; 
    5044        bool lgNegPop=false; 
    51         static double SaveHe23S_photorate=0.; 
    5245        valarray<int32> ipiv(numlevels_local);  
    5346        /* this block of variables will be obtained and freed here */ 
    …… 
    7164        /* now do the 2D array */ 
    7265        z.alloc(numlevels_local,numlevels_local); 
    73  
    74         /* >>chng 03 may 01, move defn of these two up here, before entering matrix */ 
    75         /* >>chng 01 jun 21, add */ 
    76         /* charge transfer with helium itself, onto everything heavier */ 
    77         if( nelem == ipHELIUM && ipISO == ipHE_LIKE ) 
    78         { 
    79                 /* >>chng 04 mar 09, add this average of old and new */ 
    80                 /* take average of old and new 23S photo rate - can be dominated by Lya and a bit noisy  
    81                  * until ots rates have settled down */ 
    82                 if( nzone>1 ) 
    83                 { 
    84                         double frac=0.5; 
    85                         iso.gamnc[ipISO][nelem][ipHe2s3S] = frac*iso.gamnc[ipISO][nelem][ipHe2s3S] + 
    86                                 (1.-frac)*SaveHe23S_photorate; 
    87                 } 
    88                 SaveHe23S_photorate = iso.gamnc[ipISO][nelem][ipHe2s3S]; 
    89         } 
    9066 
    9167        /* >>chng 05 aug 17, must use real electron density for collisional ionization of H 
    …… 
    215191                iso.pop_ion_ov_neut[ipISO][nelem] = iso.xIonSimple[ipISO][nelem]; 
    216192                StatesElem[ipISO][nelem][0].Pop =  1./SDIV(iso.pop_ion_ov_neut[ipISO][nelem]); 
    217                 for( n=1; n < numlevels_local; n++ ) 
     193                for( long n=1; n < numlevels_local; n++ ) 
    218194                { 
    219195                        StatesElem[ipISO][nelem][n].Pop =  0.; 
    …… 
    223199        else 
    224200        { 
     201                ASSERT( NISO == 2 ); 
    225202                long SpinUsed[NISO] = {2,1}; 
    226203                long indexNmaxP = 
    …… 
    284261 
    285262                                z[level][level] += RadDecay + pump_down + coll_down; 
    286                                 z[level][ipLo] = - (RadDecay + pump_down + coll_down);                                  
     263                                z[level][ipLo] = - (RadDecay + pump_down + coll_down); 
    287264 
    288265                                if( level == indexNmaxP ) 
    …… 
    300277                                        iso.qTot2S[ipISO][nelem] += coll_down/dense.EdenHCorr; 
    301278                                } 
    302                                 if( (ipLo == 1) && (ipISO==ipH_LIKE || (StElm[level].S==1)) ) 
     279                                if( (ipLo == 1) && (ipISO==ipH_LIKE || StElm[level].S==1) ) 
    303280                                { 
    304281                                        iso.qTot2S[ipISO][nelem] += coll_up/dense.EdenHCorr; 
    …… 
    411388                        Save_creation[ipLo] = creation[ipLo]; 
    412389 
    413                 if( trace.lgTrace &&  
    414                         trace.lgIsoTraceFull[ipISO] && (nelem == trace.ipIsoTrace[ipISO]) ) 
    415                 { 
    416                         /** \todo 2 generalize this, can probably then remove the trace below. */ 
    417                         if( ipISO == ipH_LIKE ) 
    418                                 PrtHydroTrace2(ipH_LIKE, nelem, SaveZ); 
    419                         else if( ipISO == ipHE_LIKE ) 
    420                         { 
    421                                 fprintf( ioQQQ, "  pop level     others => (iso_level)\n" ); 
    422                                 for( n=0; n < numlevels_local; n++ ) 
     390                if( trace.lgTrace && trace.lgIsoTraceFull[ipISO] && (nelem == trace.ipIsoTrace[ipISO]) ) 
     391                { 
     392                        fprintf( ioQQQ, "  pop level     others => (iso_level)\n" ); 
     393                        for( long n=0; n < numlevels_local; n++ ) 
     394                        { 
     395                                fprintf( ioQQQ, "  %s %s %2ld", iso.chISO[ipISO], elementnames.chElementNameShort[nelem], n ); 
     396                                for( long j=0; j < numlevels_local; j++ ) 
    423397                                { 
    424                                         fprintf( ioQQQ, "       HeI%2ld", n ); 
    425                                         for( j=0; j < numlevels_local; j++ ) 
    426                                         { 
    427                                                 fprintf( ioQQQ,"\t%.9e", z[j][n] ); 
    428                                         } 
    429                                         fprintf( ioQQQ, "\n" ); 
    430                                 } 
    431                                 fprintf(ioQQQ," recomb          "); 
    432                                 for( n=0; n < numlevels_local; n++ ) 
    433                                 { 
    434                                         fprintf( ioQQQ,"\t%.9e", creation[n] ); 
     398                                        fprintf( ioQQQ,"\t%.9e", z[j][n] ); 
    435399                                } 
    436400                                fprintf( ioQQQ, "\n" ); 
    437                                 fprintf(ioQQQ," recomb ct %.2e co %.2e hectr %.2e hctr %.2e\n", 
    438                                         atmdat.HeCharExcRecTotal, 
    439                                         findspecies("CO")->hevmol , 
    440                                         atmdat.HeCharExcRecTo[nelem][nelem-1]*dense.xIonDense[ipHELIUM][0] , 
    441                                         atmdat.HCharExcRecTo[nelem][nelem-1]*dense.xIonDense[ipHYDROGEN][0] ); 
    442                         } 
     401                        } 
     402                        fprintf(ioQQQ," recomb          "); 
     403                        for( long n=0; n < numlevels_local; n++ ) 
     404                        { 
     405                                fprintf( ioQQQ,"\t%.9e", creation[n] ); 
     406                        } 
     407                        fprintf( ioQQQ, "\n" ); 
     408                        fprintf(ioQQQ," recomb ct %.2e co %.2e hectr %.2e hctr %.2e\n", 
     409                                atmdat.HeCharExcRecTotal, 
     410                                findspecies("CO")->hevmol , 
     411                                atmdat.HeCharExcRecTo[nelem][nelem-ipISO]*dense.xIonDense[ipHELIUM][0] , 
     412                                atmdat.HCharExcRecTo[nelem][nelem-ipISO]*dense.xIonDense[ipHYDROGEN][0] ); 
    443413                } 
    444414 
    …… 
    463433                        double BigSoln= 0.; 
    464434                        error[level] = 0.; 
    465                         for( n=ipH1s; n < numlevels_local; n++ ) 
     435                        for( long n=ipH1s; n < numlevels_local; n++ ) 
    466436                        { 
    467437                                /* remember the largest value of the soln matrix to div by below */ 
    …… 
    534504                                        level, 
    535505                                        StatesElem[ipISO][nelem][level].Pop ); 
    536  
    537506                        } 
    538507 
    …… 
    643612        } 
    644613 
    645         if( ipISO == ipHE_LIKE ) 
    646         { 
    647  
    648                 if( nzone>0 && nelem==ipHELIUM ) 
    649                 { 
    650                         /* find fraction of He0 destructions due to photoionization of 2 ^S */ 
    651                         double ratio; 
    652                         double rateOutOf2TripS = StatesElem[ipHE_LIKE][nelem][ipHe2s3S].Pop * iso.RateLevel2Cont[ipHE_LIKE][nelem][ipHe2s3S]; 
     614        if( ipISO == ipHE_LIKE && nelem==ipHELIUM && nzone>0 ) 
     615        { 
     616                /* find fraction of He0 destructions due to photoionization of 2^3S */ 
     617                double ratio; 
     618                double rateOutOf2TripS = StatesElem[ipISO][nelem][ipHe2s3S].Pop * iso.RateLevel2Cont[ipISO][nelem][ipHe2s3S]; 
     619                if( rateOutOf2TripS > SMALLFLOAT ) 
     620                { 
     621                        ratio = rateOutOf2TripS / 
     622                                ( rateOutOf2TripS + StatesElem[ipISO][nelem][ipHe1s1S].Pop*ionbal.RateIonizTot[nelem][nelem-ipISO] ); 
     623                } 
     624                else 
     625                { 
     626                        ratio = 0.; 
     627                } 
     628                if( ratio > he.frac_he0dest_23S ) 
     629                { 
     630                        /* remember zone where this happended and fraction, and frac due to photoionization */ 
     631                        he.nzone = nzone; 
     632                        he.frac_he0dest_23S = ratio; 
     633                        rateOutOf2TripS = StatesElem[ipISO][nelem][ipHe2s3S].Pop *iso.gamnc[ipISO][nelem][ipHe2s3S]; 
    653634                        if( rateOutOf2TripS > SMALLFLOAT ) 
    654635                        { 
    655                                 ratio = rateOutOf2TripS / 
    656                                         ( rateOutOf2TripS + StatesElem[ipHE_LIKE][nelem][ipHe1s1S].Pop*ionbal.RateIonizTot[nelem][nelem-ipISO] ); 
     636                                he.frac_he0dest_23S_photo = rateOutOf2TripS / 
     637                                        ( rateOutOf2TripS + StatesElem[ipISO][nelem][ipHe1s1S].Pop*ionbal.RateIonizTot[nelem][nelem-ipISO] ); 
    657638                        } 
    658639                        else 
    659640                        { 
    660                                 ratio = 0.; 
    661                         } 
    662                         if( ratio > he.frac_he0dest_23S ) 
    663                         { 
    664                                 /* remember zone where this happended and fraction, and frac due to photoionization */ 
    665                                 he.nzone = nzone; 
    666                                 he.frac_he0dest_23S = ratio; 
    667                                 rateOutOf2TripS = StatesElem[ipHE_LIKE][nelem][ipHe2s3S].Pop *iso.gamnc[ipHE_LIKE][nelem][ipHe2s3S]; 
    668                                 if( rateOutOf2TripS > SMALLFLOAT ) 
    669                                 { 
    670                                         he.frac_he0dest_23S_photo = rateOutOf2TripS / 
    671                                                 ( rateOutOf2TripS + StatesElem[ipHE_LIKE][nelem][ipHe1s1S].Pop*ionbal.RateIonizTot[nelem][nelem-ipISO] ); 
    672                                 } 
    673                                 else 
    674                                 { 
    675                                         he.frac_he0dest_23S_photo = 0.; 
    676                                 } 
     641                                he.frac_he0dest_23S_photo = 0.; 
    677642                        } 
    678643                } 
    …… 
    699664        return; 
    700665} 
    701  
    702 STATIC void PrtHydroTrace2( long int ipISO, long int nelem, const multi_arr<double,2,C_TYPE>& SaveZ ) 
    703 { 
    704         long ipHi , j , n; 
    705         double collider; 
    706  
    707         DEBUG_ENTRY( "PrtHydroTrace2()" ); 
    708  
    709         /* >>chng 05 aug 17, use eden as collider for H and corrected eden for heavier 
    710          * nuclei - in hot PDR H is collisionally ionized, but not by other H0 since 
    711          * collision is homonuclear */ 
    712         if( nelem==ipHYDROGEN ) 
    713         { 
    714                 /* special version for H0 onto H0 */ 
    715                 collider = dense.EdenHontoHCorr; 
    716         } 
    717         else 
    718         { 
    719                 collider = dense.EdenHCorr; 
    720         } 
    721  
    722         if( nelem == ipHYDROGEN ) 
    723         { 
    724                 double sum=0 , HRateDestGnd[10]; 
    725                 /* >>chng 06 aug 29, all of these from numLevels_max to _local. */ 
    726                 for(ipHi=1; ipHi<iso.numLevels_local[ipISO][nelem]; ++ipHi ) 
    727                 { 
    728                         if( Transitions[ipISO][nelem][ipHi][ipH1s].Emis->Aul <= iso.SmallA ) 
    729                                 continue; 
    730  
    731                         sum += Transitions[ipISO][nelem][ipHi][ipH1s].Emis->pump; 
    732                 } 
    733                 fprintf(ioQQQ,"PrtHydroTrace2 pump\t%.2f", fnzone); 
    734                 for(ipHi=1; ipHi<iso.numLevels_local[ipISO][nelem]; ++ipHi ) 
    735                 { 
    736                         if( Transitions[ipISO][nelem][ipHi][ipH1s].Emis->Aul <= iso.SmallA ) 
    737                                 continue; 
    738  
    739                         fprintf(ioQQQ,"\t%.3f",Transitions[ipISO][nelem][ipHi][ipH1s].Emis->pump/SDIV(sum)); 
    740                 } 
    741                 fprintf(ioQQQ,"\n"); 
    742  
    743                 /* remember ground state destruction rate */ 
    744                 HRateDestGnd[0] = iso.RateLevel2Cont[ipISO][nelem][ipH1s]; 
    745                 HRateDestGnd[1] = iso.gamnc[ipISO][nelem][ipH1s]/ 
    746                         iso.RateLevel2Cont[ipISO][nelem][ipH1s]; 
    747                 /* >>chng 05 aug 17, from EdeHCorr to eden as per above comment */ 
    748                 HRateDestGnd[2] = iso.ColIoniz[ipISO][nelem][ipH1s]* 
    749                         collider/iso.RateLevel2Cont[ipISO][nelem][ipH1s]; 
    750                 HRateDestGnd[3] = secondaries.csupra[nelem][nelem]/iso.RateLevel2Cont[ipISO][nelem][ipH1s]; 
    751                 HRateDestGnd[4] = secondaries.Hx12[ipISO][nelem][ipH2p]* 
    752                         9./iso.RateLevel2Cont[ipISO][nelem][ipH1s]; 
    753                 HRateDestGnd[5] = atmdat.HCharExcIonTotal/ 
    754                         iso.RateLevel2Cont[ipISO][nelem][ipH1s]; 
    755                 HRateDestGnd[6] = atmdat.HCharExcIonTotal; 
    756                 HRateDestGnd[7] = sum/iso.RateLevel2Cont[ipISO][nelem][ipH1s]; 
    757                 HRateDestGnd[8] = sum; 
    758  
    759                 fprintf(ioQQQ," grnd dest fracs\t%.2f",fnzone); 
    760                 for( j=ipH1s; j < 9; j++ ) 
    761                 { 
    762                         fprintf( ioQQQ,"\t"); 
    763                         fprintf( ioQQQ,PrintEfmt("%8.1e", HRateDestGnd[j] )); 
    764                 } 
    765                 fprintf(ioQQQ,"\thcoldc\t%e\n", dense.EdenHCorr ); 
    766                 fprintf( ioQQQ, "\n" ); 
    767         } 
    768  
    769         fprintf( ioQQQ, "  pop level     others => (iso_level)\n" ); 
    770         /* >>chng 06 aug 29, from numLevels_max to _local */ 
    771         for( n=ipH1s; n < iso.numLevels_local[ipISO][nelem]; n++ ) 
    772         { 
    773                 fprintf( ioQQQ, "       HII%2ld", n ); 
    774                 for( j=ipH1s; j < iso.numLevels_local[ipISO][nelem]; j++ ) 
    775                 { 
    776                         fprintf( ioQQQ," "); 
    777                         /*fprintf( ioQQQ,PrintEfmt("%8.1e", SaveZ[j][n] ) );*/ 
    778                         fprintf( ioQQQ,"%.9e\t", SaveZ[j][n] ); 
    779                 } 
    780                 fprintf( ioQQQ, "\n" ); 
    781         } 
    782  
    783         return; 
    784 } 
    785 /*lint +e662 Possible access of out-of-bounds pointer*/