Changeset e7678b2 in sasmodels

Timestamp:

Feb 29, 2016 6:21:55 AM (8 years ago)

Author:

piotr

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

73860b6

Parents:

deac08c

Message:

Code review from PAK

Location:

sasmodels

Files:

• kernel_template.c (modified) (1 diff)
• models/core_shell_bicelle.c (modified) (4 diffs)
• models/core_shell_bicelle.py (modified) (1 diff)
• models/core_shell_ellipsoid.c (modified) (2 diffs)
• models/core_shell_ellipsoid_xt.c (modified) (3 diffs)
• models/flexible_cylinder.c (modified) (3 diffs)
• models/flexible_cylinder.py (modified) (1 diff)
• models/flexible_cylinder_ex.c (modified) (4 diffs)
• models/flexible_cylinder_ex.py (modified) (1 diff)
• models/lib/J1c.c (added)
• models/lib/Si.c (modified) (1 diff)
• models/lib/wrc_cyl.c (modified) (12 diffs)

sasmodels/kernel_template.c

@@ -88 +88 @@
 #  define M_PI_4 0.7853981633974483
 #endif
+#ifndef M_E
+#  define M_E 2.718281828459045091
+#endif
 
 // Non-standard function library

sasmodels/models/core_shell_bicelle.c

@@ -29 +29 @@
 }
 
+inline double sinc(double x) {return x==0 ? 1.0 : sin(x)/x;}
+
 static double
 bicelle_kernel(double qq,
@@ -41 +43 @@
               double dum)
 {
-        double dr1,dr2,dr3;
-        double besarg1,besarg2;
-        double vol1,vol2,vol3;
-        double sinarg1,sinarg2;
-        double t1,t2,t3;
-        double retval,si1,si2,be1,be2;
+    double si1,si2,be1,be2;
 
-        dr1 = rhoc-rhoh;
-        dr2 = rhor-rhosolv;
-        dr3=  rhoh-rhor;
-        vol1 = M_PI*rad*rad*(2.0*length);
-        vol2 = M_PI*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick);
-        vol3= M_PI*(rad)*(rad)*(2.0*length+2.0*facthick);
-        besarg1 = qq*rad*sin(dum);
-        besarg2 = qq*(rad+radthick)*sin(dum);
-        sinarg1 = qq*length*cos(dum);
-        sinarg2 = qq*(length+facthick)*cos(dum);
+    const double dr1 = rhoc-rhoh;
+    const double dr2 = rhor-rhosolv;
+    const double dr3 = rhoh-rhor;
+    const double vol1 = M_PI*rad*rad*(2.0*length);
+    const double vol2 = M_PI*(rad+radthick)*(rad+radthick)*2.0*(length+facthick);
+    const double vol3 = M_PI*rad*rad*2.0*(length+facthick);
+    double sn,cn;
+    SINCOS(dum, sn, cn);
+    double besarg1 = qq*rad*sn;
+    double besarg2 = qq*(rad+radthick)*sn;
+    double sinarg1 = qq*length*cn;
+    double sinarg2 = qq*(length+facthick)*cn;
 
-        if(besarg1 == 0) {
-                be1 = 0.5;
-        } else {
-                be1 = J1(besarg1)/besarg1;
-        }
-        if(besarg2 == 0) {
-                be2 = 0.5;
-        } else {
-                be2 = J1(besarg2)/besarg2;
-        }
-        if(sinarg1 == 0) {
-                si1 = 1.0;
-        } else {
-                si1 = sin(sinarg1)/sinarg1;
-        }
-        if(sinarg2 == 0) {
-                si2 = 1.0;
-        } else {
-                si2 = sin(sinarg2)/sinarg2;
-        }
-        t1 = 2.0*vol1*dr1*si1*be1;
-        t2 = 2.0*vol2*dr2*si2*be2;
-        t3 = 2.0*vol3*dr3*si2*be1;
+    be1 = J1c(besarg1);
+    be2 = J1c(besarg2);
+    si1 = sinc(sinarg1);
+    si2 = sinc(sinarg2);
 
-        retval = ((t1+t2+t3)*(t1+t2+t3))*sin(dum);
-        return(retval);
+    const double t = vol1*dr1*si1*be1 +
+                     vol2*dr2*si2*be2 +
+                     vol3*dr3*si2*be1;
+
+    const double retval = t*t*sn;
+
+    return(retval);
 
 }
@@ -99 +84 @@
                    double rhosolv)
 {
+    // set up the integration end points
+    const double uplim = M_PI/4;
+    const double halfheight = length/2.0;
 
+    double summ = 0.0;
+    for(int i=0;i<N_POINTS_76;i++) {
+        double zi = (Gauss76Z[i] + 1.0)*uplim;
+        double yyy = Gauss76Wt[i] * bicelle_kernel(qq, rad, radthick, facthick,
+                             halfheight, rhoc, rhoh, rhor,rhosolv, zi);
+        summ += yyy;
+    }
 
-        double answer,halfheight;
-        double lolim,uplim,summ,yyy,zi;
-        int nord,i;
-
-        // set up the integration end points
-        nord = 76;
-        lolim = 0.0;
-        uplim = M_PI/2;
-        halfheight = length/2.0;
-
-        summ = 0.0;
-        i=0;
-        for(i=0;i<nord;i++) {
-                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
-                yyy = Gauss76Wt[i] * bicelle_kernel(qq, rad, radthick, facthick,
-                                     halfheight, rhoc, rhoh, rhor,rhosolv, zi);
-                summ += yyy;
-        }
-
-        // calculate value of integral to return
-        answer = (uplim-lolim)/2.0*summ;
-        return(answer);
+    // calculate value of integral to return
+    double answer = uplim*summ;
+    return(answer);
 }
 
@@ -138 +114 @@
           double phi)
 {
-    double cyl_x, cyl_y;
-    double alpha, cos_val;
-    double answer;
-
     //convert angle degree to radian
-    theta *= M_PI/180.0;
-    phi *= M_PI/180.0;
+    theta *= M_PI_180;
+    phi *= M_PI_180;
 
     // Cylinder orientation
-    cyl_x = cos(theta) * cos(phi);
-    cyl_y = sin(theta);
+    const double cyl_x = cos(theta) * cos(phi);
+    const double cyl_y = sin(theta);
 
     // Compute the angle btw vector q and the axis of the cylinder
-    cos_val = cyl_x*q_x + cyl_y*q_y;
-    alpha = acos( cos_val );
+    const double cos_val = cyl_x*q_x + cyl_y*q_y;
+    const double alpha = acos( cos_val );
 
     // Get the kernel
-    answer = bicelle_kernel(q, radius, rim_thickness, face_thickness,
+    double answer = bicelle_kernel(q, radius, rim_thickness, face_thickness,
                            length/2.0, core_sld, face_sld, rim_sld,
                            solvent_sld, alpha) / fabs(sin(alpha));

sasmodels/models/core_shell_bicelle.py

@@ -91 +91 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/J1.c", "lib/gauss76.c", "core_shell_bicelle.c"]
+source = ["lib/Si.c", "lib/J1.c", "lib/J1c.c", "lib/gauss76.c", "core_shell_bicelle.c"]
 
 demo = dict(scale=1, background=0,

sasmodels/models/core_shell_ellipsoid.c

@@ -44 +44 @@
           double solvent_sld)
 {
-        double delpc,delps;
-        double uplim,lolim;             //upper and lower integration limits
-        double summ,zi,yyy,answer; //running tally of integration
 
-        lolim = 0.0;
-        uplim = 1.0;
+    //upper and lower integration limits
+    const double lolim = 0.0;
+    const double uplim = 1.0;
 
-        summ = 0.0;      //initialize intergral
+    double summ = 0.0;   //initialize intergral
 
-        delpc = core_sld - shell_sld; //core - shell
-        delps = shell_sld - solvent_sld; //shell - solvent
+    const double delpc = core_sld - shell_sld;    //core - shell
+    const double delps = shell_sld - solvent_sld; //shell - solvent
 
-        for(int i=0;i<76;i++) {
-                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
-                yyy = Gauss76Wt[i] * gfn4(zi,
-                                          equat_core,
-                                          polar_core,
-                                          equat_shell,
-                                          polar_shell,
-                                          delpc,
-                                          delps,
-                                          q);
-                summ += yyy;
-        }
+    for(int i=0;i<N_POINTS_76;i++) {
+        double zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
+        double yyy = Gauss76Wt[i] * gfn4(zi,
+                                  equat_core,
+                                  polar_core,
+                                  equat_shell,
+                                  polar_shell,
+                                  delpc,
+                                  delps,
+                                  q);
+        summ += yyy;
+    }
 
-        answer = (uplim-lolim)/2.0*summ;
+    double answer = (uplim-lolim)/2.0*summ;
 
-        //convert to [cm-1]
-        answer *= 1.0e-4;
+    //convert to [cm-1]
+    answer *= 1.0e-4;
 
-        return answer;
+    return answer;
 }
 
@@ -89 +87 @@
           double phi)
 {
-    double cyl_x, cyl_y;
-    double cos_val;
-    double answer;
-    double sldcs,sldss;
-
     //convert angle degree to radian
-    theta = theta * M_PI/180.0;
-    phi = phi * M_PI/180.0;
+    theta = theta * M_PI_180;
+    phi = phi * M_PI_180;
 
 
     // ellipsoid orientation, the axis of the rotation is consistent with the ploar axis.
-    cyl_x = cos(theta) * cos(phi);
-    cyl_y = sin(theta);
+    const double cyl_x = cos(theta) * cos(phi);
+    const double cyl_y = sin(theta);
 
-    sldcs = core_sld - shell_sld;
-    sldss = shell_sld- solvent_sld;
+    const double sldcs = core_sld - shell_sld;
+    const double sldss = shell_sld- solvent_sld;
 
     // Compute the angle btw vector q and the
     // axis of the cylinder
-    cos_val = cyl_x*q_x + cyl_y*q_y;
+    const double cos_val = cyl_x*q_x + cyl_y*q_y;
 
     // Call the IGOR library function to get the kernel: MUST use gfn4 not gf2 because of the def of params.
-    answer = gfn4(cos_val,
+    double answer = gfn4(cos_val,
                   equat_core,
                   polar_core,

sasmodels/models/core_shell_ellipsoid_xt.c

@@ -30 +30 @@
                    double x_polar_shell)
 {
-        double equat_shell, polar_shell;
-    equat_shell = equat_core + t_shell;
-    polar_shell = equat_core*x_core + t_shell*x_polar_shell;
+    const double equat_shell = equat_core + t_shell;
+    const double polar_shell = equat_core*x_core + t_shell*x_polar_shell;
     double vol = 4.0*M_PI/3.0*equat_shell*equat_shell*polar_shell;
     return vol;
@@ -47 +46 @@
           double solvent_sld)
 {
-        double delpc,delps;
-        double uplim,lolim;             //upper and lower integration limits
-        double summ,zi,yyy,answer; //running tally of integration
-        double polar_core, equat_shell, polar_shell;
+    const double lolim = 0.0;
+    const double uplim = 1.0;
 
-        lolim = 0.0;
-        uplim = 1.0;
+    double summ = 0.0;   //initialize intergral
 
-        summ = 0.0;      //initialize intergral
-
-        delpc = core_sld - shell_sld; //core - shell
-        delps = shell_sld - solvent_sld; //shell - solvent
+    const double delpc = core_sld - shell_sld; //core - shell
+    const double delps = shell_sld - solvent_sld; //shell - solvent
 
 
-    polar_core = equat_core*x_core;
-    equat_shell = equat_core + t_shell;
-    polar_shell = equat_core*x_core + t_shell*x_polar_shell;
+    const double polar_core = equat_core*x_core;
+    const double equat_shell = equat_core + t_shell;
+    const double polar_shell = equat_core*x_core + t_shell*x_polar_shell;
 
-        for(int i=0;i<76;i++) {
-                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
-                yyy = Gauss76Wt[i] * gfn4(zi,
-                                          equat_core,
+    for(int i=0;i<N_POINTS_76;i++) {
+        double zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
+        double yyy = Gauss76Wt[i] * gfn4(zi,
+                                  equat_core,
                                   polar_core,
                                   equat_shell,
                                   polar_shell,
-                                          delpc,
-                                          delps,
-                                          q);
-                summ += yyy;
-        }
+                                  delpc,
+                                  delps,
+                                  q);
+        summ += yyy;
+    }
 
-        answer = (uplim-lolim)/2.0*summ;
+    double answer = (uplim-lolim)/2.0*summ;
+    //convert to [cm-1]
+    answer *= 1.0e-4;
 
-        //convert to [cm-1]
-        answer *= 1.0e-4;
-
-        return answer;
+    return answer;
 }
 
@@ -98 +91 @@
           double phi)
 {
-    double cyl_x, cyl_y;
-    double cos_val;
-    double answer;
-    double sldcs,sldss;
-        double polar_core, equat_shell, polar_shell;
-
     //convert angle degree to radian
-    theta = theta * M_PI/180.0;
-    phi = phi * M_PI/180.0;
-
+    theta = theta * M_PI_180;
+    phi = phi * M_PI_180;
 
     // ellipsoid orientation, the axis of the rotation is consistent with the ploar axis.
-    cyl_x = cos(theta) * cos(phi);
-    cyl_y = sin(theta);
+    const double cyl_x = cos(theta) * cos(phi);
+    const double cyl_y = sin(theta);
 
-    sldcs = core_sld - shell_sld;
-    sldss = shell_sld- solvent_sld;
+    const double sldcs = core_sld - shell_sld;
+    const double sldss = shell_sld- solvent_sld;
 
     // Compute the angle btw vector q and the
     // axis of the cylinder
-    cos_val = cyl_x*q_x + cyl_y*q_y;
+    const double cos_val = cyl_x*q_x + cyl_y*q_y;
 
-    polar_core = equat_core*x_core;
-    equat_shell = equat_core + t_shell;
-    polar_shell = equat_core*x_core + t_shell*x_polar_shell;
+    const double polar_core = equat_core*x_core;
+    const double equat_shell = equat_core + t_shell;
+    const double polar_shell = equat_core*x_core + t_shell*x_polar_shell;
 
     // Call the IGOR library function to get the kernel:
     // MUST use gfn4 not gf2 because of the def of params.
-    answer = gfn4(cos_val,
+    double answer = gfn4(cos_val,
                   equat_core,
                   polar_core,

sasmodels/models/flexible_cylinder.c

@@ -22 +22 @@
 {
 
-        double Pi = 4.0*atan(1.0);
-
-        double cont = sld-solvent_sld;
-        double qr = q*radius;
-        double flex = Sk_WR(q,length,kuhn_length);
-    double crossSect = (2.0*J1(qr)/qr)*(2.0*J1(qr)/qr);
-
-        flex *= crossSect;
-        flex *= Pi*radius*radius*length;
-        flex *= cont*cont;
-        flex *= 1.0e-4;
-
-        return flex;
+    const double cont = sld-solvent_sld;
+    const double qr = q*radius;
+    //const double crossSect = (2.0*J1(qr)/qr)*(2.0*J1(qr)/qr);
+    const double crossSect = J1c(qr);
+    double flex = Sk_WR(q,length,kuhn_length);
+    flex *= crossSect*crossSect;
+    flex *= M_PI*radius*radius*length;
+    flex *= cont*cont;
+    flex *= 1.0e-4;
+    return flex;
 }
 
@@ -45 +42 @@
 {
 
-        double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld);
-
-        return result;
+    double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld);
+    return result;
 }
 
@@ -57 +53 @@
             double solvent_sld)
 {
-        double q;
-        q = sqrt(qx*qx+qy*qy);
-        double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld);
+    double q;
+    q = sqrt(qx*qx+qy*qy);
+    double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld);
 
-        return result;
+    return result;
 }

sasmodels/models/flexible_cylinder.py

@@ -90 +90 @@
     ]
 # pylint: enable=bad-whitespace, line-too-long
-source = ["lib/J1.c", "lib/wrc_cyl.c", "flexible_cylinder.c"]
+source = ["lib/J1.c", "lib/J1c.c", "lib/wrc_cyl.c", "flexible_cylinder.c"]
 
 demo = dict(scale=1.0, background=0.0001,

sasmodels/models/flexible_cylinder_ex.c

@@ -17 +17 @@
 elliptical_crosssection(double q, double a, double b)
 {
-    double uplim,lolim,Pi,summ,arg,zi,yyy,answer;
-    int i,nord=76;
+    double summ=0.0;
 
-    Pi = 4.0*atan(1.0);
-    lolim=0.0;
-    uplim=Pi/2.0;
-    summ=0.0;
+    for(int i=0;i<N_POINTS_76;i++) {
+        double zi = ( Gauss76Z[i] + 1.0 )*M_PI/4.0;
+        double sn, cn;
+        SINCOS(zi, sn, cn);
+        double arg = q*sqrt(a*a*sn*sn+b*b*cn*cn);
+        double yyy = pow((double)J1c(arg),2);
+        yyy *= Gauss76Wt[i];
+        summ += yyy;
+    }
 
-    for(i=0;i<nord;i++) {
-                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
-                arg = q*sqrt(a*a*sin(zi)*sin(zi)+b*b*cos(zi)*cos(zi));
-                yyy = pow((2.0 * J1(arg) / arg),2);
-                yyy *= Gauss76Wt[i];
-                summ += yyy;
-    }
-    answer = (uplim-lolim)/2.0*summ;
-    answer *= 2.0/Pi;
-    return(answer);
+    summ /= 2.0;
+    return(summ);
 
 }
@@ -47 +43 @@
 {
 
-        double Pi,flex,crossSect, cont;
+    double flex,crossSect, cont;
 
-        Pi = 4.0*atan(1.0);
-        cont = sld - solvent_sld;
-        crossSect = elliptical_crosssection(q,radius,(radius*axis_ratio));
+    cont = sld - solvent_sld;
+    crossSect = elliptical_crosssection(q,radius,(radius*axis_ratio));
 
-        flex = Sk_WR(q,length,kuhn_length);
-        flex *= crossSect;
-        flex *= Pi*radius*radius*axis_ratio*axis_ratio*length;
-        flex *= cont*cont;
-        flex *= 1.0e-4;
+    flex = Sk_WR(q,length,kuhn_length);
+    flex *= crossSect;
+    flex *= M_PI*radius*radius*axis_ratio*axis_ratio*length;
+    flex *= cont*cont;
+    flex *= 1.0e-4;
 
-        return flex;
+    return flex;
 }
 
@@ -71 +66 @@
 {
 
-        double result = flexible_cylinder_ex_kernel(q,
-                        length,
-                        kuhn_length,
-                        radius,
-                        axis_ratio,
-                        sld,
-                        solvent_sld);
+    double result = flexible_cylinder_ex_kernel(q,
+                    length,
+                    kuhn_length,
+                    radius,
+                    axis_ratio,
+                    sld,
+                    solvent_sld);
 
-        return result;
+    return result;
 }
 
@@ -90 +85 @@
             double solvent_sld)
 {
-        double q;
-        q = sqrt(qx*qx+qy*qy);
-        double result = flexible_cylinder_ex_kernel(q,
-                        length,
-                        kuhn_length,
-                        radius,
-                        axis_ratio,
-                        sld,
-                        solvent_sld);
+    double q;
+    q = sqrt(qx*qx+qy*qy);
+    double result = flexible_cylinder_ex_kernel(q,
+                    length,
+                    kuhn_length,
+                    radius,
+                    axis_ratio,
+                    sld,
+                    solvent_sld);
 
-        return result;
+    return result;
 }

sasmodels/models/flexible_cylinder_ex.py

@@ -113 +113 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/J1.c", "lib/gauss76.c", "lib/wrc_cyl.c", "flexible_cylinder_ex.c"]
+source = ["lib/J1.c", "lib/J1c.c", "lib/gauss76.c", "lib/wrc_cyl.c", "flexible_cylinder_ex.c"]
 
 demo = dict(scale=1.0, background=0.0001,

sasmodels/models/lib/Si.c

@@ -4 +4 @@
 double Si(double x)
 {
-        double out;
-        double pi = 4.0*atan(1.0);
+    double out;
 
-        if (x >= pi*6.2/4.0){
-                double out_sin = 0.0;
-                double out_cos = 0.0;
-                out = pi/2.0;
+    if (x >= M_PI*6.2/4.0){
+        double out_sin = 0.0;
+        double out_cos = 0.0;
+        out = M_PI/2.0;
 
-                // Explicitly writing factorial values triples the speed of the calculation
-                out_cos = 1./x - 2./pow(x,3) + 24./pow(x,5) - 720./pow(x,7);
-                out_sin = 1./pow(x,2) - 6./pow(x,4) + 120./pow(x,6) - 5040./pow(x,8);
+        // Explicitly writing factorial values triples the speed of the calculation
+        out_cos = 1./x - 2./pow(x,3) + 24./pow(x,5) - 720./pow(x,7);
+        out_sin = 1./pow(x,2) - 6./pow(x,4) + 120./pow(x,6) - 5040./pow(x,8);
 
-                out -= cos(x) * out_cos;
-                out -= sin(x) * out_sin;
-                return out;
-        }
+        out -= cos(x) * out_cos;
+        out -= sin(x) * out_sin;
+        return out;
+    }
 
-        // Explicitly writing factorial values triples the speed of the calculation
-        out = x - pow(x, 3)/18. + pow(x,5)/600. - pow(x,7)/35280. + pow(x,9)/3265920. - pow(x,11)/439084800.;
+    // Explicitly writing factorial values triples the speed of the calculation
+    out = x - pow(x, 3)/18. + pow(x,5)/600. - pow(x,7)/35280. + pow(x,9)/3265920. - pow(x,11)/439084800.;
 
-        return out;
+    return out;
 }

sasmodels/models/lib/wrc_cyl.c

@@ -5 +5 @@
 AlphaSquare(double x)
 {
+    // Potentially faster. Needs proper benchmarking.
+    // add native_powr to kernel_template
+    //double t = native_powr( (1.0 + (x/3.12)*(x/3.12) +
+    //     (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) );
+    //return t;
+
     return pow( (1.0 + (x/3.12)*(x/3.12) +
          (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) );
@@ -13 +19 @@
 Rgsquarezero(double q, double L, double b)
 {
-    return (L*b/6.0) * (1.0 - 1.5*(b/L) + 1.5*pow((b/L),2) -
-         0.75*pow((b/L),3)*(1.0 - exp(-2.0*(L/b))));
+    const double r = b/L;
+    return (L*b/6.0) *
+           (1.0 - r*1.5  + 1.5*r*r - 0.75*r*r*r*(1.0 - exp(-2.0/r)));
 }
 
@@ -31 +38 @@
 }
 
-static double
+static inline double
 sech_WR(double x)
 {
@@ -40 +47 @@
 a1long(double q, double L, double b, double p1, double p2, double q0)
 {
-    double yy,C,C1,C2,C3,C4,C5,miu,Rg2;
-    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15;
-    double E,pi;
-    double b2,b3,b4,q02,q03,q04,q05,Rg22;
-
-    pi = 4.0*atan(1.0);
-    E = 2.718281828459045091;
+    double C;
+    const double onehalf = 1.0/2.0;
 
     if( L/b > 10.0) {
@@ -54 +56 @@
     }
 
-    C1 = 1.22;
-    C2 = 0.4288;
-    C3 = -1.651;
-    C4 = 1.523;
-    C5 = 0.1477;
-    miu = 0.585;
-
-    Rg2 = Rgsquare(q,L,b);
-    Rg22 = Rg2*Rg2;
-    b2 = b*b;
-    b3 = b*b*b;
-    b4 = b3*b;
-    q02 = q0*q0;
-    q03 = q0*q0*q0;
-    q04 = q03*q0;
-    q05 = q04*q0;
-
-    t1 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 +
-        (7.0*b2)/(15.0*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))));
-
-    t2 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) +
-        (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) -
-        tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
-
-    t3 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) +
-        C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) +
-        C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))));
-
-    t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
-
-    t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) -
-        b*p2*pow(q0,((-1.0) - p1 - p2))));
-
-    t6 = (b*C*(((-((14.0*b3)/(15.0*q03*Rg2))) +
-        (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) +
-        (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 +
-        (7.0*b2)/(15.0*q02*Rg2)))*Rg2)/b)));
-
-    t7 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) +
-        C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) +
-        C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) +
-        (sqrt(Rg2)*q0)/b)/C5),2));
-
-    t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) +
-        (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2));
-
-    t9 = (2.0*b4*(((2.0*q0*Rg2)/b -
-        (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) -
-        tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
-
-    t10 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) +
-        (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) +
-        (sqrt(Rg2)*q0)/b)/C5))))));
-
-    t11 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) -
-        3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) -
-        2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) -
-        1.0/miu)))/miu));
-
-    t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
-
-    t13 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 +
-        (7.0*b2)/(15.0*q02* Rg2))) +
-        (7.0*b2)/(15.0*q02*Rg2))));
-
-    t14 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) +
-        (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) +
-        (sqrt(Rg2)*q0)/b)/C5))))));
-
-    t15 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) +
-        C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) +
-        C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))));
-
-
-    yy = (pow(q0,p1)*(((-((b*pi)/(L*q0))) +t1/L +t2/(q04*Rg22) +
-        1.0/2.0*t3*t4)) + (t5*((pow(q0,(p1 - p2))*
-        (((-pow(q0,(-p1)))*(((b2*pi)/(L*q02) +t6/L +t7/(2.0*C5) -
-        t8/(C5*q04*Rg22) + t9/(q04*Rg22) -t10/(q05*Rg22) + 1.0/2.0*t11*t12)) -
-        b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) + t13/L +
-        t14/(q04*Rg22) + 1.0/2.0*t15*((1.0 + tanh(((-C4) +
-        (sqrt(Rg2)*q0)/b)/C5)))))))))));
+    const double C1 = 1.22;
+    const double C2 = 0.4288;
+    const double C3 = -1.651;
+    const double C4 = 1.523;
+    const double C5 = 0.1477;
+    const double miu = 0.585;
+
+    const double Rg2 = Rgsquare(q,L,b);
+    const double Rg22 = Rg2*Rg2;
+    const double Rg = sqrt(Rg2);
+    const double Rgb = Rg*q0/b;
+
+    const double b2 = b*b;
+    const double b3 = b*b*b;
+    const double b4 = b3*b;
+    const double q02 = q0*q0;
+    const double q03 = q0*q0*q0;
+    const double q04 = q03*q0;
+    const double q05 = q04*q0;
+
+    const double Rg02 = Rg2*q02;
+
+    const double t1 = (b*C*((4.0/15.0 - pow((double)M_E,(-(Rg02/b2))) *
+         ((11.0/15.0 + (7.0*b2)/(15.0*Rg02))) +
+         (7.0*b2)/(15.0*Rg02))));
+
+    const double t2 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) +
+         Rg02/b2))*((1.0 + onehalf*(((-1.0) -
+         tanh((-C4 + Rgb/C5)))))));
+
+    const double t3 = ((C3*pow(Rgb,((-3.0)/miu)) +
+         C2*pow(Rgb,((-2.0)/miu)) +
+         C1*pow(Rgb,((-1.0)/miu))));
+
+    const double t4 = ((1.0 + tanh(((-C4) + Rgb)/C5)));
+
+    const double t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) -
+         b*p2*pow(q0,((-1.0) - p1 - p2))));
+
+    const double t6 = (b*C*(((-((14.0*b3)/(15.0*q03*Rg2))) +
+         (14.0*b3*pow((double)M_E,(-(Rg02/b2))))/(15.0*q03*Rg2) +
+         (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*((11.0/15.0 +
+         (7.0*b2)/(15.0*Rg02)))*Rg2)/b)));
+
+    const double t7 = (Rg*((C3*pow(((Rgb)),((-3.0)/miu)) +
+         C2*pow(((Rgb)),((-2.0)/miu)) +
+         C1*pow(((Rgb)),((-1.0)/miu))))*pow(sech_WR(((-C4) +
+         Rgb)/C5),2));
+
+    const double t8 = (b4*Rg*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) +
+         Rg02/b2))*pow(sech_WR(((-C4) + Rgb)/C5),2));
+
+    const double t9 = (2.0*b4*(((2.0*q0*Rg2)/b -
+         (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))*((1.0 + onehalf*(((-1.0) -
+         tanh(((-C4) + Rgb)/C5))))));
+
+    const double t10 = (8.0*b4*b*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) +
+         Rg02/b2))*((1.0 + onehalf*(((-1.0) - tanh(((-C4) +
+         Rgb)/C5))))));
+
+    const double t11 = (((-((3.0*C3*Rg*pow(((Rgb)),((-1.0) -
+          3.0/miu)))/miu)) - (2.0*C2*Rg*pow(((Rgb)),((-1.0) -
+          2.0/miu)))/miu - (C1*Rg*pow(((Rgb)),((-1.0) -
+          1.0/miu)))/miu));
+
+    const double t12 = ((1.0 + tanh(((-C4) + Rgb)/C5)));
+
+    const double t13 = (b*C*((4.0/15.0 - pow((double)M_E,(-(Rg02/b2)))*((11.0/15.0 +
+          (7.0*b2)/(15.0*q02* Rg2))) +
+          (7.0*b2)/(15.0*Rg02))));
+
+    const double t14 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) +
+          Rg02/b2))*((1.0 + onehalf*(((-1.0) - tanh(((-C4) +
+          Rgb)/C5))))));
+
+    const double t15 = ((C3*pow(((Rgb)),((-3.0)/miu)) +
+        C2*pow(((Rgb)),((-2.0)/miu)) +
+        C1*pow(((Rgb)),((-1.0)/miu))));
+
+
+    double yy = (pow(q0,p1)*(((-((b*M_PI)/(L*q0))) +t1/L +t2/(q04*Rg22) +
+        onehalf*t3*t4)) + (t5*((pow(q0,(p1 - p2))*
+        (((-pow(q0,(-p1)))*(((b2*M_PI)/(L*q02) +t6/L +t7/(2.0*C5) -
+        t8/(C5*q04*Rg22) + t9/(q04*Rg22) -t10/(q05*Rg22) + onehalf*t11*t12)) -
+        b*p1*pow(q0,((-1.0) - p1))*(((-((b*M_PI)/(L*q0))) + t13/L +
+        t14/(q04*Rg22) + onehalf*t15*((1.0 + tanh(((-C4) +
+        Rgb)/C5)))))))))));
 
     return (yy);
@@ -142 +150 @@
 a2long(double q, double L, double b, double p1, double p2, double q0)
 {
-    double yy,C1,C2,C3,C4,C5,miu,C,Rg2;
-    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi;
-    double E,b2,b3,b4,q02,q03,q04,q05,Rg22;
-
-    pi = 4.0*atan(1.0);
-    E = 2.718281828459045091;
+    double C;
+    const double onehalf = 1.0/2.0;
+
     if( L/b > 10.0) {
         C = 3.06/pow((L/b),0.44);
@@ -154 +159 @@
     }
 
-    C1 = 1.22;
-    C2 = 0.4288;
-    C3 = -1.651;
-    C4 = 1.523;
-    C5 = 0.1477;
-    miu = 0.585;
-
-    Rg2 = Rgsquare(q,L,b);
-    Rg22 = Rg2*Rg2;
-    b2 = b*b;
-    b3 = b*b*b;
-    b4 = b3*b;
-    q02 = q0*q0;
-    q03 = q0*q0*q0;
-    q04 = q03*q0;
-    q05 = q04*q0;
-
-    t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) -
-        b*p2*pow(q0,((-1.0) - p1 - p2)) ));
-
-    t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) +
-        (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) +
-        (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 +
-        (7*b2)/(15.0*q02*Rg2)))*Rg2)/b)))/L;
-
-    t3 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) +
-        C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) +
-        C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*
-        pow(sech_WR(((-C4) +(sqrt(Rg2)*q0)/b)/C5),2.0))/(2.0*C5);
-
-    t4 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) +
-        (q02*Rg2)/b2))*pow(sech_WR(((-C4) +
-        (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22);
-
-    t5 = (2.0*b4*(((2.0*q0*Rg2)/b -
-        (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*
-        ((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) +
-        (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22);
-
-    t6 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) +
-        (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) +
-        (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22);
-
-    t7 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) -
-        3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),
-        ((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),
-        ((-1.0) - 1.0/miu)))/miu));
-
-    t8 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
-
-    t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 +
-        (7.0    *b2)/(15*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))))/L;
-
-    t10 = (2.0*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) +
-        (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) +
-        (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22);
-
-    yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) +
-        t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))*
-        (((-((b*pi)/(L*q0))) + t9 + t10 +
-        1.0/2.0*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) +
-        C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) +
-        C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*
-        ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))));
+    const double C1 = 1.22;
+    const double C2 = 0.4288;
+    const double C3 = -1.651;
+    const double C4 = 1.523;
+    const double C5 = 0.1477;
+    const double miu = 0.585;
+
+    const double Rg2 = Rgsquare(q,L,b);
+    const double Rg22 = Rg2*Rg2;
+    const double b2 = b*b;
+    const double b3 = b*b*b;
+    const double b4 = b3*b;
+    const double q02 = q0*q0;
+    const double q03 = q0*q0*q0;
+    const double q04 = q03*q0;
+    const double q05 = q04*q0;
+    const double Rg = sqrt(Rg2);
+    const double Rgb = Rg*q0/b;
+    const double Rg02 = Rg2*q02;
+
+    const double t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) -
+         b*p2*pow(q0,((-1.0) - p1 - p2)) ));
+
+    const double t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) +
+         (14.0*b3*pow((double)M_E,(-(Rg02/b2))))/(15.0*q03*Rg2) +
+         (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*((11.0/15.0 +
+         (7*b2)/(15.0*Rg02)))*Rg2)/b)))/L;
+
+    const double t3 = (Rg*((C3*pow(((Rgb)),((-3.0)/miu)) +
+         C2*pow(((Rgb)),((-2.0)/miu)) +
+         C1*pow(((Rgb)),((-1.0)/miu))))*
+         pow(sech_WR(((-C4) +Rgb)/C5),2.0))/(2.0*C5);
+
+    const double t4 = (b4*Rg*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) +
+         Rg02/b2))*pow(sech_WR(((-C4) +
+         Rgb)/C5),2))/(C5*q04*Rg22);
+
+    const double t5 = (2.0*b4*(((2.0*q0*Rg2)/b -
+         (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))*
+         ((1.0 + onehalf*(((-1.0) - tanh(((-C4) +
+         Rgb)/C5))))))/(q04*Rg22);
+
+    const double t6 = (8.0*b4*b*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) +
+         Rg02/b2))*((1.0 + onehalf*(((-1) - tanh(((-C4) +
+         Rgb)/C5))))))/(q05*Rg22);
+
+    const double t7 = (((-((3.0*C3*Rg*pow(((Rgb)),((-1.0) -
+         3.0/miu)))/miu)) - (2.0*C2*Rg*pow(((Rgb)),
+         ((-1.0) - 2.0/miu)))/miu - (C1*Rg*pow(((Rgb)),
+         ((-1.0) - 1.0/miu)))/miu));
+
+    const double t8 = ((1.0 + tanh(((-C4) + Rgb)/C5)));
+
+    const double t9 = (b*C*((4.0/15.0 - pow((double)M_E,(-(Rg02/b2)))*((11.0/15.0 +
+         (7.0*b2)/(15*Rg02))) + (7.0*b2)/(15.0*Rg02))))/L;
+
+    const double t10 = (2.0*b4*(((-1) + pow((double)M_E,(-(Rg02/b2))) +
+          Rg02/b2))*((1.0 + onehalf*(((-1) - tanh(((-C4) +
+          Rgb)/C5))))))/(q04*Rg22);
+
+    const double yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*M_PI)/(L*q02) +
+         t2 + t3 - t4 + t5 - t6 + onehalf*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))*
+         (((-((b*M_PI)/(L*q0))) + t9 + t10 +
+         onehalf*((C3*pow(((Rgb)),((-3.0)/miu)) +
+         C2*pow(((Rgb)),((-2.0)/miu)) +
+         C1*pow(((Rgb)),((-1.0)/miu))))*
+         ((1.0 + tanh(((-C4) + Rgb)/C5))))))))));
 
     return (yy);
@@ -224 +232 @@
 //
 static double
-a1short(double q, double L, double b, double p1short, double p2short, double q0)
-{
-    double yy,Rg2_sh;
-    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3;
-    double pi;
-
-    E = 2.718281828459045091;
-    pi = 4.0*atan(1.0);
-    Rg2_sh = Rgsquareshort(q,L,b);
-    Rg2_sh2 = Rg2_sh*Rg2_sh;
-    b3 = b*b*b;
-    t1 = ((q0*q0*Rg2_sh)/(b*b));
-    Et1 = pow(E,t1);
-    Emt1 =pow(E,-t1);
-    q02 = q0*q0;
-    q0p = pow(q0,(-4.0 + p1short) );
-
-    yy = ((1.0/(L*((p1short - p2short))*Rg2_sh2)*
+a_short(double q, double L, double b, double p1short,
+        double p2short, double factor, double pdiff, double q0)
+{
+    const double Rg2_sh = Rgsquareshort(q,L,b);
+    const double Rg2_sh2 = Rg2_sh*Rg2_sh;
+    const double b3 = b*b*b;
+    const double t1 = ((q0*q0*Rg2_sh)/(b*b));
+    const double Et1 = exp(t1);
+    const double Emt1 = 1.0/Et1;
+    const double q02 = q0*q0;
+    const double q0p = pow(q0,(-4.0 + p1short) );
+
+    double yy = ((factor/(L*(pdiff)*Rg2_sh2)*
         ((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short +
         2.0*b3*Et1*L*p2short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh -
-        2.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 +
-        Et1*p2short*pi*q02*q0*Rg2_sh2))))));
+        2.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*M_PI*q02*q0*Rg2_sh2 +
+        Et1*p2short*M_PI*q02*q0*Rg2_sh2))))));
 
     return(yy);
 }
+static double
+a1short(double q, double L, double b, double p1short, double p2short, double q0)
+{
+
+    double factor = 1.0;
+    return a_short(q, L, b, p1short, p2short, factor, p1short - p2short, q0);
+}
 
 static double
 a2short(double q, double L, double b, double p1short, double p2short, double q0)
 {
-    double yy,Rg2_sh;
-    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p;
-    double pi;
-
-    E = 2.718281828459045091;
-    pi = 4.0*atan(1.0);
-    Rg2_sh = Rgsquareshort(q,L,b);
-    Rg2_sh2 = Rg2_sh*Rg2_sh;
-    t1 = ((q0*q0*Rg2_sh)/(b*b));
-    Et1 = pow(E,t1);
-    Emt1 =pow(E,-t1);
-    q02 = q0*q0;
-    q0p = pow(q0,(-4.0 + p2short) );
-
-    //E is the number e
-    yy = ((-(1.0/(L*((p1short - p2short))*Rg2_sh2)*
-        ((b*Emt1*q0p*((8.0*b*b*b*L - 8.0*b*b*b*Et1*L - 2.0*b*b*b*L*p1short +
-        2.0*b*b*b*Et1*L*p1short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh -
-        2.0*b*Et1*L*p1short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 +
-        Et1*p1short*pi*q02*q0*Rg2_sh2)))))));
-
-    return (yy);
+    double factor = -1.0;
+    return a_short(q, L, b, p2short, p1short, factor, p1short-p2short, q0);
 }
 
@@ -301 +291 @@
 {
     // ORIGINAL
-    double result = 2.0*(exp(-arg) + arg -1.0)/(pow((arg),2));
+    double result = 2.0*(exp(-arg) + arg -1.0)/(arg*arg);
 
     // CONVERSION 1 from http://herbie.uwplse.org/
@@ -333 +323 @@
 Sexv(double q, double L, double b)
 {
-    double yy,C1,C2,C3,miu,Rg2;
-    C1=1.22;
-    C2=0.4288;
-    C3=-1.651;
-    miu = 0.585;
-
-    Rg2 = Rgsquare(q,L,b);
-
-    yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) +
-        w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) +
-        C2*pow((q*sqrt(Rg2)),(-2.0/miu)) +
-        C3*pow((q*sqrt(Rg2)),(-3.0/miu)));
-
+
+    const double C1=1.22;
+    const double C2=0.4288;
+    const double C3=-1.651;
+    const double miu = 0.585;
+    const double qRg = q*sqrt(Rgsquare(q,L,b));
+    const double x = pow(qRg, -1.0/miu);
+
+    double yy = (1.0 - w_WR(qRg))*Sdebye(q,L,b) +
+            w_WR(qRg)*x*(C1 + x*(C2 + x*C3));
     return (yy);
 }
@@ -353 +340 @@
 Sexvnew(double q, double L, double b)
 {
-    double yy,C1,C2,C3,miu;
-    double del=1.05,C_star2,Rg2;
-
-    C1=1.22;
-    C2=0.4288;
-    C3=-1.651;
-    miu = 0.585;
+    double yy;
+
+    const double C1 =1.22;
+    const double C2 =0.4288;
+    const double C3 =-1.651;
+    const double miu = 0.585;
+    const double del=1.05;
+    const double qRg = q*sqrt(Rgsquare(q,L,b));
+    const double x = pow(qRg, -1.0/miu);
+
 
     //calculating the derivative to decide on the corection (cutoff) term?
     // I have modified this from WRs original code
-
-    if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) {
-        C_star2 = 0.0;
-    } else {
-        C_star2 = 1.0;
-    }
-
-    Rg2 = Rgsquare(q,L,b);
-
-    yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) +
-        C_star2*w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) +
-        C2*pow((q*sqrt(Rg2)),(-2.0/miu)) + C3*pow((q*sqrt(Rg2)),(-3.0/miu)));
+    const double qdel = (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q);
+    const double C_star2 =(qdel >= 0.0) ? 0.0 : 1.0;
+
+    yy = (1.0 - w_WR(qRg))*Sdebye(q,L,b) +
+         C_star2*w_WR(qRg)*
+         x*(C1 + x*(C2 + x*C3));
 
     return (yy);
@@ -382 +366 @@
 double Sk_WR(double q, double L, double b)
 {
-  //
-  double p1,p2,p1short,p2short,q0;
-  double C,ans,q0short,Sexvmodify,pi;
-
-  pi = 4.0*atan(1.0);
-
-  p1 = 4.12;
-  p2 = 4.42;
-  p1short = 5.36;
-  p2short = 5.62;
-  q0 = 3.1;
-
-  q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0);
-
-
-  if(L/b > 10.0) {
-    C = 3.06/pow((L/b),0.44);
-  } else {
-    C = 1.0;
-  }
-  //
-
-  //
-  if( L > 4*b ) { // Longer Chains
-    if (q*b <= 3.1) {   //Modified by Yun on Oct. 15,
-
-      Sexvmodify = Sexvnew(q, L, b);
-
-      ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0*u_WR(q,L,b)) -
+    //
+    const double p1 = 4.12;
+    const double p2 = 4.42;
+    const double p1short = 5.36;
+    const double p2short = 5.62;
+    const double q0 = 3.1;
+
+    double q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0);
+    double Sexvmodify, ans;
+
+    const double C =(L/b > 10.0) ? 3.06/pow((L/b),0.44) : 1.0;
+
+    if( L > 4*b ) { // Longer Chains
+       if (q*b <= 3.1) {   //Modified by Yun on Oct. 15,
+         Sexvmodify = Sexvnew(q, L, b);
+         ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0*u_WR(q,L,b)) -
             (11.0/15.0 + 7.0/(15.0*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L);
 
-    } else { //q(i)*b > 3.1
-      ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) +
-            a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + pi/(q*L);
+       } else { //q(i)*b > 3.1
+         ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) +
+               a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + M_PI/(q*L);
+       }
+    } else { //L <= 4*b Shorter Chains
+       if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) {
+         if (q*b<=0.01) {
+            ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0;
+         } else {
+            ans = Sdebye1(q,L,b);
+         }
+       } else {  //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3)
+         ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p1short)) +
+               a2short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p2short)) +
+               M_PI/(q*L);
+       }
     }
-  } else { //L <= 4*b Shorter Chains
-    if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) {
-      if (q*b<=0.01) {
-        ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0;
-      } else {
-        ans = Sdebye1(q,L,b);
-      }
-    } else {  //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3)
-      ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p1short)) +
-            a2short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p2short)) +
-            pi/(q*L);
-    }
-  }
 
   return(ans);