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/mole_h2_coll.cpp

    r1815 r2034  
    2121static bool **lgDefn_H2He_coll; 
    2222 
    23 STATIC realnum H2_CollidRateEvalOne( long iVibHi, long iRotHi,long iVibLo, 
    24                                                          long iRotLo, long ipHi , long ipLo , long  nColl ) 
     23/* compute rate coefficient for a single quenching collision */ 
     24STATIC realnum H2_CollidRateEvalOne(  
     25        /*returns collision rate coefficient, cm-3 s-1 for quenching collision 
     26         * from this upper state */ 
     27         long iVibHi, long iRotHi,long iVibLo, 
     28         /* to this lower state */ 
     29        long iRotLo, long ipHi , long ipLo ,  
     30        /* colliders are H, He, H2(ortho), H2(para), and H+ */ 
     31        long  nColl ) 
    2532{ 
    2633        double fitted; 
    2734        realnum rate; 
    28         double t = phycon.te/1000. + 1.; 
    29         double t2 = POW2(t); 
     35        double t3Plus1 = phycon.te/1000. + 1.; 
     36        double t3Plus1Squared = POW2(t3Plus1); 
    3037        /* these are fits to the existing collision data  
    3138         * used to create g-bar rates */ 
    …… 
    6875                } 
    6976 
    70                 /* >>chng 07 apr 2, GS 
    71                 * this is option to use g-bar g bar guess of rate coefficient for  
    72                 * collision with He, this does not change ortho & para  
    73                 * code had always applied g-bar, which was incorrect 
    74                 * turn mole.lgColl_gbar on/off with atom h2 gbar on off */ 
     77                /* use g-bar g bar guess of rate coefficient for  
     78                 * collision with He, this does not change ortho & para  
     79                 * turn mole.lgColl_gbar on/off with atom h2 gbar on off */ 
    7580                else if( mole.lgColl_gbar  &&  
    7681                        (H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVibLo][iRotLo]==0) ) 
    …… 
    139144        else if( CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][0]!= 0 ) 
    140145        { 
    141                 /* evaluate collision rates for those with real collision data */ 
     146                /* these are the fits from  
     147                 *>>refer       H2      coll    Le Bourlot, J., Pineau des Forets,  
     148                 *>>refercon    G., & Flower, D.R. 1999, MNRAS, 305, 802 
     149                 * evaluate collision rates for those with real collision data */ 
    142150 
    143151                double r = CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][0] +  
    144                         CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][1]/t +  
    145                         CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][2]/t2; 
     152                        CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][1]/t3Plus1 +  
     153                        CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][2]/t3Plus1Squared; 
    146154 
    147155                rate = (realnum)pow(10.,r)*mole.lgColl_deexec_Calc; 
    …… 
    201209} 
    202210 
    203 /*H2_CollidRateEvalAll - set H2 collision rates */ 
     211/*H2_CollidRateEvalAll - set H2 collision rate coefficients */ 
    204212void H2_CollidRateEvalAll( void ) 
    205213{ 
    …… 
    302310 
    303311        /* >>chng 02 nov 13, the Sun and Dalgarno rates diverge to + inf below this temp */ 
    304         /* >>chng 05 feb 09, do not return zero when t < 100 - instead, don't let T fall below 100 */ 
     312        /* >>chng 05 feb 09, do not return zero when T < 100 - instead, don't let T fall below 100 */ 
    305313        double excit1; 
    306314        double te_used = MAX2( 100. , phycon.te ); 
    …… 
    657665{ 
    658666 
    659         double k, t, b2, c2, f2, t2; 
     667        double k, t3Plus1, b2, c2, f2, t3; 
    660668 
    661669        DEBUG_ENTRY( "H2_He_coll()" ); 
    …… 
    695703                temp = MIN2( 1e4 , temp ); 
    696704 
    697                 t2 = temp/1e3;  
    698                 /* t = T*10^-3 + 1*/ 
    699                 t = t2+1;  
    700                 b2 = H2_He_coll_fit_par[init][final][1]/(H2_He_coll_fit_par[init][final][3]*t2+1); 
    701                 c2 = H2_He_coll_fit_par[init][final][2]/(t*t); 
     705                t3 = temp/1e3;  
     706                /* t3Plus1 = T/1000 + 1*/ 
     707                t3Plus1 = t3+1;  
     708                b2 = H2_He_coll_fit_par[init][final][1]/(H2_He_coll_fit_par[init][final][3]*t3+1); 
     709                c2 = H2_He_coll_fit_par[init][final][2]/(t3Plus1*t3Plus1); 
    702710                { 
    703711                        enum {DEBUG_LOC=false}; 
    …… 
    707715                                        init,final, 
    708716                                        H2_He_coll_fit_par[init][final][5], 
    709                                         H2_He_coll_fit_par[init][final][6]*t2+1,  
     717                                        H2_He_coll_fit_par[init][final][6]*t3+1,  
    710718                                        H2_He_coll_fit_par[init][final][7] 
    711                                 /*pow(H2_He_coll_fit_par[init][final][6]*t2+1,H2_He_coll_fit_par[init][final][7])*/ 
     719                                /*pow(H2_He_coll_fit_par[init][final][6]*t3+1,H2_He_coll_fit_par[init][final][7])*/ 
    712720                                ); 
    713721                        } 
    714722                } 
    715723                /* this is log of f2 - see whether it is within bounds */ 
    716                 sum1 = H2_He_coll_fit_par[init][final][7] * log10( H2_He_coll_fit_par[init][final][6]*t2+1. ); 
     724                sum1 = H2_He_coll_fit_par[init][final][7] * log10( H2_He_coll_fit_par[init][final][6]*t3+1. ); 
    717725                /* this protects against overflow */ 
    718726                if( fabs(sum1)< 38. ) 
    …… 
    721729                        * bug in the pgcc compiler that caused h2_pdr_leiden_f1 to throw an fpe 
    722730                        * the changed code is equivalent */ 
    723                         /*f2 = H2_He_coll_fit_par[init][final][5]/pow(H2_He_coll_fit_par[init][final][6]*t2+1.,H2_He_coll_fit_par[init][final][7]);*/ 
     731                        /*f2 = H2_He_coll_fit_par[init][final][5]/pow(H2_He_coll_fit_par[init][final][6]*t3+1.,H2_He_coll_fit_par[init][final][7]);*/ 
    724732                        f2 = H2_He_coll_fit_par[init][final][5]/pow(10. , sum1 ); 
    725733                }