Changeset 2109 for branches/c08_branch/source/grains.cpp

Timestamp:

06/13/08 22:15:29 (7 months ago)

Author:

peter

Message:

Merge r2108 from mainline - speed up grain code.

Files:

• branches/c08_branch/source/grains.cpp (modified) (14 diffs)

branches/c08_branch/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);
                                 }
                         }