Index: sasmodels/models/barbell.c
===================================================================
--- sasmodels/models/barbell.c	(revision 11ca2ab20045d2ca8709faf56c10caaeaaa02734)
+++ sasmodels/models/barbell.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -39,5 +39,5 @@
     // translate dx in [-1,1] to dx in [lower,upper]
     const double integral = total*zm;
-    const double bell_fq = 2*M_PI*cube(radius_bell)*integral;
+    const double bell_fq = 2.0*M_PI*cube(radius_bell)*integral;
     return bell_fq;
 }
Index: sasmodels/models/binary_hard_sphere.c
===================================================================
--- sasmodels/models/binary_hard_sphere.c	(revision 87bc7073c5a3916988bead7a177b1e03724d5daf)
+++ sasmodels/models/binary_hard_sphere.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -49,6 +49,6 @@
     // /* do form factor calculations  */
     
-    v1 = 4.0*M_PI/3.0*r1*r1*r1;
-    v2 = 4.0*M_PI/3.0*r2*r2*r2;
+    v1 = M_4PI_3*r1*r1*r1;
+    v2 = M_4PI_3*r2*r2*r2;
     
     n1 = phi1/v1;
@@ -70,6 +70,6 @@
     sc1 = sph_j1c(qr1);
     sc2 = sph_j1c(qr2);
-    b1 = r1*r1*r1*(rho1-rhos)*4.0/3.0*M_PI*sc1;
-    b2 = r2*r2*r2*(rho2-rhos)*4.0/3.0*M_PI*sc2;
+    b1 = r1*r1*r1*(rho1-rhos)*M_4PI_3*sc1;
+    b2 = r2*r2*r2*(rho2-rhos)*M_4PI_3*sc2;
     inten = n1*b1*b1*psf11;
     inten += sqrt(n1*n2)*2.0*b1*b2*psf12;
Index: sasmodels/models/capped_cylinder.c
===================================================================
--- sasmodels/models/capped_cylinder.c	(revision 5bddd89a3ed5b121f80ca5d5663daeb28ce822d8)
+++ sasmodels/models/capped_cylinder.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -45,5 +45,5 @@
     // translate dx in [-1,1] to dx in [lower,upper]
     const double integral = total*zm;
-    const double cap_Fq = 2*M_PI*cube(radius_cap)*integral;
+    const double cap_Fq = 2.0*M_PI*cube(radius_cap)*integral;
     return cap_Fq;
 }
@@ -56,5 +56,5 @@
     const double bj = sas_J1c(q*radius*sin_alpha);
     const double si = sinc(q*half_length*cos_alpha);
-    const double cyl_Fq = 2.*M_PI*radius*radius*half_length*bj*si;
+    const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si;
     const double Aq = cap_Fq + cyl_Fq;
     return Aq;
Index: sasmodels/models/core_shell_parallelepiped.c
===================================================================
--- sasmodels/models/core_shell_parallelepiped.c	(revision 22221346d45e39a667fe231ce7ca5ecfc79b0ac5)
+++ sasmodels/models/core_shell_parallelepiped.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -89,8 +89,8 @@
             double tb = (a_scaled+2.0*thick_rim_b)/length_b;
     
-	        double arg1 = (0.5*mudum*a_scaled) * sin(0.5*M_PI*uu);
-	        double arg2 = (0.5*mudum) * cos(0.5*M_PI*uu);
-	        double arg3=  (0.5*mudum*ta) * sin(0.5*M_PI*uu);
-	        double arg4=  (0.5*mudum*tb) * cos(0.5*M_PI*uu);
+	        double arg1 = (0.5*mudum*a_scaled) * sin(M_PI_2*uu);
+	        double arg2 = (0.5*mudum) * cos(M_PI_2*uu);
+	        double arg3=  (0.5*mudum*ta) * sin(M_PI_2*uu);
+	        double arg4=  (0.5*mudum*tb) * cos(M_PI_2*uu);
 
 	        if(arg1==0.0){
Index: sasmodels/models/core_shell_sphere.c
===================================================================
--- sasmodels/models/core_shell_sphere.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
+++ sasmodels/models/core_shell_sphere.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -16,4 +16,4 @@
 double form_volume(double radius, double thickness)
 {
-    return 4.0 * M_PI / 3.0 * pow((radius + thickness), 3);
+    return M_4PI_3 * cube(radius + thickness);
 }
Index: sasmodels/models/flexible_cylinder_elliptical.c
===================================================================
--- sasmodels/models/flexible_cylinder_elliptical.c	(revision 87bc7073c5a3916988bead7a177b1e03724d5daf)
+++ sasmodels/models/flexible_cylinder_elliptical.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -18,5 +18,5 @@
 
     for(int i=0;i<N_POINTS_76;i++) {
-        double zi = ( Gauss76Z[i] + 1.0 )*M_PI/4.0;
+        double zi = ( Gauss76Z[i] + 1.0 )*M_PI_4;
         double sn, cn;
         SINCOS(zi, sn, cn);
Index: sasmodels/models/fractal_core_shell.c
===================================================================
--- sasmodels/models/fractal_core_shell.c	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/fractal_core_shell.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -13,5 +13,5 @@
 double form_volume(double radius, double thickness)
 {
-    return 4.0 * M_PI / 3.0 * pow((radius + thickness), 3);
+    return M_4PI_3 * cube(radius + thickness);
 }
 
Index: sasmodels/models/fuzzy_sphere.py
===================================================================
--- sasmodels/models/fuzzy_sphere.py	(revision 9a4811ac2bdb93459d640d111cba9f1016660790)
+++ sasmodels/models/fuzzy_sphere.py	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -86,5 +86,5 @@
 # This should perhaps be volume normalized?
 form_volume = """
-    return 1.333333333333333*M_PI*radius*radius*radius;
+    return M_4PI_3*cube(radius);
     """
 
Index: sasmodels/models/hayter_msa.c
===================================================================
--- sasmodels/models/hayter_msa.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
+++ sasmodels/models/hayter_msa.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -23,9 +23,7 @@
 	double SIdiam, diam, Kappa, cs, IonSt;
 	double  Perm, Beta;
-	double pi, charge;
+	double charge;
 	int ierr;
 	
-	pi = M_PI;
-
 	diam=2*radius_effective;		//in A
 
@@ -38,5 +36,5 @@
 	charge=zz*Elcharge;		//in Coulomb (C)
 	SIdiam = diam*1.0E-10;		//in m
-	Vp=4.0*pi/3.0*(SIdiam/2.0)*(SIdiam/2.0)*(SIdiam/2.0);	//in m^3
+	Vp=M_4PI_3*cube(SIdiam/2.0);	//in m^3
 	cs=csalt*6.022E23*1.0E3;	//# salt molecules/m^3
 	
@@ -50,5 +48,5 @@
 	Kappa=sqrt(2*Beta*IonSt/Perm);     //Kappa calc from Ionic strength
 									   //	Kappa=2/SIdiam					// Use to compare with HP paper
-	gMSAWave[5]=Beta*charge*charge/(pi*Perm*SIdiam*pow((2.0+Kappa*SIdiam),2));
+	gMSAWave[5]=Beta*charge*charge/(M_PI*Perm*SIdiam*pow((2.0+Kappa*SIdiam),2));
 	
 	//         Finally set up dimensionless parameters 
Index: sasmodels/models/hollow_rectangular_prism.c
===================================================================
--- sasmodels/models/hollow_rectangular_prism.c	(revision 87bc7073c5a3916988bead7a177b1e03724d5daf)
+++ sasmodels/models/hollow_rectangular_prism.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -34,7 +34,7 @@
    //Integration limits to use in Gaussian quadrature
     double v1a = 0.0;
-    double v1b = 0.5 * M_PI;  //theta integration limits
+    double v1b = M_PI_2;  //theta integration limits
     double v2a = 0.0;
-    double v2b = 0.5 * M_PI;  //phi integration limits
+    double v2b = M_PI_2;  //phi integration limits
     
     //Order of integration
@@ -88,5 +88,5 @@
     // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization)
     // The factor 2 is due to the different theta integration limit (pi/2 instead of pi)
-    answer *= (2.0/M_PI);
+    answer /= M_PI_2;
 
     // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization.
Index: sasmodels/models/hollow_rectangular_prism_thin_walls.c
===================================================================
--- sasmodels/models/hollow_rectangular_prism_thin_walls.c	(revision 87bc7073c5a3916988bead7a177b1e03724d5daf)
+++ sasmodels/models/hollow_rectangular_prism_thin_walls.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -26,7 +26,7 @@
    //Integration limits to use in Gaussian quadrature
     double v1a = 0.0;
-    double v1b = 0.5 * M_PI;  //theta integration limits
+    double v1b = M_PI_2;  //theta integration limits
     double v2a = 0.0;
-    double v2b = 0.5 * M_PI;  //phi integration limits
+    double v2b = M_PI_2;  //phi integration limits
     
     //Order of integration
@@ -71,5 +71,5 @@
     // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization)
     // The factor 2 is due to the different theta integration limit (pi/2 instead of pi)
-    answer *= (2.0/M_PI);
+    answer /= M_PI_2;
 
     // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization.
Index: sasmodels/models/lib/core_shell.c
===================================================================
--- sasmodels/models/lib/core_shell.c	(revision 7d4b2ae2af901813f3ff60e60aa0eae98f4dfd5e)
+++ sasmodels/models/lib/core_shell.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -17,5 +17,5 @@
     const double core_contrast = core_sld - shell_sld;
     const double core_bes = sph_j1c(core_qr);
-    const double core_volume = 4.0 * M_PI / 3.0 * radius * radius * radius;
+    const double core_volume = M_4PI_3 * cube(radius);
     double f = core_volume * core_bes * core_contrast;
 
@@ -24,5 +24,5 @@
     const double shell_contrast = shell_sld - solvent_sld;
     const double shell_bes = sph_j1c(shell_qr);
-    const double shell_volume = 4.0 * M_PI / 3.0 * pow((radius + thickness), 3);
+    const double shell_volume = M_4PI_3 * cube(radius + thickness);
     f += shell_volume * shell_bes * shell_contrast;
     return f * f * 1.0e-4;
Index: sasmodels/models/lib/gfn.c
===================================================================
--- sasmodels/models/lib/gfn.c	(revision 177c1a138481a9833b4929af687e86583b5c16d5)
+++ sasmodels/models/lib/gfn.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -12,5 +12,4 @@
 {
     // local variables
-    const double pi43=4.0/3.0*M_PI;
     const double aa = crmaj;
     const double bb = crmin;
@@ -20,10 +19,10 @@
     //const double siq = (uq == 0.0 ? 1.0 : 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq);
     const double siq = sph_j1c(uq);
-    const double vc = pi43*aa*aa*bb;
+    const double vc = M_4PI_3*aa*aa*bb;
     const double gfnc = siq*vc*delpc;
 
     const double ut2 = (trmin*trmin*xx*xx + trmaj*trmaj*(1.0-xx*xx));
     const double ut= sqrt(ut2)*qq;
-    const double vt = pi43*trmaj*trmaj*trmin;
+    const double vt = M_4PI_3*trmaj*trmaj*trmin;
     //const double sit = (ut == 0.0 ? 1.0 : 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut);
     const double sit = sph_j1c(ut);
@@ -33,4 +32,4 @@
     const double result = tgfn*tgfn;
 
-    return (result);
+    return result;
 }
Index: sasmodels/models/linear_pearls.c
===================================================================
--- sasmodels/models/linear_pearls.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
+++ sasmodels/models/linear_pearls.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -19,5 +19,5 @@
 {
     // Pearl volume
-    double pearl_vol = 4.0 /3.0 * M_PI * pow(radius, 3.0);
+    double pearl_vol = M_4PI_3 * cube(radius);
     // Return total volume
     return num_pearls * pearl_vol;;
@@ -35,5 +35,5 @@
     double contrast_pearl = pearl_sld - solvent_sld;
     //each volume
-    double pearl_vol = 4.0 /3.0 * M_PI * pow(radius, 3.0);
+    double pearl_vol = M_4PI_3 * cube(radius);
     //total volume
     double tot_vol = num_pearls * pearl_vol;
Index: sasmodels/models/mass_fractal.c
===================================================================
--- sasmodels/models/mass_fractal.c	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/mass_fractal.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -34,5 +34,5 @@
 double form_volume(double radius){
 
-    return 1.333333333333333*M_PI*radius*radius*radius;
+    return M_4PI_3*radius*radius*radius;
 }
 
Index: sasmodels/models/mass_surface_fractal.c
===================================================================
--- sasmodels/models/mass_surface_fractal.c	(revision 30fbe2eda4d6f97acc840bd5e605de3829caf3c7)
+++ sasmodels/models/mass_surface_fractal.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -37,5 +37,5 @@
 double form_volume(double radius){
 
-    return 1.333333333333333*M_PI*radius*radius*radius;
+    return M_4PI_3*radius*radius*radius;
 }
 
Index: sasmodels/models/multilayer_vesicle.c
===================================================================
--- sasmodels/models/multilayer_vesicle.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
+++ sasmodels/models/multilayer_vesicle.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -19,5 +19,5 @@
 
         // layer 1
-        voli = 4.0*M_PI/3.0*ri*ri*ri;
+        voli = M_4PI_3*ri*ri*ri;
         fval += voli*sldi*sph_j1c(ri*q);
 
@@ -25,5 +25,5 @@
 
         // layer 2
-        voli = 4.0*M_PI/3.0*ri*ri*ri;
+        voli = M_4PI_3*ri*ri*ri;
         fval -= voli*sldi*sph_j1c(ri*q);
 
Index: sasmodels/models/parallelepiped.c
===================================================================
--- sasmodels/models/parallelepiped.c	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/parallelepiped.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -71,6 +71,6 @@
             double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
             double mudum = mu * sqrt(1.0-sigma*sigma);
-	        double arg1 = 0.5 * mudum * cos(0.5*M_PI*uu);
-	        double arg2 = 0.5 * mudum * a_scaled * sin(0.5*M_PI*uu);
+	        double arg1 = 0.5 * mudum * cos(M_PI_2*uu);
+	        double arg2 = 0.5 * mudum * a_scaled * sin(M_PI_2*uu);
             if(arg1==0.0) {
 	        tmp1 = 1.0;
Index: sasmodels/models/pearl_necklace.c
===================================================================
--- sasmodels/models/pearl_necklace.c	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/pearl_necklace.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -23,7 +23,4 @@
 	double num_strings = num_pearls - 1.0;
 	
-	//Pi
-	double pi = 4.0*atan(1.0);
-	
 	// center to center distance between the neighboring pearls
 	double A_s = edge_sep + 2.0 * radius;
@@ -36,6 +33,6 @@
 	
 	// each volume
-	double string_vol = edge_sep * pi * thick_string * thick_string / 4.0;
-	double pearl_vol = 4.0 / 3.0 * pi * radius * radius * radius;
+	double string_vol = edge_sep * M_PI * thick_string * thick_string / 4.0;
+	double pearl_vol = M_4PI_3 * radius * radius * radius;
 
 	//total volume
@@ -110,9 +107,8 @@
 	double total_vol;
 
-	double pi = 4.0*atan(1.0);
 	double number_of_strings = num_pearls - 1.0;
 	
-	double string_vol = edge_sep * pi * thick_string * thick_string / 4.0;
-	double pearl_vol = 4.0 / 3.0 * pi * radius * radius * radius;
+	double string_vol = edge_sep * M_PI * thick_string * thick_string / 4.0;
+	double pearl_vol = M_4PI_3 * radius * radius * radius;
 
 	total_vol = number_of_strings * string_vol;
Index: sasmodels/models/rectangular_prism.c
===================================================================
--- sasmodels/models/rectangular_prism.c	(revision 87bc7073c5a3916988bead7a177b1e03724d5daf)
+++ sasmodels/models/rectangular_prism.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -75,5 +75,5 @@
     // The factor 2 appears because the theta integral has been defined between 
     // 0 and pi/2, instead of 0 to pi.
-    answer *= (2.0/M_PI); //Form factor P(q)
+    answer /= M_PI_2; //Form factor P(q)
 
     // Multiply by contrast^2 and volume^2
Index: sasmodels/models/sc_paracrystal.c
===================================================================
--- sasmodels/models/sc_paracrystal.c	(revision 0b717c523aa8eedc5d395057a1a2919398332ca7)
+++ sasmodels/models/sc_paracrystal.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -55,5 +55,5 @@
 
 	double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5);
-	integrand *= 2.0/M_PI;
+	integrand /= M_PI_2;
 
 	return(integrand);
@@ -68,5 +68,5 @@
 {
 	const double va = 0.0;
-	const double vb = M_PI/2.0; //orientation average, outer integral
+	const double vb = M_PI_2; //orientation average, outer integral
 
     double summ=0.0;
Index: sasmodels/models/star_polymer.c
===================================================================
--- sasmodels/models/star_polymer.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
+++ sasmodels/models/star_polymer.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -6,11 +6,11 @@
 {
 
-    double u_2 = radius2 * pow(q,2);
+    double u_2 = radius2 * q * q;
     double v = u_2 * arms / (3.0 * arms - 2.0);
 
-    double term1 = v - 1.0 + exp(-v);
-    double term2 = ((arms - 1.0)/2.0)* pow((1.0 - exp(-v)),2.0);
+    double term1 = v + expm1(-v);
+    double term2 = ((arms - 1.0)/2.0) * square(expm1(-v));
 
-    return (2.0 * (term1 + term2)) / (arms * pow(v,2.0));
+    return (2.0 * (term1 + term2)) / (arms * v * v);
 
 }
Index: sasmodels/models/triaxial_ellipsoid.c
===================================================================
--- sasmodels/models/triaxial_ellipsoid.c	(revision 11ca2ab20045d2ca8709faf56c10caaeaaa02734)
+++ sasmodels/models/triaxial_ellipsoid.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -10,5 +10,5 @@
 double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar)
 {
-    return 1.333333333333333*M_PI*radius_equat_minor*radius_equat_major*radius_polar;
+    return M_4PI_3*radius_equat_minor*radius_equat_major*radius_polar;
 }
 
Index: sasmodels/models/vesicle.c
===================================================================
--- sasmodels/models/vesicle.c	(revision 2c74c11fc4544381ecf9949bd1618c67ff624c6d)
+++ sasmodels/models/vesicle.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -8,8 +8,5 @@
 {
     //note that for the vesicle model, the volume is ONLY the shell volume
-    double volume;
-    volume =4.*M_PI*(radius+thickness)*(radius+thickness)*(radius+thickness)/3;
-    volume -=4.*M_PI*radius*radius*radius/3.;
-    return volume;
+    return M_4PI_3*(cube(radius+thickness) - cube(radius));
 }
 
@@ -32,15 +29,15 @@
     // core first, then add in shell
     contrast = sld_solvent-sld;
-    vol = 4.0*M_PI/3.0*radius*radius*radius;
-    f = vol*sph_j1c(q*radius)*contrast;
+    vol = M_4PI_3*cube(radius);
+    f = vol * sph_j1c(q*radius) * contrast;
  
     //now the shell. No volume normalization as this is done by the caller
     contrast = sld-sld_solvent;
-    vol = 4.0*M_PI/3.0*(radius+thickness)*(radius+thickness)*(radius+thickness);
-    f += vol*sph_j1c(q*(radius+thickness))*contrast;
+    vol = M_4PI_3*cube(radius+thickness);
+    f += vol * sph_j1c(q*(radius+thickness)) * contrast;
 
     //rescale to [cm-1]. 
-    f2 = volfraction*f*f*1.0e-4;
+    f2 = volfraction * f*f*1.0e-4;
     
-    return(f2);
+    return f2;
 }
Index: sasmodels/models/vesicle.py
===================================================================
--- sasmodels/models/vesicle.py	(revision e77872ef7445f58a9eb6ad1eb0015b7b63d53fac)
+++ sasmodels/models/vesicle.py	(revision 3a48772f136eb2795903c63707474e7527b4a325)
@@ -116,6 +116,6 @@
     '''
 
-    whole = 4. * pi * (radius + thickness) ** 3. / 3.
-    core = 4. * pi * radius ** 3. / 3.
+    whole = 4./3. * pi * (radius + thickness)**3
+    core = 4./3. * pi * radius**3
     return whole, whole - core
 
