Changeset 2108

Timestamp:

06/13/08 11:23:57 (5 months ago)

Author:

peter

Message:

source/grains_mie.cpp:
source/grains.cpp:
source/grainvar.h:

Speed up the code for initializing the grain electron yields. No change in functionality.

Location:

trunk/source

Files:

• grains.cpp (modified) (14 diffs)
• grains_mie.cpp (modified) (3 diffs)
• grainvar.h (modified) (1 diff)

trunk/source/grains.cpp

@@ -33 +33 @@
 inline double ASINH(double x)
 {
-        if( abs(x) <= 2.3e-4 )
-                return x - pow3(x)/6.;
-        else if( abs(x) <= 8.e-3 )
+        if( abs(x) <= 8.e-3 )
                 return ((0.075*pow2(x) - 1./6.)*pow2(x) + 1.)*x;
         else if( abs(x) <= 1./sqrt(DBL_EPSILON) )
@@ -47 +45 @@
         {
                 if( x < 0. )
-                        return -log(-2.*x);
+                        return -(log(-x)+LN_TWO);
                 else
-                        return log(2.*x);
+                        return log(x)+LN_TWO;
         }
 }
@@ -145 +143 @@
                 return 1.e-7;
         else
-                return 3.e-6*pow((double)gv.bin[nd]->dustp[0],-0.85)*pow(e*EVRYD*1.e-3,1.5);
+                return 3.e-6*gv.bin[nd]->eec*sqrt(pow3(e*EVRYD*1.e-3));
 }
         
@@ -192 +190 @@
 STATIC void UpdatePot2(long,long);
 /* Helper function to calculate primary and secondary yields and the average electron energy at infinity */
-STATIC void Yfunc(long,long,double,double,double,double,double,/*@out@*/double*,/*@out@*/double*,
+inline void Yfunc(long,long,double,double,double,double,double,/*@out@*/double*,/*@out@*/double*,
                   /*@out@*/double*,/*@out@*/double*);
 /* This calculates the y0 function for band electrons (Sect. 4.1.3/4.1.4 of WDB06) */
@@ -203 +201 @@
 STATIC double y1psa(long,long,double);
 /* This calculates the y2 function for primary and Auger electrons (Eq. 8 of WDB06) */
-STATIC double y2pa(double,double,long,double*);
+inline double y2pa(double,double,long,double*);
 /* This calculates the y2 function for secondary electrons (Eq. 20-21 of WDB06) */
-STATIC double y2s(double,double,long,double*);
+inline double y2s(double,double,long,double*);
 /* find highest ionization stage with non-zero population */
 STATIC long HighestIonStage(void);
@@ -3216 +3214 @@
         }
 
+        double GrainPot = chrg2pot(Zg,nd);
+
         if( nfill > ipLo )
         {
@@ -3241 +3241 @@
                         /* calculate yield for band structure */
                         Ehi = Ehi_band = anu[i] - ThresInfVal - Emin;
-                        Wcorr = ThresInfVal + Emin - chrg2pot(Zg,nd);
+                        Wcorr = ThresInfVal + Emin - GrainPot;
                         Eel = anu[i] - DustWorkFcn;
                         yzero = y0b( nd, nz, i );
@@ -3252 +3252 @@
 
                         /* add in yields for inner shell ionization */
-                        for( ns=1; ns < gv.bin[nd]->nShells; ns++ )
+                        unsigned long nsmax = gv.bin[nd]->nShells;
+                        for( ns=1; ns < nsmax; ns++ )
                         {
                                 sptr = gv.bin[nd]->sd[ns];
@@ -3268 +3269 @@
 
                                 /* add in Auger yields */
-                                for( n=0; n < sptr->nData; n++ )
+                                long nmax = sptr->nData;
+                                for( n=0; n < nmax; n++ )
                                 {
                                         double max = sptr->AvNr[n]*sptr->p[i];
-                                        Ehi = sptr->Ener[n] - chrg2pot(Zg,nd);
+                                        Ehi = sptr->Ener[n] - GrainPot;
                                         Eel = sptr->Ener[n];
                                         Yfunc(nd,nz,sptr->y01A[n][i],max,Elo,Ehi,Eel,&Yp1,&Ys1,&Ehp,&Ehs);
@@ -3432 +3434 @@
 
 /* Helper function to calculate primary and secondary yields and the average electron energy at infinity */
-STATIC void Yfunc(long nd,
+inline void Yfunc(long nd,
                   long nz,
                   double y01,
@@ -3472 +3474 @@
                 /* this is Eq. 18 of WDB06 */
                 /* Eel may be negative near threshold -> set yield to zero */
-                f3 = max(Eel,0.)/(eps*elec_esc_length(Eel,nd)*(1./gv.bin[nd]->AvRadius + 1.e7));
+                f3 = max(Eel,0.)/(eps*elec_esc_length(Eel,nd)*gv.bin[nd]->eyc);
                 *Ys = y2sec*f3*min(y01,maxval);
         }
@@ -3614 +3616 @@
 
 /* This calculates the y2 function for primary and Auger electrons (Eq. 8 of WDB06) */
-STATIC double y2pa(double Elo,
+inline double y2pa(double Elo,
                    double Ehi,
                    long Zg,
@@ -3658 +3660 @@
 
 /* This calculates the y2 function for secondary electrons (Eqs. 20-21 of WDB06) */
-STATIC double y2s(double Elo,
+inline double y2s(double Elo,
                   double Ehi,
                   long Zg,
@@ -3713 +3715 @@
                                 {
                                         N0 = alpha*(1./sqR0 - 1./sqRh);
-                                        E0 = alpha*ETILDE*(ASINH(x) + ASINH(yl) - yh/sqRh);
+                                        E0 = alpha*ETILDE*(ASINH(x*sqR0 + yl*sqRh) - yh/sqRh);
                                 }
                         }

trunk/source/grains_mie.cpp

@@ -1088 +1088 @@
         mie_next_data(chFile,io2,chLine,&dl);
         mie_read_realnum(chFile,chLine,&gv.bin[nd]->dustp[0],true,dl);
+        /* constant needed in the evaluation of the electron escape length */
+        gv.bin[nd]->eec = pow((double)gv.bin[nd]->dustp[0],-0.85);
 
         /* molecular weight of grain molecule (amu) */
@@ -1108 +1110 @@
         mie_next_data(chFile,io2,chLine,&dl);
         mie_read_realnum(chFile,chLine,&gv.bin[nd]->AvRadius,true,dl);
+        gv.bin[nd]->eyc = 1./gv.bin[nd]->AvRadius + 1.e7;
 
         /* average grain area <4pi*a^2> for entire size distr (cm^2) */
@@ -1279 +1282 @@
                                         gv.bin[nd2]->dustp[j] = gv.bin[nd]->dustp[j];
                                 }
+                                gv.bin[nd2]->eec = gv.bin[nd]->eec;
+                                gv.bin[nd2]->eyc = gv.bin[nd]->eyc;
                                 gv.bin[nd2]->DustWorkFcn = gv.bin[nd]->DustWorkFcn;
                                 gv.bin[nd2]->BandGap = gv.bin[nd]->BandGap;

trunk/source/grainvar.h

@@ -240 +240 @@
         /** general information on the grains */
         char chDstLab[13];      /**< label for the species */
+        double eec;             /**< pow(dustp[0],-0.85), needed for electron esacpe length */
+        double eyc;            /**< 1./AvRadius + 1.e7, needed for electron yield */
         realnum dustp[5],       /**< 0 = specific weight (g/cm^3), 1 = mol. weight (amu), 2 = default abundance,
                                  *   3 = default depletion, 4 = fraction of the mass in this grain bin */