Changeset 2127 for trunk/source/helike_recom.cpp

Timestamp:

06/25/08 11:17:40 (5 months ago)

Author:

rporter

Message:

/trunk/source - a bunch of relatively minor changes, essentially all refactoring and comments. One change to command parser. "atom XX-like error generation" now accepts "PESSimistic" option, which chooses the worst of two sets of data in the Monte Carlo calculation.

Files:

• trunk/source/helike_recom.cpp (modified) (10 diffs)

trunk/source/helike_recom.cpp

@@ -4 +4 @@
 /*cross_section - calculates the photoionization cross_section for a given level and photon energy*/
 /*radrecomb - calculates radiative recombination coefficients. */
-/*RecomInt - Integral in milne relation.  Called by qg32.       */
 /*He_cross_section returns cross section (cm^-2), 
  * given EgammaRyd, the photon energy in Ryd,
@@ -31 +30 @@
 STATIC double X1Int( double u );
 STATIC double X2Int( double u );
-
-STATIC double cross_section(double EgammaRyd);
-
-static double RecomInt(double EE);
+STATIC double cross_section(double EgammaRyd, long nelem, long ipLev);
 
 #if 0
@@ -46 +42 @@
 #endif
 
-static double kTRyd,EthRyd,Xn_S59; 
-static long int ipLev,globalZ;
+static double EthRyd,Xn_S59; 
 
 /* Here are several structures containing fits to cross-sections calculated in the following: */
@@ -1512 +1507 @@
 
         /* These corrections ensure helium cross-sections agree with the highly accurate 
-         * values given by Hummer and Storey (1998).  The index is [ipLev - 1]  */
+         * values given by Hummer and Storey (1998).  The index is [ipLevel - 1]        */
         /* >>refer He   PCS     Hummer, D.G., & Storey, P.J. 1998, MNRAS 297, 1073      */  
         double PreferredThres[NUM_PREF_TCS] = {
@@ -1527 +1522 @@
                 13.54,26.58,55.71,-1.00,-1.00,-1.00,-1.00,-1.00,-1.00,58.07};
 
-        globalZ = nelem;
-
-        ipLev = ipLevel;
-
-        EthRyd = iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLev];
+        EthRyd = iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLevel];
 
         /* make sure this routine not called within collapsed high levels */
@@ -1537 +1528 @@
 
         /* this is a resolved level */
-        cs = (1.e-18)*cross_section( EgammaRyd );
+        cs = (1.e-18)*cross_section( EgammaRyd, nelem, ipLevel );
 
         /* Scale to preferred threshold values */
@@ -1546 +1537 @@
                 if( ipCurrentLevel!=ipLevel )
                 {
-                        ThreshCS = cross_section( EthRyd );
+                        ThreshCS = cross_section( EthRyd, nelem, ipLevel );
                         ipCurrentLevel = ipLevel;
                 }
@@ -1560 +1551 @@
 
 /*cross_section calculates the photoionization cross_section for a given level and photon energy*/
-STATIC double cross_section(double EgammaRyd)
+STATIC double cross_section(double EgammaRyd, long nelem, long ipLev)
 {
         /* These fit parameters (E0, sigma, y_a, P, y_w, yzero, and yone) all come from the following work: */
@@ -1627 +1618 @@
 
         double pcs,Egamma,lnE,a,b,c,y,F,x;
-        long nelem = globalZ;
-
         int ipISO=ipHE_LIKE;
 
@@ -2038 +2027 @@
 #endif
 
-/*helike_radrecomb_from_cross_section calculates radiative recombination coefficients. */
-double helike_radrecomb_from_cross_section(double temp, long nelem, long ipLo)
-{
-        double alpha,RecomIntegral=0.,b,E1,E2,step,OldRecomIntegral,TotChangeLastFive;
-        double change[5] = {0.,0.,0.,0.,0.};
-
-        EthRyd = iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLo];
-
-        /* Factors outside integral in Milne relation.  */
-        b = MILNE_CONST * StatesElem[ipHE_LIKE][nelem][ipLo].g * pow(temp,-1.5) / 4.;
-        /* kT in Rydbergs.      */
-        kTRyd = temp / TE1RYD;
-        globalZ = nelem;
-        ipLev = ipLo;
-
-        /* Begin integration.   */
-        /* First define characteristic step */
-        E1 = EthRyd;
-        step = MIN2( 0.25*kTRyd, 0.5*E1 );
-        E2 = E1 + step;
-        /* Perform initial integration, from threshold to threshold + step.     */
-        RecomIntegral = qg32( E1, E2, RecomInt);
-        /* Repeat the integration, adding each new result to the total, 
-         * except that the step size is doubled every time, since values away from 
-         * threshold tend to fall off more slowly.      */
-        do
-        {
-                OldRecomIntegral = RecomIntegral;
-                E1 = E2;
-                step *= 1.25;
-                E2 = E1 + step;
-                RecomIntegral += qg32( E1, E2, RecomInt);
-                change[4] = change[3];
-                change[3] = change[2];
-                change[2] = change[1];
-                change[1] = change[0];
-                change[0] = (RecomIntegral - OldRecomIntegral)/RecomIntegral;
-                TotChangeLastFive = change[0] + change[1] + change[2] + change[3] + change[4];
-        /* Continue integration until the upper limit exceeds 100kTRyd, an arbitrary
-         * point at which the integrand component exp(electron energy/kT) is very small,
-         * making neglible cross-sections at photon energies beyond that point,
-         * OR when the last five steps resulted in less than a 1 percent change.        */
-        } while ( ((E2-EthRyd) < 100.*kTRyd) && ( TotChangeLastFive > 0.0001) );
-
-        /* Calculate recombination coefficient. */
-        alpha = b * RecomIntegral;
-
-        alpha = MAX2( alpha, SMALLDOUBLE );
-
-        return alpha;
-}
-
-/*RecomInt, used in comput milne relation for he rec - the energy is photon Rydbergs.   */
-static double RecomInt(double ERyd)
-{
-        double x1, temp;
-
-        /* Milne relation integrand     */
-        x1 = ERyd * ERyd * exp(-1.0 * ( ERyd - EthRyd ) / kTRyd);
-        temp = He_cross_section( ERyd , ipLev, globalZ );
-        x1 *= temp;
-
-        return x1;
-}
-
 /* >>refer      He-like RR      Seaton, M.J. 1959, MNRAS 119, 81S */
 double Recomb_Seaton59( long nelem, double temp, long n)