Index: sasmodels/models/bcc_paracrystal.c
===================================================================
--- sasmodels/models/bcc_paracrystal.c	(revision 0b717c523aa8eedc5d395057a1a2919398332ca7)
+++ sasmodels/models/bcc_paracrystal.c	(revision 4962519cc1910422a467e4f643dba223623e0e09)
@@ -33,5 +33,5 @@
 
 	const double temp10 = exp((-1.0/8.0)*temp1*(temp7*temp7 + temp8*temp8 + temp9*temp9));
-	result = pow(1.0-(temp10*temp10),3)*temp6
+	result = cube(1.0 - (temp10*temp10))*temp6
 	    / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10)
 	      * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10)
Index: sasmodels/models/core_shell_sphere.py
===================================================================
--- sasmodels/models/core_shell_sphere.py	(revision 9a4811ac2bdb93459d640d111cba9f1016660790)
+++ sasmodels/models/core_shell_sphere.py	(revision 4962519cc1910422a467e4f643dba223623e0e09)
@@ -95,6 +95,6 @@
     """
     return (1, 1)
-    whole = 4.0 * pi / 3.0 * pow((radius + thickness), 3)
-    core = 4.0 * pi / 3.0 * radius * radius * radius
+    whole = 4.0/3.0 * pi * (radius + thickness)**3
+    core = 4.0/3.0 * pi * radius**3
     return whole, whole - core
 
Index: sasmodels/models/dab.py
===================================================================
--- sasmodels/models/dab.py	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/dab.py	(revision 4962519cc1910422a467e4f643dba223623e0e09)
@@ -59,6 +59,6 @@
 
 Iq = """
-    double numerator   = pow(cor_length, 3);
-    double denominator = pow(1 + pow(q*cor_length,2), 2);
+    double numerator   = cube(cor_length);
+    double denominator = square(1 + square(q*cor_length));
     
     return numerator / denominator ;
Index: sasmodels/models/fcc_paracrystal.c
===================================================================
--- sasmodels/models/fcc_paracrystal.c	(revision 0b717c523aa8eedc5d395057a1a2919398332ca7)
+++ sasmodels/models/fcc_paracrystal.c	(revision 4962519cc1910422a467e4f643dba223623e0e09)
@@ -32,5 +32,5 @@
 
 	const double temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9)));
-	result = pow((1.0-(temp10*temp10)),3)*temp6
+	result = cube(1.0-(temp10*temp10))*temp6
 	    / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10)
 	      * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10)
Index: sasmodels/models/flexible_cylinder_elliptical.c
===================================================================
--- sasmodels/models/flexible_cylinder_elliptical.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
+++ sasmodels/models/flexible_cylinder_elliptical.c	(revision 4962519cc1910422a467e4f643dba223623e0e09)
@@ -21,8 +21,7 @@
         double sn, cn;
         SINCOS(zi, sn, cn);
-        double arg = q*sqrt(a*a*sn*sn+b*b*cn*cn);
-        double yyy = pow((double)sas_J1c(arg),2);
-        yyy *= Gauss76Wt[i];
-        summ += yyy;
+        double arg = q*sqrt(a*a*sn*sn + b*b*cn*cn);
+        double yyy = sas_J1c(arg);
+        summ += Gauss76Wt[i] * yyy * yyy;
     }
 
Index: sasmodels/models/fractal_core_shell.py
===================================================================
--- sasmodels/models/fractal_core_shell.py	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/fractal_core_shell.py	(revision 4962519cc1910422a467e4f643dba223623e0e09)
@@ -100,6 +100,6 @@
         @param thickness: shell thickness
     """
-    whole = 4.0 * pi / 3.0 * pow((radius + thickness), 3)
-    core = 4.0 * pi / 3.0 * radius * radius * radius
+    whole = 4.0/3.0 * pi * (radius + thickness)**3
+    core = 4.0/3.0 * pi * radius**3
     return whole, whole-core
 
Index: sasmodels/models/hayter_msa.c
===================================================================
--- sasmodels/models/hayter_msa.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
+++ sasmodels/models/hayter_msa.c	(revision 4962519cc1910422a467e4f643dba223623e0e09)
@@ -48,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/(M_PI*Perm*SIdiam*pow((2.0+Kappa*SIdiam),2));
+	gMSAWave[5]=Beta*charge*charge/(M_PI*Perm*SIdiam*square(2.0+Kappa*SIdiam));
 	
 	//         Finally set up dimensionless parameters 
Index: sasmodels/models/lamellar_stack_paracrystal.c
===================================================================
--- sasmodels/models/lamellar_stack_paracrystal.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
+++ sasmodels/models/lamellar_stack_paracrystal.c	(revision 4962519cc1910422a467e4f643dba223623e0e09)
@@ -67,5 +67,5 @@
 	double Snq;
 
-	Snq = an/( (double)Nlayers*pow((1.0+ww*ww-2.0*ww*cos(qval*davg)),2) );
+	Snq = an/( (double)Nlayers*square(1.0+ww*ww-2.0*ww*cos(qval*davg)) );
 	
 	return(Snq);
Index: sasmodels/models/lib/wrc_cyl.c
===================================================================
--- sasmodels/models/lib/wrc_cyl.c	(revision ba32cddebeff04b66d828ca58b1843244d4425fb)
+++ sasmodels/models/lib/wrc_cyl.c	(revision 4962519cc1910422a467e4f643dba223623e0e09)
@@ -14,6 +14,5 @@
     //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) );
+    return pow(1.0+square(x/3.12)+cube(x/8.67), 0.176/3.0);
 }
 
@@ -23,6 +22,5 @@
 {
     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)));
+    return (L*b/6.0) * (1.0 - r*(1.5 + r*(1.5 + r*0.75*expm1(-2.0/r))));
 }
 
@@ -41,5 +39,5 @@
 }
 
-static inline double
+static double
 sech_WR(double x)
 {
@@ -51,5 +49,4 @@
 {
     double C;
-    const double onehalf = 1.0/2.0;
 
     if( L/b > 10.0) {
@@ -86,5 +83,5 @@
 
     const double t2 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) +
-         Rg02/b2))*((1.0 + onehalf*(((-1.0) -
+         Rg02/b2))*((1.0 + 0.5*(((-1.0) -
          tanh((-C4 + Rgb/C5)))))));
 
@@ -112,9 +109,9 @@
 
     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) -
+         (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))*((1.0 + 0.5*(((-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) +
+         Rg02/b2))*((1.0 + 0.5*(((-1.0) - tanh(((-C4) +
          Rgb)/C5))))));
 
@@ -131,5 +128,5 @@
 
     const double t14 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) +
-          Rg02/b2))*((1.0 + onehalf*(((-1.0) - tanh(((-C4) +
+          Rg02/b2))*((1.0 + 0.5*(((-1.0) - tanh(((-C4) +
           Rgb)/C5))))));
 
@@ -140,9 +137,9 @@
 
     double yy = (pow(q0,p1)*(((-((b*M_PI)/(L*q0))) +t1/L +t2/(q04*Rg22) +
-        onehalf*t3*t4)) + (t5*((pow(q0,(p1 - p2))*
+        0.5*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)) -
+        t8/(C5*q04*Rg22) + t9/(q04*Rg22) -t10/(q05*Rg22) + 0.5*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) +
+        t14/(q04*Rg22) + 0.5*t15*((1.0 + tanh(((-C4) +
         Rgb)/C5)))))))))));
 
@@ -154,5 +151,4 @@
 {
     double C;
-    const double onehalf = 1.0/2.0;
 
     if( L/b > 10.0) {
@@ -201,9 +197,9 @@
     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) +
+         ((1.0 + 0.5*(((-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) +
+         Rg02/b2))*((1.0 + 0.5*(((-1) - tanh(((-C4) +
          Rgb)/C5))))))/(q05*Rg22);
 
@@ -219,11 +215,11 @@
 
     const double t10 = (2.0*b4*(((-1) + pow((double)M_E,(-(Rg02/b2))) +
-          Rg02/b2))*((1.0 + onehalf*(((-1) - tanh(((-C4) +
+          Rg02/b2))*((1.0 + 0.5*(((-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))*
+         t2 + t3 - t4 + t5 - t6 + 0.5*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))*
          (((-((b*M_PI)/(L*q0))) + t9 + t10 +
-         onehalf*((C3*pow(((Rgb)),((-3.0)/miu)) +
+         0.5*((C3*pow(((Rgb)),((-3.0)/miu)) +
          C2*pow(((Rgb)),((-2.0)/miu)) +
          C1*pow(((Rgb)),((-1.0)/miu))))*
Index: sasmodels/models/linear_pearls.c
===================================================================
--- sasmodels/models/linear_pearls.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
+++ sasmodels/models/linear_pearls.c	(revision 4962519cc1910422a467e4f643dba223623e0e09)
@@ -43,14 +43,6 @@
     double separation = edge_sep + 2.0 * radius;
 
-    double x=q*radius;
-
-    // Try Taylor on x*xos(x)
-	// double out_cos = x - pow(x,3)/2 + pow(x,5)/24 - pow(x,7)/720 + pow(x,9)/40320;
-    // psi -= x*out_cos;
-
     //sine functions of a pearl
-    double psi = sin(q * radius);
-    psi -= x * cos(x);
-    psi /= pow((q * radius), 3.0);
+    double psi = sph_j1c(q * radius);
 
     // N pearls contribution
@@ -61,7 +53,5 @@
     }
     // form factor for num_pearls
-    double form_factor = n_contrib;
-    form_factor *= pow((m_s*psi*3.0), 2.0);
-    form_factor /= (tot_vol * 1.0e4);
+    double form_factor = 1.0e-4 * n_contrib * square(m_s*psi) / tot_vol;
 
     return form_factor;
Index: sasmodels/models/linear_pearls.py
===================================================================
--- sasmodels/models/linear_pearls.py	(revision 40a87fa5096610e24c3ead10299f609e99f093f2)
+++ sasmodels/models/linear_pearls.py	(revision 4962519cc1910422a467e4f643dba223623e0e09)
@@ -63,5 +63,5 @@
 single = False
 
-source = ["linear_pearls.c"]
+source = ["lib/sph_j1c.c", "linear_pearls.c"]
 
 demo = dict(scale=1.0, background=0.0,
Index: sasmodels/models/pearl_necklace.py
===================================================================
--- sasmodels/models/pearl_necklace.py	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/pearl_necklace.py	(revision 4962519cc1910422a467e4f643dba223623e0e09)
@@ -112,5 +112,5 @@
     """
     tot_vol = volume(radius, edge_sep, thick_string, num_pearls)
-    rad_out = pow((3.0*tot_vol/4.0/pi), 0.33333)
+    rad_out = (tot_vol/(4.0/3.0*pi)) ** (1./3.)
     return rad_out
 
Index: sasmodels/models/porod.py
===================================================================
--- sasmodels/models/porod.py	(revision 40a87fa5096610e24c3ead10299f609e99f093f2)
+++ sasmodels/models/porod.py	(revision 4962519cc1910422a467e4f643dba223623e0e09)
@@ -41,5 +41,5 @@
     """
     with errstate(divide='ignore'):
-        return power(q, -4)
+        return q**-4
 
 Iq.vectorized = True  # Iq accepts an array of q values
Index: sasmodels/models/sc_paracrystal.c
===================================================================
--- sasmodels/models/sc_paracrystal.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
+++ sasmodels/models/sc_paracrystal.c	(revision 4962519cc1910422a467e4f643dba223623e0e09)
@@ -49,11 +49,10 @@
 	double da = d_factor*dnn;
 	double temp1 = qq*qq*da*da;
-	double temp2 = pow( 1.0-exp(-1.0*temp1) ,3);
+	double temp2 = cube(-expm1(-temp1));
 	double temp3 = qq*dnn;
 	double temp4 = 2.0*exp(-0.5*temp1);
 	double temp5 = exp(-1.0*temp1);
 
-	double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5);
-	integrand /= M_PI_2;
+	double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5)/M_PI_2;
 
 	return(integrand);
Index: sasmodels/models/teubner_strey.py
===================================================================
--- sasmodels/models/teubner_strey.py	(revision b3f2a24980756e0d0b1e7b07be74b8df6c8674c5)
+++ sasmodels/models/teubner_strey.py	(revision 4962519cc1910422a467e4f643dba223623e0e09)
@@ -71,5 +71,5 @@
 
 import numpy as np
-from numpy import inf,power,pi
+from numpy import inf, pi
 
 name = "teubner_strey"
@@ -93,7 +93,7 @@
     drho2 = (sld-sld_solvent)*(sld-sld_solvent)
     k = 2.0*pi*xi/d
-    a2 = power(1.0+power(k,2.0),2.0)
-    c1 = -2.0*xi*xi*power(k,2.0)+2*xi*xi
-    c2 = power(xi,4.0)
+    a2 = (1.0 + k**2)**2
+    c1 = 2.0*xi**2 * (1.0 - k**2)
+    c2 = xi**4
     prefactor = 8.0*pi*volfraction*(1.0-volfraction)*drho2*c2/xi
     #k2 = (2.0*pi/d)*(2.0*pi/d)
